#!/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()
