#!/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) if S.exist('electron_density'): tvec_electron_density, electron_density= S['electron_density'] electron_density = interp(tvec, tvec_electron_density, electron_density, left=0, right=0) 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([tvec, 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()