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 import *

from matplotlib.pyplot import *
from scipy.optimize import leastsq
import time
from scipy.sparse import spdiags, eye
from scipy.sparse.linalg import spsolve, splu
from scipy.signal import medfilt, fftconvolve
from scipy.stats.mstats import mquantiles

print "------------basicdiagn --------"

# get time vector 
tvec, tmp = Shot()['NIstandard']

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


def save_config():
    names = [
    'ToroidalMagneticFieldCoilInductance',
    'TororoidalMagneticFieldCapacitor',
    'BreakDownElectricFieldCapacitor',
    'CurrentDriveElectricFieldCapacitor',
    'StabilizationCapacitor',
    'Zeff'
    ]
    for i in names:
	exec('saveconst("'+i+'", '+i+')')





#def getHXR():
    ##tvec, data = Shot()['any', 'hxr']
    #tvec, data = Shot()['hxr']

    #data[data < 0] = 0

    #sigma = len(tvec)/800   # !!! smoothing parameter !!!!
    #ker = arange(2*sigma, dtype=double)
    #ker = hstack([ker[::-1], ker[1:]])/sigma
    #ker = exp( - ker **2)
    #ker /= sum(ker)
    #data = fftconvolve(data, ker, mode="same")
    #save_adv('HXR_smooth', tvec[::4], data[::4])


    #return data




def getBtoroidal():
    tvec, dBt = Shot()['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 = Shot()['NIstandard', 'uloop']

    #m_Uloop = mean(Uloop)
    Uloop *= UloopCalibration

    ind_zero =  (tvec < TororoidalMagneticFieldTrigger  + 1e-4)
    Uloop[ind_zero] = 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
    #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 ) ]

    #print Uloop[ind_bias]

    #print TororoidalMagneticFieldTrigger, nanmin([CdFieldTrigger, StFieldTrigger])

    Uloop_bias = median(Uloop[ind_bias])

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

def getIrogowski():

    tvec, dIrog = Shot()['NIstandard', 'irog']

    #dIrog[abs(tvec - TororoidalMagneticFieldTrigger) < 1e-4 ] = 0    # cut off triggers

    Ibias = mean(dIrog[tvec <   max(nanmin([TororoidalMagneticFieldTrigger, BdFieldTrigger, CdFieldTrigger, StFieldTrigger])-1e-4 , 2e-3)] )
    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(diff(Irog),0.99)[0]  # maximum
    IrogMax = mean(Irog)

    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()

    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 < PlasmaStart - 1e-4)
	else:
	    ind = (tvec >  I_start) & (tvec < I_start + time_lag)
	if all(~ind):  # no intersect
	    ind = (tvec >  CdFieldTrigger) & (tvec < CdFieldTrigger + time_lag)
	if all(~ind):  # no intersect
	    ind = (tvec >  BdFieldTrigger) & (tvec < BdFieldTrigger + time_lag)

	if  all(~ind):
	    print "!! Something wrong with triggers !!!"

	ind = where(ind)[0][::max(1, sum(ind)/100)]  # take maximaly 100 points


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

	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] > 0.5 ): # failed convergence
	    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] > 0.8 ): # failed convergence
	print "failed convergence in chamber properties solver", x0, x
	x = x0  # if everything fails ...

    Ich = fit(x, result = True)

    saveconst('ChamberInduktance', x[0]*dt)
    saveconst('ChamberResistance', x[1])

    #print 'Time of fitting %g' % (time.time() - t0)



    Ipl = Irog - Ich

    #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

    #plot(Ipl)
    #show()

    save_adv_par('Iplasma', tvec, Ipl)
    save_adv_par('Ich', tvec, Ich)
    return Ipl, Ich

def getPhotod(PlasmaStart, channel, name):

    tvec, Photod = Shot()['NIstandard', channel]
    bias = median(Photod[tvec < nanmin([CdFieldTrigger, BdFieldTrigger, StFieldTrigger]) -1e-4 ])   # remove noise from tyristors
    Photod -= bias
    #Photod[ind_wrong] = medfilt(Photod, 51)[ind_wrong]
    Photod *= sign(mean(Photod))  # invert  signal
    save_adv_par(name, tvec, Photod)
    return Photod


#def getPhotodHalpha(PlasmaStart):

    #tvec, PhotodHalpha = Shot()['NIstandard', 'haphoto']

    #bias = median(PhotodHalpha[tvec < nanmin([CdFieldTrigger, BdFieldTrigger, StFieldTrigger]) -1e-4 ])   # remove bias from tyristors
    #PhotodHalpha -= bias

    ##PhotodHalpha[ind_wrong] = medfilt(PhotodHalpha, 51)[ind_wrong]

    #PhotodHalpha *= sign(mean(PhotodHalpha))  # invert  signal

    #save_adv_par('PhotodHalpha', tvec, PhotodHalpha)

    #return PhotodHalpha

def PlasmaDetect(Iplasma):
    polar = sign(mean(Iplasma))
    ReversedCD = (polar == -1)
    saveconst("ReversedCD", ReversedCD)

    N_steps = len(Iplasma)

    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
    ####################
    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]]
    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" )

    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):

    MeanPhotod = median(Photod[(tvec > PlasmaStart) & (tvec < PlasmaEnd)])

    saveconst("PhotodMean", 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

    medElectronTemp = medfilt(ElectronTempTime, 2*(sum(ind)/1000)+1)

    maxT = mquantiles(medElectronTemp[(tvec > PlasmaStart) & (tvec < PlasmaEnd)],0.99)[0]

    save_adv_par("ElectronTemp", tvec, ElectronTempTime)
    save_adv_par("ElectronMed", tvec, medElectronTemp)
    saveconst("ElectronTempMax", maxT, fmt="%.1f" )



    return ElectronTempTime, medElectronTemp, maxT


def getBreakDownVoltage(Uloop, Btor, Ipla, PlasmaStart, PlasmaEnd):

    ind = (tvec > PlasmaStart) & (tvec < PlasmaStart + (PlasmaStart + PlasmaEnd)/10)

    #print PlasmaStart, PlasmaEnd
    #plot(Uloop)
    #show()


    break_ind  = where(ind)[0][argmax(Uloop[ind])]
    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.1 :
	PlasmaStatus += "Failure (dIdt_rogMax < 0.1V);"
    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)

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