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]

Navigation