Source code :: main

#!/usr/bin/python2
# -*- coding: utf-8 -*-


""" CREATED: 7/2012
    AUTHOR: MICHAL ODSTRČIL
"""

print "importing modules "


import matplotlib 
matplotlib.rcParams['backend'] = 'Agg'
matplotlib.rc('font',  size='10')
matplotlib.rc('text', usetex=True)  # FIXME !! nicer but slower !!!

from numpy import *
from pygolem_lite.config import *
from pygolem_lite.modules import *
from pygolem_lite import Shot
import sys
from scipy.stats.mstats import mquantiles
import time

print "importing modules done "


def plot_data(file_type):

    #saveconst('status', 1)

    S = Shot()
    shot = S.shot_num
    
    plot_params = dict( figsize = (9,7), file_type = file_type)
    plasma = S['plasma']    

    if plasma:
      start = S['plasma_start']    
      end = S['plasma_end']    
    else:
      start = 0
      end = 40e-3
      
      
      
    # !!!! try to load data during the second basicdiagn reload after all diagnostics !!!!!
    t0 = time.time()
    max_HXR = 0
    HXR_const = 0.005
    if S.exist('hxr_smooth'):
	[tvec_HXR, HXR] =  S['hxr_smooth'] 
	HXR *= HXR_const# 25  #  max(mquantiles(HXR,0.995)*1.2,1)   ##  !! do not renormalize only noise (data < 0.2)
	max_HXR = mquantiles(HXR[(tvec_HXR > start) & (tvec_HXR < end)] , 0.99)

    max_photo = 0.3
    for sig in ['photodiode', 'photodiode_alpha', 'photodiode_other']:
	try:
	    tvec, photodiode = S[sig]
	    max_photo = max(max_photo, mquantiles(photodiode[(tvec > start) & (tvec < end)] , 0.99) )
	except:
	    pass
    ph_range = [0, max( max_photo*1.2 ,max_HXR*1.2 )] # do not allow smaller Yrange than 0-1

    shot_title = 'Golem shot No:' + str(shot)
    
    data = [
	get_data('loop_voltage', 'Loop voltage', 'U [V]', ylim = [0,None]),
	get_data('toroidal_field', 'Toroidal mag. field', 'B$_t$ [T]' , ylim=[0,None],  reduction = True),
	get_data('plasma_current', 'Plasma current', 'I$_{p}$ [kA]' ,  data_rescale = 1e-3, reduction = True) if plasma else \
	get_data('rogowski_current', 'Chamber current', 'I$_{ch}$ [kA]', data_rescale = 1e-3,  reduction = True) ,
	[get_data('photodiode', 'Visible', 'Intensity [a.u.]' , ylim = ph_range),
	get_data('photodiode_alpha', 'H$_\\alpha$', 'Intensity [a.u.]' , ylim = ph_range)]+
	[get_data('photodiode_other', 'Other', 'Intensity [a.u.]' , ylim = ph_range) ] +
	[get_data('hxr_smooth', 'HXR', 'Intensity [a.u.]' , ylim = ph_range, data_rescale = HXR_const,  reduction = True) if S.exist('hxr_smooth') and S['hxr_mean'] > 0.01 else None ], 
	[get_data('electron_density', 'Electron density ($n_e$)', 'n$_{e}$ [10$^{19}\cdot$m$^{-3}$]' ,   data_rescale= 1e-19 ) \
	if (plasma and S.exist('electron_density') and S['electron_density_mean'] > 1e16 and S['electron_density:reliability'] < 0.25) else None ]
	]
	
	
    print "start plotting"
    multiplot(data, shot_title, 'graphpres',   **plot_params)    
      
    print "end plotting"
    from time import gmtime, strftime
    print strftime("%Y-%m-%d %H:%M:%S", gmtime())


	
    # GPf
    data = [
	get_data('loop_voltage', 'Loop voltage', 'U [V]'),
	get_data('toroidal_field', 'Toroidal mag. field', 'B$_t$ [T]', reduction = True),
	get_data('rogowski_current', 'Total current', 'I$_{ch}$ [kA]' , data_rescale = 1e-3, reduction = True),
	[get_data('photodiode', 'Visible', 'Intensity [a.u.]' ),
	get_data('photodiode_alpha', 'H$_\\alpha', 'Intensity [a.u.]' )] ,
	]
    
    paralel_multiplot(data, shot_title, 'graphprint',  **plot_params)



    data = [
	get_data('loop_voltage', 'Loop voltage', 'U [V]',  xlim = [0,None]),
	get_data('toroidal_field', 'Toroidal mag. field', 'B$_t$ [T]' , xlim = [0,None], reduction = True),
	[get_data('plasma_current', 'Plasma current', 'I [kA]' , xlim = [0,None], data_rescale = 1e-3, reduction = True), 
	get_data('chamber_current', 'Chamber current', 'I [kA]',xlim = [0,None], data_rescale = 1e-3, reduction = True), 
	get_data('rogowski_current', 'Total current', 'I [kA]',xlim = [0,None], data_rescale = 1e-3, reduction = True)],  
	[get_data('photodiode', 'Visible', 'Intensity [a.u.]', xlim = [0,None]),
	get_data('photodiode_alpha', 'H$_\\alpha', 'Intensity [a.u.]' , xlim = [0,None])] ,
	]

    paralel_multiplot(data, shot_title + ' - final data ', 'graphpresfull',   **plot_params)

   
    #GPitegrated
    data = [
	get_data('loop_voltage', 'Loop voltage', 'U [V]',  xlim = [0,None]),
	get_data('toroidal_field', 'Toroidal mag. field', 'B$_t$ [T]',  xlim = [0,None], reduction = True),
	[get_data('rogowski_current', 'Total current', 'I$_p$+I$_{ch}$ [kA]' ,  xlim = [0,None], data_rescale = 1e-3, reduction = True),  \
	get_data('plasma_current', 'Plasma current', 'I$_{p}$ [kA]' , xlim = [0,None], data_rescale = 1e-3, reduction = True)] if plasma else \
	get_data('rogowski_current', 'Total current', 'I$_{ch}$ [kA]' ,  xlim = [0,None], data_rescale = 1e-3, reduction = True) ,
	[get_data('photo', 'Visible', '[DAS V]' , xlim = [0,None] ),
	get_data('haphoto', 'H$_\\alpha', '[DAS V]',  xlim = [0,None] )] 
	]

    paralel_multiplot(data, shot_title + " - integrated data", 'graphpresi',  **plot_params)

    ####raw

    data = [
	get_data('uloop', 'Loop voltage', '[DAS V]', xlim = [0,None] ),
	get_data('btor', 'Derivative  of mag. field', 'dB$_t$/dt [DAS V]', xlim = [0,None] ),
	[get_data('irog', 'Raw Rogowski signal', 'dI$_{p+ch}$/dt [DAS V]',  xlim = [0,None] ),
	get_data('PlasmaDetect', 'Derivative of I$_{p}$', 'dI$_{p}$/dt [a.u.]',  xlim = [0,None] )], 
	[get_data('photo', 'Visible', '[DAS V]' , xlim = [0,None] ),
	get_data('haphoto', 'H$_\\alpha', '[DAS V]',  xlim = [0,None] )]
	]

    paralel_multiplot(data, shot_title + " - raw data", 'graphpresb',  **plot_params)


    # icon
    data = [
	get_data('loop_voltage', '', '', xlabel = ""),
	[get_data('rogowski_current', '', '', xlabel = "" )] +\
	[get_data('plasma_current','', '', xlabel = "" )] if plasma else \
	[get_data('rogowski_current', '', '', xlabel = "" )]
	]
	
    paralel_multiplot(data, '', 'icon', (4,3), 40)


    # GPic
    data = [
	get_data('loop_voltage', '', 'U [V]',   xlim=[0,None], ylim = [0,None], xlabel=""),
	[get_data('rogowski_current', '', 'I$_p$+I$_{ch}$ [kA]' ,  xlim=[0,None],  ylim = [0,None], xlabel = "", data_rescale = 1e-3, reduction = True), 
	get_data('plasma_current', '', 'I$_{p}$ [kA]' , xlim=[0,None],  ylim = [0,None], xlabel = "",  data_rescale = 1e-3, reduction = True)] \
	if plasma else get_data('rogowski_current', '', 'I$_{ch}$ [kA]' , xlim=[0,None],   ylim = [0,None], xlabel = "", data_rescale = 1e-3, reduction = True)
	]
    
    paralel_multiplot(data, "", 'graphic', (9,1.5), 100,  'horizontal' )

    print ' Time: %g' %  ( time.time() - t0)

    saveconst('status', 0)

def prepare_data():
    
    from basic_diagn import *
    
    t = time.time()

    Aktual_PfeifferMerkaVakua = loadconst("Aktual_PfeifferMerkaVakua")
    
    getDate()

    save_config()
    
    [Btor, dBtor, BtMax, BtMean] = getBtoroidal()
    
    [Uloop, UloopMax, UloopMean, ReversedCD] = getUloop()
   
    [Irog, dIdt_rogMax, IrogMax, I_start] = getIrogowski(ReversedCD)

    [Ipla, Ich] = getIplasma(Uloop,Irog, I_start)
    
    Plasma, PlasmaStart, PlasmaEnd, PlasmaTimeLength = PlasmaDetect(Ipla, dIdt_rogMax)
    
    if Plasma:
	[Ipla, Ich] = getIplasma(Uloop,Irog,I_start, PlasmaStart, PlasmaEnd )  # second iteration (hopefully better)
    
    for name, fname in zip(['photovrb', 'haphoto', 'photo'], ['PhotodVrba',  'PhotodHalpha','Photod']):
	try:
	    photodiode = getPhotod(PlasmaStart, PlasmaEnd, name, fname)
	    getMeanPhotod(photodiode, PlasmaStart, PlasmaEnd, fname)
	except Exception, e:
	    print "Photodiode " + name + " failed " , str(e)

    MeanBt = getMeanBt(Btor,PlasmaStart, PlasmaEnd )
    TotalCharge  = getTotalCharge(Ipla, PlasmaStart,PlasmaEnd )

    MeanUloop = getMeanUloop(Uloop, PlasmaStart,PlasmaEnd )
    MeanIpla = getMeanCurrent(Ipla, PlasmaStart,PlasmaEnd )
    OhmicHeatingPower = getOhmicHeatingPower(MeanUloop,MeanIpla)
    
    PlasmaStatus  = Failures(Plasma, UloopMax, dIdt_rogMax, MeanUloop, BtMax, MeanBt, PlasmaStart, PlasmaEnd  )
        
    if Plasma:
	Qedge = getQedge(MeanBt,MeanIpla)
	ElectronTemperature = getMeanElectronTemperature(MeanUloop, MeanIpla)	
	Umax, Btime, Bbreak, Bipla = getBreakDownVoltage(Uloop, Btor, Ipla, PlasmaStart, PlasmaEnd)
	StateEqElectronDensity = getStateEqElectronDensity(Aktual_PfeifferMerkaVakua)
	ElectronConfinementTimeFirstApprox = getElectronConfinementTimeFirstApprox(MeanUloop,MeanIpla, StateEqElectronDensity, ElectronTemperature )
	ElectronTempTime, medElectronTemp, maxT = getElectronTemperature(Uloop, Ipla, PlasmaStart, PlasmaEnd)
	GreenwaldDensity = getGreenwaldDensity(Ipla)
	Qedge = getQedgeTime(Btor,Ipla,   PlasmaStart,PlasmaEnd)
	tranges = array([0.05, 0.1,0.15, 0.2, 0.3, 0.4])*1e-3
	nt = len(tranges)
	err = zeros(nt)
	for i in range(nt):
	   rate, err[i] =  getBreakDownRate(Ipla, Ich, PlasmaStart, PlasmaEnd, tranges[i])
	   print "win %g rate %g err %g" % (tranges[i], rate, err[i])
	err[isnan(err)] = inf
	getBreakDownRate(Ipla, Ich, PlasmaStart, PlasmaEnd, tranges[argmin(err)])

    ChamberResistance = getChamberResistance(Plasma)

    OhmicHeatingPower, OhmicHeatingChamber = getOhmicHeatingPowerTime(Ipla, Ich, Uloop)

    MagneticFlux = getMagneticFlux(Uloop)
    Saturation = getTransformatorSaturation(Uloop, Plasma, PlasmaEnd)
    EnergyBalance(Ipla, Irog, Uloop, PlasmaStart, PlasmaEnd)

    try:
	BreakdownProba()
    except Exception, e:
	print "BreakdownProba failed",  str(e)
    
    print "time of basic diagn generation" , time.time() - t

    
def main():
    if sys.argv[1] ==  "acquisition":
	prepare_data()
    elif sys.argv[1] ==  "plots":  
	plot_data('png')
	#plot_data('svgz')
	saveconst('status', 0)
    elif  sys.argv[1] ==  "postanalysis":  
	plot_data('png')  # replot and include diagnostics 
	plot_data('svgz')
	saveconst('status', 0)
  
  


if __name__ == "__main__":
    main()