Source code :: basic_diagn

#!/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 --------"

# get time vector 
S = Shot()
tvec, tmp = S['NIstandard']
shot = S.shot_num
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


### 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

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['NIstandard', '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

    BtMax = amax(Bt)      # max
    BtMean = mean(Bt)     # mean

    saveconst('BtMax', BtMax, fmt="%.3f" )
    #saveconst('BtMean', BtMean)
    save_adv_par('Btoroidal', tvec, Bt)
    save_adv_par('dBdt_toroidal', tvec, dBt)
    saveconst('noise_level', sBt)


    return Bt, dBt, BtMax, BtMean

def getUloop():
    tvec, Uloop = S['NIstandard', '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)
    if ReversedCD:
	Uloop *= -1
    #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[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_par('Uloop', tvec, Uloop)
    return Uloop, UloopMax, UloopMean, ReversedCD

def getIrogowski(ReversedCD):

    tvec, dIrog = S['NIstandard', '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_par('Irogowski', tvec, Irog)
    save_adv_par('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_par('Iplasma', tvec, Ipl)
    save_adv_par('Ich', tvec, Ich)
    saveconst('plasma_current_decay', plasma_current_decay)

    return Ipl, Ich

def getPhotod(PlasmaStart,PlasmaEnd,  channel, name):

    tvec, Photod = S['NIstandard', 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 = (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

    
    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_par(name, tvec, Photod)
    return Photod


def PlasmaDetect(Iplasma, dIdt_rogMax):


    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

    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_par("ElectronTemp", tvec, medElectronTemp)
    #save_adv_par("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(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_par("GreenwaldDensity", tvec, N_gw)
    return N_gw

def getOhmicHeatingPowerTime(Ipla, Ich,  Uloop):
    Power =   Ipla * Uloop
    PowerCh = Ich * Uloop
    save_adv_par("OhmicHeatingPower", tvec, Power)
    save_adv_par("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_par("Qedge",tvec,  QedgeTime)

    return QedgeTime

def getMagneticFlux(Uloop):
    Flux = cumsum(Uloop+abs(Uloop))/2*dt
    save_adv_par("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_par("PowerTotal", tvec,  P_total)
    save_adv_par("PowerMagnetic", tvec,  P_mag)
    save_adv_par("PowerPlasma", tvec,  P_plasma)
    save_adv_par("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