Source code :: basic_diagn
[Return]
[Download]#!/usr/bin/python2
import sys
sys.path.append('/home/michal/Desktop/Prace/Golem/web_framework')
from numpy import *
from pygolem.config import *
from pygolem.modules import *
from matplotlib.pyplot import *
from scipy.optimize import fmin, leastsq
import time
from scipy.sparse import spdiags, eye
from scipy.sparse.linalg import spsolve, splu
from scipy.signal import medfilt
from scipy.stats.mstats import mquantiles
print "OPRAVIT NACITANI !!!!"
#d = loadtxt('Nidatap.lvm')
#save('Nidatap', single(d))
data = load('Nidatap.npz')['data']
#data = loadtxt('Nidatap.lvm', usecols=[0,1,2,3,4])
tvec = data[:,0]
dt = tvec[1] - tvec[0] # time step
CdFieldTrigger = loadconst("Tcd_aktual") * 1e-6 + TororoidalMagneticFieldTrigger # [s]
BdFieldTrigger = loadconst('Tbd_aktual') * 1e-6 + TororoidalMagneticFieldTrigger # [s]
StFieldTrigger = loadconst('Tst_aktual') * 1e-6 + TororoidalMagneticFieldTrigger # [s]
BdVoltage = loadconst('Ubd_limit') # [V]
CdVoltage = loadconst('Ucd_limit') # [V]
StVoltage = loadconst('Ust_limit') # [V]
if BdVoltage == 0:
BdFieldTrigger = inf
if CdVoltage == 0:
CdFieldTrigger = inf
if StFieldTrigger == 0:
StFieldTrigger = inf
def getBtoroidal():
dBt = data[:,2] # mag field derivative
mBt = mean(dBt)
dBt -= mBt # remove bias
Bt = cumsum(dBt) * dt * Bt_calibration # integrate
BtMax = amax(Bt) # max
BtMean = mean(Bt) # mean
saveconst('BtMax', BtMax)
#saveconst('BtMean', BtMean)
savetxt_adv('Btoroidal_dp', tvec, Bt)
savetxt_adv('dBdt_toroidal_dp', tvec, dBt)
return Bt, BtMax, BtMean
def getUloop():
Uloop = data[:,1]
m_Uloop = mean(Uloop)
Uloop *= UloopCalibration
start = TororoidalMagneticFieldTrigger + 2e-3 # [s]
Uloop[tvec < start] = 0 # remove noise from tyristors
ind_wrong = abs(TororoidalMagneticFieldTrigger-tvec) < 1e-4 # [s]
Uloop[ind_wrong] = median(Uloop[ind_wrong] ) #remove noise from toroidal field trigger
UloopMax = mquantiles(Uloop,0.99)[0]
UloopMean = sum(cumsum(abs(Uloop))) / len(Uloop) ** 2 # used for detection of
saveconst('UloopMax', UloopMax)
#saveconst('UloopMean', UloopMean)
savetxt_adv('Uloop_dp', tvec, Uloop)
return Uloop, UloopMax, UloopMean
def getIrogowski():
dIrog = data[:,3]
ind_wrong = abs( CdFieldTrigger - tvec) < 1e-4 # [s]
ind_wrong |= abs(TororoidalMagneticFieldTrigger - tvec) < 1e-4 # [s]
ind_wrong |= abs(BdFieldTrigger - tvec) < 1e-4 # [s]
ind_wrong |= abs(StFieldTrigger - tvec) < 1e-4 # [s]
dIrog[ind_wrong] = 0 # cut off triggers
#plot(dIrog)
#show()
Ibias = mean(dIrog[tvec < max(amin([TororoidalMagneticFieldTrigger, BdFieldTrigger, CdFieldTrigger, StFieldTrigger])-1e-4 , 2e-3)] )
dIrog -= Ibias
Irog = cumsum(dIrog) * dt
Irog *= RogowskiCalibration
Irog[tvec < min(CdFieldTrigger, BdFieldTrigger) + 1e-3] = 0 # remove noise
dIdt_rogMax = mquantiles(diff(Irog),0.99)[0] # maximum
IrogMax = mean(Irog)
saveconst('dIdt_rogMax', dIdt_rogMax)
saveconst('IrogMax', IrogMax)
savetxt_adv('Irogowski_dp', tvec, Irog)
savetxt_adv('dIdt_rogowski_dp', tvec, dIrog)
#plot(Irog)
#show()
return Irog, dIdt_rogMax, IrogMax
def getIplasma(Uloop,Irog ):
ind = (tvec > CdFieldTrigger + 0.1e-3) & (tvec < CdFieldTrigger + 1.5e-3)
if all(~ind): # no intersect
ind = (tvec > BdFieldTrigger + 0.1e-3) & (tvec < BdFieldTrigger + 1.5e-3)
npix = len(Uloop)
diag0=spdiags(ones(npix), 0,npix,npix,format='dia')
diag1=spdiags(ones(npix), -1,npix,npix,format='dia')
def fit(param, result = False):
Ich = spsolve(diag0*(param[0]+param[1]) + diag1*-param[0],Uloop)
if result:
return Ich
return (Irog - Ich)[ind]
[a,b] = [1.66, 0.0103] # material constants
x = leastsq(fit, array([a,b]), xtol=1e-3)
if x[1] == 0: # failed convergence
x[0] = [a,b]
Ich = fit(x[0], result = True)
Ipl = Irog - Ich
savetxt_adv('Iplasma_dp', tvec, Ipl)
savetxt_adv('Ich_dp', tvec, Ich)
return Ipl, Ich
def getPhotod(PlasmaStart):
Photod = data[:,4]
bias = median(Photod[tvec < min(CdFieldTrigger, BdFieldTrigger, StFieldTrigger) -2e-4 ]) # remove noise from tyristors
Photod -= bias
Photod[tvec < nanmax([CdFieldTrigger, BdFieldTrigger, PlasmaStart]) + 2e-4 ] = 0
Photod *= sign(mean(Photod)) # invert signal
Photod[Photod < 0] = 0
savetxt_adv('Photod_dp', tvec, Photod)
return Photod
def getPhotodHalpha(PlasmaStart):
PhotodHalfa = data[:,4]
bias = median(PhotodHalfa[tvec < min(CdFieldTrigger, BdFieldTrigger, StFieldTrigger) -2e-4 ]) # remove noise from tyristors
PhotodHalfa -= bias
PhotodHalfa[tvec < nanmax([CdFieldTrigger, BdFieldTrigger, PlasmaStart]) + 2e-4 ] = 0
PhotodHalfa *= sign(mean(PhotodHalfa)) # invert signal
PhotodHalfa[PhotodHalfa < 0] = 0
savetxt_adv('PhotodHalpha_dp', tvec, PhotodHalfa)
return PhotodHalfa
def PlasmaDetect(Iplasma):
polar = sign(mean(Iplasma))
ReversedCD = (polar == -1)
saveconst("ReversedCD", ReversedCD)
N_steps = len(Iplasma)
dt = tvec[1] - tvec[0]
PlasmaStart=tvec[0]+1.5e-3
downsample = max(len(Iplasma)/1000, 1)
d = medfilt(Iplasma[::downsample], int((3e-4/dt/downsample)/2)*2+1)
d = medfilt( diff(d)/dt/downsample, int((3e-4/dt/downsample)/2)*2+1)
################
min_dI = 1e5 # sensitivity of plasma detection
####################
try:
PlasmaStart = tvec[::downsample][where(d > max(0.2*amax(d), min_dI))[0][0]]
PlasmaEnd = tvec[::downsample][where(d < min(0.2*amin(d), -min_dI))[0][-1]]
except:
PlasmaStart = nan
PlasmaEnd = nan
Plasma = int(~ isnan(PlasmaStart) and ~ isnan(PlasmaEnd))
PlasmaTimeLength = PlasmaStart - PlasmaEnd
saveconst("PlasmaStartAdvanced", PlasmaStart - 0.1e-3)
saveconst("PlasmaStart", PlasmaStart)
saveconst("PlasmaEnd", PlasmaEnd)
saveconst("PlasmaEndDelayed", PlasmaEnd * 1.1)
saveconst("Plasma", Plasma )
saveconst("PlasmaTimeLength", PlasmaTimeLength)
return Plasma, PlasmaStart, PlasmaEnd, PlasmaTimeLength
def getMeanBt(Btor,PlasmaStart, PlasmaEnd ):
MeanBt = mean(Btor[(tvec > PlasmaStart) & (tvec < PlasmaEnd)])
saveconst('BtMean', MeanBt)
saveconst('ReversedBt', int(MeanBt > 0))
return MeanBt
def getMeanPhotod(Photod, PlasmaStart, PlasmaEnd):
MeanPhotod = median(Photo[(tvec > PlasmaStart) & (tvec < PlasmaEnd)])
saveconst("MeanPhotod", MeanPhotod)
return MeanPhotod
def getTotalCharge(Ipla, PlasmaStart,PlasmaEnd ):
TotalCharge = (Ipla + abs(Ipla)) / 2 * dt
TotalCharge = cumsum(TotalCharge[(tvec > PlasmaStart) & (tvec < PlasmaEnd)])
TotalCharge = mean(TotalCharge)
saveconst('TotalCharge', TotalCharge)
return TotalCharge
def getMeanCurrent(Ipla, PlasmaStart,PlasmaEnd ):
MeanIpla = median(Ipla[(tvec > PlasmaStart) & (tvec < PlasmaEnd)])
saveconst("IplaMean", MeanIpla)
return MeanIpla
def getMeanUloop(Uloop, PlasmaStart,PlasmaEnd ):
MeanUloop = median(Uloop[(tvec > PlasmaStart) & (tvec < PlasmaEnd)])
saveconst("UloopMean", MeanUloop)
return MeanUloop
def getOhmicHeatingPower(MeanUloop,MeanPlasmaCurrent):
OhmicHeatingPower = abs(MeanUloop*MeanPlasmaCurrent)
saveconst("OhmicHeatingPower", OhmicHeatingPower)
return OhmicHeatingPower
def getQedge(MeanBt,MeanPlasmaCurrent):
Qedge = 90.3*abs(MeanBt/(MeanPlasmaCurrent))
saveconst("Qedge", Qedge)
return Qedge
def getElectronTemperature(MeanUloop, MeanPlasmaCurrent):
ElectronTemperature = 89.8* (abs( MeanPlasmaCurrent / MeanUloop )) ** (2/3) #Ref: PhD Brotankova eq. 3.21 p. 26
saveconst("ElectronTemperature", ElectronTemperature)
return ElectronTemperature
def getBreakDownVoltage(Uloop, Btor, Ipla, PlasmaStart, PlasmaEnd):
break_ind = argmax(Uloop[(tvec > PlasmaStart) & (tvec < (PlasmaStart + PlasmaEnd)/2)])
Umax = Uloop[break_ind]
Btime = tvec[break_ind]
Bbreak = Btor[break_ind]
Bipla = Btor[break_ind]
saveconst("BreakDownVoltage", Umax)
saveconst("BreakDownTime", Btime)
saveconst("BreakDownBt", Bbreak)
saveconst("BreakDownIp", Bipla)
return Umax, Btime, Bbreak, Bipla
def getStateEqElectronDensity(Aktual_PfeifferMerkaVakua):
StateEqElectronDensity = (Aktual_PfeifferMerkaVakua/1000)/(kB * RoomTemperature)
saveconst("StateEqElectronDensity", StateEqElectronDensity)
return StateEqElectronDensity
def getElectronConfinementTimeFirstApprox(MeanUloop,MeanPlasmaCurrent, StateEqElectronDensity, ElectronTemperature ):
ElectronConfinementTimeFirstApprox = 3/8* PlasmaVolume * (ElectronTemperature* eV * StateEqElectronDensity) / abs( MeanPlasmaCurrent * MeanUloop) #Ref: PhD Brotankova
saveconst("ElectronConfinementTimeFirstApprox", ElectronConfinementTimeFirstApprox)
return ElectronConfinementTimeFirstApprox
def getChamberResistance(Plasma):
if Plasma:
Ibd = loadconst("BreakDownIp")
Ubd = loadconst("BreakDownVoltage")
else:
Ubd = loadconst("UloopMax")
Ibd = loadconst("IrogowskiMax")
ChamberResistance = abs( Ubd / Ibd )
saveconst("ChamberResistance", ChamberResistance)
return ChamberResistance
def Failures(Plasma, UloopMax, dIdt_rogMax, MeanUloop, BtMax, MeanBt, PlasmaStart, PlasmaEnd ):
if abs(UloopMax) < 1:
Plasma = "Failure (UloopMax < 1V)"
if abs(dIdt_rogMax) < 0.1 :
Plasma = "Failure (dIdt_rogMax < 0.1V)"
if abs(MeanUloop) < 0.3 :
Plasma = "Failure(MeanUloop) < 0.3V)"
if abs(BtMax) < 0.01 :
Plasma = "Failure(BtMax) < 0.01T)"
if abs(MeanBt) < 0.01 :
Plasma = "Failure(MeanBt) < 0.01T)"
saveconst("Plasma", Plasma)
return Plasma
def ElectronTemperatureTime(Uloop, Iplasma, PlasmaStart, PlasmaEnd):
ElectronTempTime = 89.8*abs(Uloop*0.001/Iplasma)**(2/3)
medElectronTemp = medfilt(ElectronTempTime, 50)
maxT = mquantiles(medElectronTemp[(tvec > PlasmaStart) & (tvec < PlasmaEnd)],0.99)[0]
savetxt_adv("ElectronTempTime", tvec, ElectronTempTime)
savetxt_adv("medElectronTemp", tvec, medElectronTemp)
saveconst("maxElectronTemperatureT", maxT)
return ElectronTempTime, medElectronTemp, maxT
def getGreenwaldDensity(Ipla):
N_gw = Ipla/(pi*LimiterPosition**2)
savetxt_adv("GreenwaldDensityTime", tvec, N_gw)
return N_gw
def getOhmicHeatingPowerTime(Ipla, Uloop):
Power = Ipla * Uloop
savetxt_adv("OhmicHeatingPowerTime", tvec, Power)
return Power
def getQedgeTime(Btor,Ipla):
QedgeTime = 90.3*abs(Btor/Ipla)
savetxt_adv("QedgeTime",tvec, QedgeTime)
return QedgeTime
def getMagneticFlux(Uloop):
Flux = cumsum(Uloop+abs(Uloop))/2*dt
savetxt_adv("TotalMagneticFlux", tvec, Flux)
return Flux[Return]