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
from IPython import embed as keyboard
from  multiprocessing import Process


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.025 # 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 - filetype: %s" % file_type
    multiplot(data, shot_title, 'graphpres',   **plot_params)    
      
    print "end basic graphs 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 plotting all : %.2gs' %  ( time.time() - t0)

    saveconst('status', 0)

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

    Aktual_PfeifferMerkaVakua = loadconst("Aktual_PfeifferMerkaVakua")

    getDate()  #  421 us per loop

    save_config()  # 967 us per loop

    [Btor, dBtor, BtMax, BtMean] = getBtoroidal() # 17.2 ms per loop
    [Uloop, UloopMax, UloopMean, ReversedCD] = getUloop() #  17.2 ms per loop  
    [Irog, dIdt_rogMax, IrogMax, I_start] = getIrogowski(ReversedCD) # 27.1 ms per loop
    [Ipla, Ich] = getIplasma(Uloop,Irog, I_start)  # 111 ms per loop
    
    Plasma, PlasmaStart, PlasmaEnd, PlasmaTimeLength = PlasmaDetect(Ipla, dIdt_rogMax, Uloop)  # 3.41 ms per loop
    
    if Plasma:
	[Ipla, Ich] = getIplasma(Uloop,Irog,I_start, PlasmaStart, PlasmaEnd )  # second iteration (hopefully better) # 112 ms per loop
    
    for name, fname in zip(['haphoto', 'photo'], ['PhotodHalpha','Photod']):
	try:
	    Process(target=getPhotod, args = (PlasmaStart, PlasmaEnd, name, fname) ).start() # 38.6 ms per loop
	except Exception, e:
	    print "Photodiode " + name + " failed " , str(e)


    MeanBt = getMeanBt(Btor,PlasmaStart, PlasmaEnd )  #  447 us per loop
    MeanUloop = getMeanUloop(Uloop, PlasmaStart,PlasmaEnd ) # 376 us per loop
    MeanIpla = getMeanCurrent(Ipla, PlasmaStart,PlasmaEnd ) # 277 us per loop
    
    getOhmicHeatingPower(MeanUloop,MeanIpla) #  146 us per loop
    getTotalCharge(Ipla, PlasmaStart,PlasmaEnd ) # 284 us per loop

    
    Failures(Plasma, UloopMax, dIdt_rogMax, MeanUloop, BtMax, MeanBt, PlasmaStart, PlasmaEnd  )  # 201 us per loop

        
    if Plasma:
	ElectronTemperature = getMeanElectronTemperature(MeanUloop, MeanIpla)	 # 166 us per loop
	StateEqElectronDensity = getStateEqElectronDensity(Aktual_PfeifferMerkaVakua) # 314 us per loop
	getElectronConfinementTimeFirstApprox(MeanUloop,MeanIpla, StateEqElectronDensity, ElectronTemperature ) # 319 us per loop
	getElectronTemperature(Uloop, Ipla, PlasmaStart, PlasmaEnd) # 23.4 ms per loop
	getQedge(MeanBt,MeanIpla)  # 157 us per loop
	getBreakDownVoltage(Uloop, Btor, Ipla, PlasmaStart, PlasmaEnd) # 1.01 ms per loop
	getGreenwaldDensity(Ipla) #  19.8 ms per loop
	getQedgeTime(Btor,Ipla, PlasmaStart,PlasmaEnd) # 21.2 ms per loop
	#tranges = array([0.05, 0.1,0.15, 0.2, 0.3, 0.4])*1e-3
	#Process(target=find_breakdown_rate, args = (Ipla, Ich, PlasmaStart, PlasmaEnd, tranges) ).start()

    getChamberResistance(Plasma)  # 165 us per loop
    getOhmicHeatingPowerTime(Ipla, Ich, Uloop)  # 388 us per loop
    getMagneticFlux(Uloop) # 22.4 ms per loop
    getTransformatorSaturation(Uloop, Plasma, PlasmaEnd)  #  374 us per loop
    EnergyBalance(Ipla, Irog, Uloop, PlasmaStart, PlasmaEnd) # 189 ms per loop

    #keyboard()


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

def find_breakdown_rate(Ipla, Ich, PlasmaStart, PlasmaEnd, tranges):
    from basic_diagn import getBreakDownRate
    nt = len(tranges)
    err = zeros(nt)
    for i in range(nt):
	rate, err[i] =  getBreakDownRate(Ipla, Ich, PlasmaStart, PlasmaEnd, tranges[i])  # 42.5 ms per loop  * 6 !! 
	print "win %g rate %g err %g" % (tranges[i], rate, err[i])
    err[isnan(err)] = inf
    getBreakDownRate(Ipla, Ich, PlasmaStart, PlasmaEnd, tranges[argmin(err)])  # 44 ms per loop
    
def main():
    if sys.argv[1] ==  "acquisition":
	t = time.time()
	prepare_data()
	print "Time prepare_data %g s" % (time.time() - t)
	
    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()