#!/usr/bin/python2
# -*- coding: utf-8 -*-
""" CREATED: 7/2012
AUTHOR: MICHAL ODSTRČIL
"""
from numpy import *
from pygolem_lite.config import *
from pygolem_lite.modules import *
from pygolem_lite.utilities import *
from pygolem_lite import *
from matplotlib.pyplot import *
from scipy.optimize import leastsq
import time, datetime, re
from scipy.sparse import spdiags, eye
from scipy.sparse.linalg import spsolve, splu
from scipy.signal import medfilt, fftconvolve, convolve
from scipy.stats.mstats import mquantiles
print "------------basicdiagn --------"
CdFieldTrigger = loadconst("Tcd_aktual") * 1e-6 + TororoidalMagneticFieldTrigger # [s]
BdFieldTrigger = loadconst('Tbd_aktual') * 1e-6 + TororoidalMagneticFieldTrigger # [s]
StFieldTrigger = loadconst('Tst_aktual') * 1e-6 + TororoidalMagneticFieldTrigger # [s]
saveconst('Tb', TororoidalMagneticFieldTrigger) # save in SI units
saveconst('Tbd', BdFieldTrigger) # save in SI units
saveconst('Tcd', CdFieldTrigger) # save in SI units
saveconst('Tst', StFieldTrigger) # save in SI units
BVoltage = loadconst('Ub_limit') # [V]
CdVoltage = loadconst('Ucd_limit') # [V]
BdVoltage = loadconst('Ubd_limit') # [V]
StVoltage = loadconst('Ust_limit') # [V]
if BdVoltage < 10:
BdFieldTrigger = nan
if CdVoltage == 0:
CdFieldTrigger = nan
if StFieldTrigger < 10:
StFieldTrigger = nan
print "CdFieldTrigger", CdFieldTrigger, "BdFieldTrigger", BdFieldTrigger, "StFieldTrigger", StFieldTrigger
print "BVoltage", BVoltage, "CdVoltage", CdVoltage, 'BdVoltage', BdVoltage, "StVoltage", StVoltage
try:
# get time vector
S = Shot()
tvec, tmp = S['any', 'btor']
shot = S.shot_num
### data damaged by triggers !!!! do not change !!!
ind_wrong = zeros(len(tvec), dtype = bool)
###ind_wrong |= abs( CdFieldTrigger - tvec) < 1e-4 # [s]
ind_wrong |= abs(TororoidalMagneticFieldTrigger - tvec + 0.5e-4 ) < 5e-4 # [s]
ind_wrong |= abs(BdFieldTrigger - tvec + 0.5e-4) < 0.5e-4 # [s]
ind_wrong |= abs(StFieldTrigger - tvec + 0.5e-4) < 0.5e-4 # [s]
dt = mean(diff(tvec)) # time step
except Exception, e:
print "Basicdiagn failure: ", e
params = []
for i in str(S['wwwcomment']).split("="):
try:
if re.match( "\d+.*", i):
i = re.sub( "(\ *[\d\,\.]+).*", r"\1", i)
i = re.sub(",", ".", i)
params.append( float(i) )
print "parametr", float(i) , "loaded"
except Exception, e:
print "Error in loading data from comment", str(e)
for i in range(len(params)):
saveconst('param_'+str(i+1), params[i]) # save in SI units
print "===================params", params
def save_config():
names = [
'ToroidalMagneticFieldCoilInductance',
'TororoidalMagneticFieldCapacitor',
'BreakDownElectricFieldCapacitor',
'CurrentDriveElectricFieldCapacitor',
'StabilizationCapacitor',
'Zeff'
]
for i in names:
exec('saveconst("'+i+'", '+i+')')
def getDate():
s = cat( "date")
#nicely formated date
shot_date = datetime.date( 2000 + int(s[4:6]), int(s[2:4]), int(s[:2]) )
s = cat( "starttime")
shot_time = datetime.time( int(s[0:2]), int(s[2:4]), int(s[4:6]) )
save('shot_date', shot_date)
save('shot_time', shot_time)
def reverse_if_changed(signal, orientation_file, normal_orientation):
orientation = cat(orientation_file)
reverse = (orientation != normal_orientation)
signal = -signal if reverse else signal
return reverse, signal
def getBtoroidal():
tvec, dBt = S['any', 'btor'] # mag field derivative
#mBt = mean(dBt)
min_time = max(nanmin([TororoidalMagneticFieldTrigger, BdFieldTrigger, CdFieldTrigger, StFieldTrigger])-1e-4 , 2e-3)
mBt = mean(dBt[tvec < min_time] )
sBt = std(dBt[tvec < min_time] )
dBt -= mBt # remove bias
Bt = cumsum(dBt) * dt * Bt_calibration # integrate
ReversedB, Bt = reverse_if_changed(Bt, '../../../Bt_orientation', 'ACW')
saveconst("ReversedB", ReversedB)
BtMax = amax(Bt) # max
BtMean = mean(Bt) # mean
saveconst('BtMax', BtMax, fmt="%.3f" )
#saveconst('BtMean', BtMean)
save_adv('Btoroidal', tvec, Bt)
save_adv('dBdt_toroidal', tvec, dBt)
saveconst('noise_level', sBt)
return Bt, dBt, BtMax, BtMean
def getUloop():
tvec, Uloop = S['any', 'uloop']
#m_Uloop = mean(Uloop)
Uloop *= UloopCalibration
ind_zero = (tvec < TororoidalMagneticFieldTrigger + 1e-4)
Uloop[ind_zero] = 0 # remove noise from tyristors
ReversedCD, Uloop = reverse_if_changed(Uloop, '../../../CD_orientation', 'CW')
saveconst("ReversedCD", ReversedCD)
#ind_wrong = abs(TororoidalMagneticFieldTrigger-tvec) < 1e-4 # [s]
#Uloop[ind_wrong] = median(Uloop[ind_wrong] ) #remove noise from toroidal field trigger
#Uloop[ind_wrong] = 0 # medfilt(Uloop, 51)[ind_wrong] # remove signal damaged by tyristors
#plot(Uloop)
#show()
#ind_bias = [(tvec > TororoidalMagneticFieldTrigger) & ( tvec < nanmin([CdFieldTrigger, StFieldTrigger])) | (abs(tvec - TororoidalMagneticFieldTrigger) < 1e-4 ) ]
#Uloop_bias = median(Uloop[ind_bias])
#Uloop -= Uloop_bias
#plot(tvec, Uloop)
#plot(tvec[ind_bias], Uloop[ind_bias])
#show()
#Uloop = medfilt(Uloop, 5) # odstranit ?? nicim ten osklivy sum v datech !!!
Uloop[ind_zero] = 0 # remove noise from tyristors
UloopMax = mquantiles(Uloop,0.99)[0]
#UloopMean = sum(cumsum(abs(Uloop))) / len(Uloop) ** 2 # used for detection of
UloopMean = mean(abs(Uloop))
#exit()
saveconst('UloopMax', UloopMax, fmt="%.3f" )
save_adv('Uloop', tvec, Uloop)
return Uloop, UloopMax, UloopMean, ReversedCD
def getIrogowski(ReversedCD):
tvec, dIrog = S['any', 'irog']
#dIrog[abs(tvec - TororoidalMagneticFieldTrigger) < 1e-4 ] = 0 # cut off triggers
if ReversedCD:
dIrog *= -1
Ibias = mean(dIrog[tvec < max(nanmin([TororoidalMagneticFieldTrigger, BdFieldTrigger, CdFieldTrigger, StFieldTrigger])-1e-4 , 2e-3)] )
dIrog[ind_wrong & (tvec < nanmin([CdFieldTrigger, 40e-3])) ]=0
dIrog -= Ibias
Irog = cumsum(dIrog) * dt
Irog *= RogowskiCalibration
#Irog[tvec < nanmin([CdFieldTrigger, BdFieldTrigger])] = 0 # remove noise
#plot(tvec, dIrog/amax(dIrog))
#plot(tvec, convolve(medfilt(dIrog, 51),ones(200)/200, mode='same')/amax(dIrog))
#plot(tvec, convolve(Irog,ones(200)/200 , mode='same'))
#show()
try:
# try to detect cd trigger from data
ker = ones(200)/200
I_start = tvec[where( (fftconvolve(dIrog, ker, mode='same')> 0.05) | (fftconvolve(Irog,ker, mode='same') > 5) )[0][0]]
except:
I_start = nan # no I start
print "I_start", I_start
Irog[tvec < I_start - 1e-4] = 0
dIdt_rogMax = mquantiles(abs(diff(Irog)),0.98)[0] # maximum
IrogMax = mquantiles(abs(Irog),0.98)[0]
saveconst('dIdt_rogMax', dIdt_rogMax, fmt="%.3f" )
saveconst('IrogowskiMax', IrogMax, fmt="%.3f" )
save_adv('Irogowski', tvec, Irog)
save_adv('dIdt_rogowski', tvec, dIrog)
#plot(medfilt(dIrog,51))
#plot(convolve(dIrog,ones(200)/200, mode='full'))
#show()
#exit()
return Irog, dIdt_rogMax, IrogMax, I_start
def getIplasma(Uloop,Irog,I_start, PlasmaStart = None, PlasmaEnd = None ):
#plot(Uloop)
#plot(Irog)
#show()
#exit()
#print Uloop
t0 = time.time()
x0 = array([ 1.2 , ChamberResistance])
time_lag = 1.5e-3 # [ms] expected time before breakdown
if PlasmaStart is not None:
ind = (tvec > I_start) & (tvec < min(PlasmaStart - 1e-4, I_start + time_lag) )
else:
ind = (tvec > I_start) & (tvec < I_start + time_lag)
def fit(param, result = False):
Ich = zeros(len(tvec))
search_range = ( range(ind[0], ind[-1]) if not result else range(1, len(tvec)) )
#search_range = range(1, len(tvec))
for i in search_range: # implicite differencial equation I_n = ( U_n + C1/dt *I_(n-1) ) / (C2 + C1/dt)
Ich[i]=(Uloop[i]+param[0]*Ich[i-1])/(param[1]+param[0]);
if result:
return Ich
return (Irog - Ich)[ind]
conv = 0
x = nan
#time_lag = 1.5e-3 # [ms] expected time before breakdown
#while time_lag > 0.5e-3:
#if PlasmaStart is not None:
#ind = (tvec > I_start) & (tvec < min(PlasmaStart - 1e-4, I_start + time_lag) )
#else:
#ind = (tvec > I_start) & (tvec < I_start + time_lag)
#if all(~ind) and ~isnan(CdFieldTrigger): # no intersect
#ind = (tvec > CdFieldTrigger) & (tvec < CdFieldTrigger + time_lag)
#if all(~ind) and ~isnan(BdFieldTrigger): # no intersect
#ind = (tvec > BdFieldTrigger) & (tvec < BdFieldTrigger + time_lag)
#if all(~ind):
#print "!! Something wrong with triggers !!!", CdFieldTrigger, BdFieldTrigger
#ind = tvec < time_lag
#ind = where(ind)[0][::max(1, sum(ind)/100)] # take maximaly 100 points
#print "fitting from t=", tvec[ind[0]], " to t=", tvec[ind[-1]]
##print CdFieldTrigger + 0.1e-3, BdFieldTrigger + 0.1e-3, BdFieldTrigger + time_lag
##plot(tvec, Irog)
##plot(tvec[ind], Irog[ind])
###plot(Uloop)
##show()
##exit()
##all(ind == False)
#def fit(param, result = False):
#Ich = zeros(len(tvec))
#search_range = ( range(ind[0], ind[-1]) if not result else range(1, len(tvec)) )
##search_range = range(1, len(tvec))
#for i in search_range: # implicite differencial equation I_n = ( U_n + C1/dt *I_(n-1) ) / (C2 + C1/dt)
#Ich[i]=(Uloop[i]+param[0]*Ich[i-1])/(param[1]+param[0]);
#if result:
#return Ich
#return (Irog - Ich)[ind]
#x0 = array([ChamberInduktance/dt , ChamberResistance])
#x0 = array([ 1.2 , ChamberResistance])
#try:
##xxx
#x, conv = leastsq(fit, x0 , xtol=1e-3)
#print "fitting done", x, x0
#except:
#print "Fitting failed"
#conv = 0; x = [nan, nan]
#if conv == 0 or (abs(x[0]-x0[0])/x0[0] + abs(x[1]-x0[1])/x0[1] > 5 ): # failed convergence
#print "far from expected value"
#time_lag *= 0.8 # try shorter time
#else:
#time_lag = 0
#if conv == 0 or (abs(x[0]-x0[0])/x0[0] + abs(x[1]-x0[1])/x0[1] > 8 ): # failed convergence
#print "failed convergence in chamber properties solver", x0, x
#x = x0 # if everything fails ...
x = x0
Ich = fit(x, result = True)
#print Ich
saveconst('ChamberInduktance', x[0]*dt)
saveconst('ChamberResistance', x[1])
#print 'Time of fitting %g' % (time.time() - t0)
Ipl = Irog - Ich
#print Irog
#print Ipl
#plot(tvec,Uloop/amax(Uloop))
#plot(tvec,Irog/amax(Irog))
#show()
#plot(Irog)
#plot(Ich)
#plot(Ipl)
#show()
##exit()
# make a deconvolution of the Ip signal, because due to copper shielding and slow response of the rogowski coil, the signal is convoluted with exponential function,
t_exp = 0.7e-4 #roughly estimated influence of the chamber, crosschecked with signal from mirnov coils
win = 60
regularization = 1e-2
Ipl,retrofit = deconvolveExp(Ipl,t_exp,dt,win,regularization)
if PlasmaStart is not None: # second round ..
Ipl1 = median(Ipl[abs(tvec-PlasmaEnd - 1e-3) < 1e-4])
Ipl0 = median(Ipl[abs(tvec-PlasmaStart+5e-4) < 1e-4])
a = (Ipl1 - Ipl0) / (PlasmaEnd - PlasmaStart) # slope of drift
b = Ipl0 - a*PlasmaStart
#Ipl[tvec > PlasmaStart] -= a*tvec[tvec > PlasmaStart] + b
Ipl -= a*tvec + b
print " correction of drift after plasma "
ind = tvec > PlasmaEnd
dIpl = r_[0,diff(Ipl)]
print median(dIpl[ind])
dIpl -= median(dIpl[ind])
Ipl[ind]= cumsum(dIpl[ind])
Ipl[tvec < max(nanmin([TororoidalMagneticFieldTrigger, BdFieldTrigger, CdFieldTrigger, StFieldTrigger])-1e-4 , 2e-3)] = 0
plasma_current_decay = amax(diff(fftconvolve(-Ipl, ones(100)/100, mode='same'))) / dt
save_adv('Iplasma', tvec, Ipl)
save_adv('Ich', tvec, Ich)
saveconst('plasma_current_decay', plasma_current_decay)
return Ipl, Ich
def getPhotod(PlasmaStart,PlasmaEnd, channel, name):
tvec, Photod = S['any', channel]
bias = median(Photod[tvec < nanmin([CdFieldTrigger, BdFieldTrigger, StFieldTrigger]) -1e-4 ]) # remove noise from tyristors
Photod -= bias # constant light in room ...
Photod = medfilt(Photod, 3) # remove fast outliers
## remove light from room (50Hz)
PlasmaEnd += 1e-3
Photod *= sign(mean(Photod[ (tvec>PlasmaStart) & (tvec<PlasmaEnd) ])) # invert signal
try:
ind = (tvec<PlasmaStart)|(tvec>PlasmaEnd)
power_suply = ones((len(tvec), 3))
power_suply[:,1] = cos(2*pi*tvec*100)
power_suply[:,2] = sin(2*pi*tvec*100)
proj = linalg.lstsq(power_suply[ind,:],Photod[ind])[0]
corr = dot(proj.T,power_suply.T ).T
Photod -=corr
except Exception, e :
print "getPhotod: power sockets subtraction failed " , str(e)
if shot > 10527 and shot < 12897 or shot> 13105 : # only fot the new photocells
t = time.time()
if name == "Photod":
Photod,_ = deconvolveExp( Photod,1e-4,dt,300,3)
elif name == "PhotodHalpha":
Photod,_ = deconvolveExp( Photod,0.6e-3,dt,600,3)
Photod *= 1.4
print " dekonvoluce =================", time.time() - t
save_adv(name, tvec, Photod)
getMeanPhotod(Photod, PlasmaStart, PlasmaEnd, name)
return Photod
def PlasmaDetect(Iplasma, dIdt_rogMax, Uloop):
N_steps = len(Iplasma)
Iplasma = Iplasma.copy()
Iplasma[tvec < TororoidalMagneticFieldTrigger + 1e-4] = 0 # no plasma before mag. field
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)
d = convolve( d, ones(20)/20, mode="same" ) # smooth the signal
################
min_dI = 5e4 # sensitivity of plasma detection
####################
d_start = d.copy() # prevent false plasma start detect
d_start[d_start < 0] = 0
try:
# prevent to be plasma false alarm started by triggers !!!
PlasmaStart = tvec[::downsample][where( (d_start > max(0.2*amax(d) , min_dI) ) & ~ ind_wrong[downsample::downsample] )[0][0]]
PlasmaEnd = tvec[::downsample][where( d < min(0.2*amin(d), -min_dI) )[0][-1] ]
PlasmaStart += 0.25e-3 # a fix that is removing efect of detection function "d" smoothing
PlasmaEnd -= 0.25e-3
except:
PlasmaStart = nan
PlasmaEnd = nan
Plasma = int(~ isnan(PlasmaStart) and ~ isnan(PlasmaEnd)) and (PlasmaStart < PlasmaEnd)
if not Plasma:
PlasmaStart = nan
PlasmaEnd = nan
#if Plasma and mean(Uloop[argmin(abs(tvec - PlasmaStart)):argmin(abs(tvec - PlasmaEnd))]) < 0: # negative uloop in plasma
#Plasma = False; PlasmaStart = nan; PlasmaEnd = nan
PlasmaTimeLength = PlasmaEnd - PlasmaStart
saveconst("PlasmaStartAdvanced", PlasmaStart - 0.1e-3 , fmt="%.5f" )
saveconst("PlasmaStart", PlasmaStart, fmt="%.5f" )
saveconst("PlasmaEnd", PlasmaEnd, fmt="%.5f" )
saveconst("PlasmaEndDelayed", PlasmaEnd * 1.1, fmt="%.5f" )
saveconst("Plasma", Plasma )
saveconst("PlasmaTimeLength", PlasmaTimeLength, fmt="%.5f" )
save_adv("PlasmaDetect", tvec[::downsample], d/amax(abs(d))*dIdt_rogMax )
return Plasma, PlasmaStart, PlasmaEnd, PlasmaTimeLength
def getMeanBt(Btor,PlasmaStart, PlasmaEnd ):
MeanBt = mean(Btor[(tvec > PlasmaStart) & (tvec < PlasmaEnd)])
saveconst('BtMean', MeanBt, fmt="%.3f" )
saveconst('ReversedBt', int(MeanBt > 0))
return MeanBt
def getMeanPhotod(Photod, PlasmaStart, PlasmaEnd, name):
MeanPhotod = median(Photod[(tvec > PlasmaStart) & (tvec < PlasmaEnd)])
saveconst(name + "Mean", 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, fmt="%.3f" )
return TotalCharge
def getMeanCurrent(Ipla, PlasmaStart,PlasmaEnd ):
MeanIpla = median(Ipla[(tvec > PlasmaStart) & (tvec < PlasmaEnd)])
saveconst("IplaMean", MeanIpla, fmt="%.1f" )
return MeanIpla
def getMeanUloop(Uloop, PlasmaStart,PlasmaEnd ):
MeanUloop = median(Uloop[(tvec > PlasmaStart) & (tvec < PlasmaEnd)])
saveconst("UloopMean", MeanUloop, fmt="%.2f" )
return MeanUloop
def getOhmicHeatingPower(MeanUloop,MeanPlasmaCurrent):
OhmicHeatingPower = abs(MeanUloop*MeanPlasmaCurrent)
saveconst("OhmicHeatingPowerMean", OhmicHeatingPower, fmt="%.2f" )
return OhmicHeatingPower
def getQedge(MeanBt,MeanPlasmaCurrent):
Qedge = 2*pi*MeanPlasmaRadius**2 / (Mu0 * MajorRadius ) *abs(MeanBt/MeanPlasmaCurrent)
saveconst("QedgeMean", Qedge, fmt="%.1f" )
return Qedge
def getMeanElectronTemperature(MeanUloop, MeanPlasmaCurrent):
#ElectronTemperature = (MajorRadius/MeanPlasmaRadius**2 * 8 * Zeff / 1.54e3)**(2./3) * (abs( MeanPlasmaCurrent / MeanUloop )) ** (2./3) #Ref: PhD Brotankova eq. 3.21 p. 26
ElectronTemperature = (MajorRadius/MeanPlasmaRadius**2 * 8 * Zeff*(1-sqrt(Aspect))**2 / 1.54e3)**(2./3) *abs(MeanPlasmaCurrent/(MeanUloop+1e-3))**(2./3)
saveconst("ElectronTemperatureMean", ElectronTemperature, fmt="%.1f" )
return ElectronTemperature
def getElectronTemperature(Uloop, Iplasma, PlasmaStart, PlasmaEnd):
#corrected for impurities and neoclasical effects.
ind = (tvec > PlasmaStart) & (tvec < PlasmaEnd)
ElectronTempTime = (MajorRadius/MeanPlasmaRadius**2 * 8 * Zeff*(1-sqrt(Aspect))**2 / 1.54e3)**(2./3) *abs(Iplasma/(abs(Uloop)+1e-3))**(2./3)
ElectronTempTime[~ind | (ElectronTempTime > 100 ) ] = nan
ker = sum(ind)/1000
medElectronTemp = medfilt(ElectronTempTime, 2*(ker)+1)
ind = ~isnan(medElectronTemp)
medElectronTemp[ind] = fftconvolve(medElectronTemp[ind], ones(ker)/ker, mode="same")
maxT = mquantiles(medElectronTemp[(tvec > PlasmaStart) & (tvec < PlasmaEnd)],0.99)[0]
save_adv("ElectronTemp", tvec, medElectronTemp)
#save_adv("ElectronMed", tvec, medElectronTemp)
saveconst("ElectronTempMax", maxT, fmt="%.1f" )
return ElectronTempTime, medElectronTemp, maxT
def getBreakDownVoltage(Uloop, Btor, Ipla, PlasmaStart, PlasmaEnd):
dt = (PlasmaStart + PlasmaEnd)/50
ind = (tvec > PlasmaStart-dt) & (tvec < PlasmaStart + dt)
#print PlasmaStart, PlasmaEnd
#plot(Uloop)
#show()
#try:
#t0 = time.time()
#plot(medfilt(Uloop[ind], 21), '.-')
#savefig('breakdown.png')
#clf()
#print "time breakdown " , t0 - time.time()
# kompromis mezi rychlostí a stabilitou
break_ind = where(ind)[0][argmax(medfilt(Uloop[ind], 21))]
#break_ind = argmin( abs(tvec - PlasmaStart ) )
Umax = Uloop[break_ind]
Btime = tvec[break_ind]
Bbreak = Btor[break_ind]
Bipla = Btor[break_ind]
saveconst("BreakDownVoltage", Umax, fmt="%.1f" )
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_old", ChamberResistance)
return ChamberResistance
def Failures(Plasma, UloopMax, dIdt_rogMax, MeanUloop, BtMax, MeanBt, PlasmaStart, PlasmaEnd ):
PlasmaStatus = ""
if abs(UloopMax) < 1:
PlasmaStatus += "Failure (UloopMax < 1V);"
if abs(MeanUloop) < 3 and Plasma:
PlasmaStatus += "Too low Uloop (MeanUloop < 3V);"
if abs(dIdt_rogMax) < 0.01 :
PlasmaStatus += "Failure (dIdt_rogMax < 0.01V);"
if abs(MeanUloop) < 0.3 :
PlasmaStatus += "Failure (MeanUloop) < 0.3V);"
if abs(BtMax) < 0.01 :
PlasmaStatus += "Failure (BtMax) < 0.01T);"
if abs(MeanBt) < 0.01 :
PlasmaStatus += "Failure (MeanBt) < 0.01T);"
if PlasmaStatus == "":
PlasmaStatus = "OK"
print "PlasmaStatus: ", PlasmaStatus
saveconst("PlasmaStatus", PlasmaStatus)
if PlasmaStatus != "OK":
saveconst("PlasmaStartAdvanced", nan , fmt="%.5f" )
saveconst("PlasmaStart", nan, fmt="%.5f" )
saveconst("PlasmaEnd", nan, fmt="%.5f" )
saveconst("PlasmaEndDelayed", nan , fmt="%.5f" )
saveconst("PlasmaTimeLength", nan , fmt="%.5f" )
saveconst("Plasma", 0 )
return PlasmaStatus
def getGreenwaldDensity(Ipla):
N_gw = (1e-6*Ipla)/(pi*MeanPlasmaRadius**2)*1e20 #[m^-3]
save_adv("GreenwaldDensity", tvec, N_gw)
return N_gw
def getOhmicHeatingPowerTime(Ipla, Ich, Uloop):
Power = Ipla * Uloop
PowerCh = Ich * Uloop
save_adv("OhmicHeatingPower", tvec, Power)
save_adv("OhmicHeatingChamber", tvec, PowerCh)
return Power, PowerCh
def getQedgeTime(Btor,Ipla, PlasmaStart,PlasmaEnd ):
QedgeTime = 2*pi*MeanPlasmaRadius**2 / (Mu0 * MajorRadius ) *abs(Btor/(Ipla+1e-3))
QedgeTime[(tvec < PlasmaStart) | (tvec > PlasmaEnd)] = nan
QedgeTime[QedgeTime > 2* nanmedian(QedgeTime)] = nan
save_adv("Qedge",tvec, QedgeTime)
return QedgeTime
def getMagneticFlux(Uloop):
Flux = cumsum(Uloop+abs(Uloop))/2*dt
save_adv("TotalMagneticFlux", tvec, Flux)
return Flux
def getTransformatorSaturation(Uloop, Plasma, PlasmaEnd):
if Plasma == 1:
S = cumsum(Uloop[tvec < PlasmaEnd]) * dt / MaxTransformatorSaturation
else:
S = cumsum(Uloop) * dt / MaxTransformatorSaturation
S = amax(S) # maximal saturation
saveconst('TransformatorSaturation', S)
return S
def EnergyBalance(Ipla, Irog, Uloop, PlasmaStart,PlasmaEnd):
l_i = log(1.65+0.89*nu)
L_i = l_i*Mu0*MajorRadius/(4*pi)
L_e = MajorRadius*Mu0*(log(8*MajorRadius/MeanPlasmaRadius)-2) # wikipedia
L = L_i + L_e
W_e = L_e*Irog**2/2
W_i = L_i*Ipla**2/2
W_mag = W_i+W_e
win = 1000
lam= 1e0
P_mag, retrofit, chi2 = DiffFilter(W_i+W_e, dt,win, lam)
P_mag = squeeze(P_mag)
P_total = Uloop*Irog
P_input = Uloop*Ipla
L_correction = 5 #effective inductance is very different from calculated value, due to coupling with curent drive RLC circuit
sigma_plasma = abs(Ipla/(abs(Uloop)+0.1))
sigma_chamber = 1/9.2e-3 #FIXME conductivity of chamber
P_mag = -P_mag*L_correction
P_plasma = -P_input-P_mag*sigma_plasma/(sigma_plasma+sigma_chamber)
P_chamber = -(P_total-P_input)-P_mag*sigma_chamber/(sigma_plasma+sigma_chamber)
save_adv("PowerTotal", tvec, P_total)
save_adv("PowerMagnetic", tvec, P_mag)
save_adv("PowerPlasma", tvec, P_plasma)
save_adv("PowerChamber", tvec, P_chamber)
saveconst('MeanPowerPlasma', - mean( P_plasma[( tvec > PlasmaStart ) & ( tvec < PlasmaEnd)] ))
return P_total, P_mag, P_plasma, P_chamber
def BreakdownProba():
from pygolem_lite.breakdown import predict
if S['shotno'] > 15170 and S['working_gas'] == "H":
pressure = S['pressure'] / 2.4
else:
pressure = S['pressure']
bd_proba = predict(Ub=S['ub'], Ucd=S['ucd'], pressure=S['pressure'], Tcd=S['tcd'], gas_filling=S['gas_filling'], preionization=S['preionization'] > 0 )
if S['shotno'] > 15170 and S['working_gas'] != "H":
bd_proba = nan
saveconst('breakdown_probability', bd_proba)
def getBreakDownRate(Ipla, Ich, PlasmaStart, PlasmaEnd, trange = 0.4e-3):
ind = (tvec > PlasmaStart - trange/2) & (tvec < PlasmaStart + trange/2) # 200us time window
#ind = where(ind)[::max(1, sum(ind)/100)]
#print ind
down = max(1, sum(ind)/50)
#print "================ down ", down
tvec_tmp = (tvec - PlasmaStart)[ind][::down]
ipla_tmp = smooth(Ipla[ind], down/2.)[::down]
def fit(param, result = False):
ipla_exp = param[0]*exp( 1e3 * tvec_tmp / abs(param[1]) ) + param[2]
if not result:
return (ipla_tmp - ipla_exp)
else:
return ipla_exp
x0 = array([ 50 , 1, 0]) # initial guess
convergence = True
res = leastsq(fit, x0 , xtol=1e-5, full_output=1)
(param, pcov, infodict, errmsg, ier) = res
if ier not in [1,2,3,4]:
print "fitting failed"
convergence = False
print "fitting breakdown done", param, x0
s_sq = sum(fit(param)**2)/(len(ipla_tmp)-len(x0))
pcov = pcov * s_sq
#print pcov
rate = abs(param[1])*1e-3
err = sqrt(pcov[1,1])*1e-3
ind_ext = (tvec > PlasmaStart - 0.4e-3) & (tvec < PlasmaStart + 0.6e-3)
plot(tvec[ind_ext]-PlasmaStart, Ipla[ind_ext])
plot(tvec_tmp, ipla_tmp)
plot(tvec_tmp, fit(param, True))
#plot(tvec_tmp, fit(x0, True), ':')
title('Exp. Rate const= %g+-%g, window=%g' % (rate, err, trange) )
ylabel('I [A]')
xlabel('Time from breakdown [s]')
savefig('rate.png')
savetxt('data', hstack([tvec_tmp[:,None], ipla_tmp[:,None]]), fmt="%g")
if not convergence or rate > 0.01:
print "convergence failed", convergence, rate
rate = nan; err = nan
saveconst('breakdown_rate', rate)
saveconst('breakdown_rate_err', err )
return rate, err
