#!/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) ) 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 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 = mean(Bt) < 0 saveconst("ReversedB", ReversedB) Bt = -Bt if ReversedB else Bt 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 = mean(Uloop) < 0 saveconst("ReversedCD", ReversedCD) Uloop = -Uloop if ReversedCD else Uloop #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() 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 ... 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 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 Photod *= sign(mean(Photod)) # invert signal ## remove light from room (50Hz) PlasmaEnd += 1e-3 ind = (tvecPlasmaEnd) try: 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: pass #if shot > 10527: # 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 bd_proba = predict(Ub=S['ub'], Ucd=S['ucd'], pressure=S['pressure'], Tcd=S['tcd'], gas_filling=S['gas_filling'], preionization=S['preionization'] ) 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