Source code :: main
[Return]
[Download]#!/usr/bin/python2
# -*- coding: utf-8 -*-
""" CREATED: 7/2012
AUTHOR: MICHAL ODSTRČIL
"""
import matplotlib
matplotlib.rcParams['backend'] = 'Agg'
matplotlib.rc('font', size='10')
matplotlib.rc('text', usetex=True) # FIXME !! nicer but slower !!!
import pygolem_lite
from pygolem_lite import Shot
from numpy import *
from pygolem_lite.config import *
from pygolem_lite.modules import *
from matplotlib.pyplot import *
import time
from shutil import copy, move
import os
from scipy.interpolate import interp1d
def analysis():
try:
Data = Shot()
I = Data['plasma_current_mean'] * 1e-6
Bt = Data['toroidal_field_mean']
P = Data['input_power_mean']*1e-3
n = Data['electron_density_mean'] * 1e-19
M = 1 # hydrogen !!
R = MajorRadius
eps = Aspect
kappa = 1
const = 5.62e-2
print "I",I
print "B", Bt
print "P", P
print "n", n
print "R", R
print "Aspect", Aspect
t_98 = const*I**0.93*Bt**0.15*P**-0.69*n**0.41*M**0.18*R**1.97*eps**0.58*kappa**0.78
print "confinent time = ", t_98 * 1e6
saveconst('t_98', t_98)
except Exception, e:
print " confinent time failed ", e.__str__()
def plot_data(file_type):
data = Shot()
plasma = data['plasma']
start = data['plasma_start']*1e3
end = data['plasma_end']*1e3
nGW_limit = 3 # n/nGW < 3
calb_photo = 10 # guess of absolute calbration of photodiode !!
if not plasma:
t_cd, I_cd = Shot()['current_cd_coils'] # current drive coils current
#plot(t_cd, I_cd)
#savefig('I_cd.png')
#clf()
_, I_cd, _ = DiffFilter(I_cd, mean(diff(t_cd)), 1000, 1e2) # smooth the data
t_PHI, PHI = Shot()['total_magnetic_flux']
I_cd = interp(t_PHI, t_cd, squeeze(I_cd), left=0, right=0)
data = [[get_data('total_magnetic_flux', 'Total magnetic flux $\Phi$', '$\Psi$ [Vs]', xlim = [0, 40], reduction = True),
get_data([array([0, 40])*1e-3, ones(2)*MaxTransformatorSaturation], 'Max saturation', '$\Psi$ [Vs]' , ylim = [0,0.18], xlim = [0, 40], reduction = True)],
get_data( [t_PHI, I_cd] , 'Current drive I', 'I [A]', xlabel = 'Time [ms]', xlim = [0, 40], tvec_rescale=1e3),
]
paralel_multiplot(data, 'Total magnetic flux', 'magnetic_flux', (9,6), 100, 'vertical', file_type)
plot( I_cd,PHI)
#plot( PHI)
savefig('I.png')
clf()
data = get_data( [I_cd, PHI ] , 'Hysteresis curve', '$\Phi$ [Vs]', xlabel = 'I [A]',xlim=[amin(I_cd), amax(I_cd)], plot_limits = False, tvec_rescale=1)
paralel_multiplot(data, 'Hysteresis curve', 'hysteresis', (9,3), 100, 'vertical', file_type)
paralel_multiplot(data, "" , 'icon', (4,3), 40)
return
t0 = time.time()
tvec, Q95 = Shot()['Q95']
try:
ne_corr = r0/a*(pos_part(1-R/r0))**1.5
except:
tvec_n_e = tvec
try:
tvec_r, plasma_radius = Shot()['plasma_radius']
plasma_radius_f = interp1d(tvec_r, plasma_radius, fill_value = 0, bounds_error = False)
corr = (plasma_radius_f(tvec) / MeanPlasmaRadius ) ** 2 # corrected for plasma readius
corr[ corr < 0.3 ] = nan
if all(isnan(corr)):
del corr
except:
print "correction on plasma position failed"
tvec_gw, nGW = Shot()['greenwald_density']
nGW_med = max(nGW[(tvec_gw > start*1e-3) & (tvec_gw < end*1e-3)])
try:
tvec_n_e, n_e = Shot()['electron_density']
except:
print "electron_density failed"
try:
tvec_r, R_position = Shot()['plasma_position_r']
plasma_R_pos_f = interp1d(tvec_r, R_position, fill_value = 0, bounds_error = False)
a = 0.1 #[m]
r0 = plasma_radius_f(tvec_n_e)
R = plasma_R_pos_f(tvec_n_e)-0.4
pos_part = lambda x: (x+abs(x))/2
ne_corr_nu1 = r0/a*(pos_part(1-R/r0))**1.5 #peaking factor 1; profile = (1-(r/a)^2)
ne_corr_nu0_5 = r0/a*(pos_part(1-R/r0)) #peaking factor 0.5; profile = sqrt(1-(r/a)^2)
ne_corr_nu1[ne_corr_nu1<0.2] = nan
ne_corr_nu0_5[ne_corr_nu0_5<0.2] = nan
if all(isnan(ne_corr_nu0_5)):
del ne_corr_nu0_5
del ne_corr_nu1
except:
print "electron_density correction failed"
Q95[Q95 > 2* median(Q95[(tvec > start*1e-3) & (tvec < end*1e-3)])] = nan
try:
temp_spec = Shot()['spectrometr\temperature']
power_spec = Shot()['spectrometr\total_radiated_power']
except:
print "missing data spectrometer"
data_1 = [
[
get_data([tvec,Q95], 'Safety factor edge (Q$_{edge}$)', 'Q [-]' , ylim = [0,None] , xlim = [start, end], reduction = True),
get_data([tvec,Q95/3], '(Q$_{center}$) (peaking factor = 2)', 'Q [-]' , ylim = [0,None], xlim = [start, end], reduction = True),
get_data([tvec,corr*Q95], 'Q edge with radius', 'Q [-]' , ylim = [0,None] , xlim = [start, end], reduction = True, line_format='b--') if 'corr' in locals() else [],
get_data([tvec,corr*Q95/3], 'Q center with radius', 'Q [-]' , ylim = [0,None] , xlim = [start, end], reduction = True, line_format='r--') if 'corr' in locals() else []
],
[
get_data([tvec,1/Q95], 'Inverse Safety factor $\iota_{edge}$', '$\iota$ [-]' , ylim = [0,None], xlim = [start, end]),
get_data([tvec,3/Q95], ' $\iota_{center}$ (peaking factor = 2)', '$\iota$ [-]' , ylim = [0,None], xlim = [start, end], reduction = True),
get_data([tvec,1/(Q95*corr)], ' $\iota_{edge}$ with radius', '$\iota$ [-]' , ylim = [0,None], xlim = [start, end], line_format='b--', reduction = True) if 'corr' in locals() else [],
get_data([tvec,3/(Q95*corr)], ' $\iota_{center}$ with radius', '$\iota$ [-]' , ylim = [0,None], xlim = [start, end], line_format='r--', reduction = True) if 'corr' in locals() else [],
get_data([array([start,end])*1e-3, array([0.5, 0.5])], 'Q=2', '$\iota$ [-]' , ylim = [0,0.55], xlim = [start, end], reduction = True, line_format = "k:")
]
]
multiplot(data_1, 'Safety factor', 'plot_1', (9,5), 100, 'vertical', file_type)
data_2 = [[
get_data('electron_temperature_median', 'Electron temperature', 'T [eV]' , ylim = [0,None] , xlim = [start, end], reduction = True),
get_data('electron_temperature_median', '$T_e$ corrected by plasma radius', 'T [eV]' , data_rescale = 1/corr, ylim = [0,None] , xlim = [start, end], line_format='r--', reduction = True) if 'corr' in locals() else [],
get_data(temp_spec, '$T_e$ spectrometer', 'T [eV]' , ylim = [0,None] , xlim = [start, end], line_format='r--') if 'temp_spec' in locals() else [],
#get_data('electron_temperature', 'Electron temperature', 'T [eV]' , ylim = [0,None], xlim = [start, end])
]]
try:
multiplot(data_2, 'Electron temperature', 'plot_2', (9,4), 100, 'vertical', file_type)
except Exception, e:
print 'Electron temperature failed', e.__str__()
data_3 = [[
get_data('input_power_total', 'Ohmic heating - total', 'P [kW]' , xlim = [start, end], data_rescale=1e-3, reduction = True),
get_data('input_power_magnetic', 'Ohmic heating - magnetic', 'P [kW]' , xlim = [start, end], data_rescale=-1e-3, reduction = True),
get_data('input_power_plasma', 'Ohmic heating - plasma', 'P [kW]' , xlim = [start, end], data_rescale=-1e-3, reduction = True),
get_data('input_power_chamber', 'Ohmic heating - chamber', 'P [kW]' , xlim = [start, end], data_rescale=-1e-3, reduction = True),
get_data('photodiode', 'Output power (diod)', 'P [kW]' , xlim = [start, end], data_rescale= calb_photo , reduction = True, line_format = "k--"),
get_data(power_spec, 'Output power (spectrometer)', 'P [kW]' , xlim = [start, end], data_rescale=1e-3 , line_format = "r--") if 'power_spec' in locals() else [],
get_data([array([start,end])*1e-3, array([0, 0])], '', 'P [kW]' , xlim = [start, end], reduction = True, line_format = "k:")
]]
try:
multiplot(data_3, 'Input / Output power', 'plot_3', (9,4), 100, 'vertical', file_type)
except Exception, e:
print 'Input / Output power', e.__str__()
if 'n_e' in locals():
data_4 = [[
get_data('electron_density', 'Electron density ($n_e$)', 'n$_{e}$ [10$^{19}\cdot$m$^{-3}$]' , ylim = [0,None] , xlim = [start, end], data_rescale= 1e-19) if 'n_e' in locals() else [],
get_data('electron_density', '$n_e$ - position corrected, $\\nu = 1$', 'n$_{e}$ [10$^{19}\cdot$m$^{-3}$]' , ylim = [0,None] , xlim = [start, end], data_rescale= 1e-19/ne_corr_nu1, line_format='k--') if ('ne_corr_nu1' in locals()) and ('n_e' in locals() ) else [],
get_data('electron_density', '$n_e$ - position corrected, $\\nu = 0.5$', 'n$_{e}$ [10$^{19}\cdot$m$^{-3}$]' , ylim = [0,None] , xlim = [start, end], data_rescale= 1e-19/ne_corr_nu0_5, line_format='k:') if ('ne_corr_nu0_5' in locals()) and ('n_e' in locals()) else [],
get_data('greenwald_density', 'Greenwald density ($n_{GW}$)', 'n$_{GW}$ [10$^{19}\cdot$m$^{-3}$]' , ylim = [0,None] , xlim = [start, end], data_rescale= 1e-19, reduction = True),
get_data('greenwald_density', 'Greenwald density with radius ($n_{GW}$)', 'n$_{GW}$ [10$^{19}\cdot$m$^{-3}$]' , ylim = [0,None] , xlim = [start, end], data_rescale=1e-19/corr, reduction = True, line_format='r--') if 'corr' in locals() else [],
get_data([array([start,end])*1e-3, array([nGW_med, nGW_med])*nGW_limit], '~GW limit', 'n$_{e}$ [10$^{19}\cdot$m$^{-3}$]' , ylim = [0,None] , xlim = [start, end],data_rescale= 1e-19 , reduction = True, line_format = "k:") if 'n_e' in locals() else [],
]]
try:
multiplot(data_4, 'Density', 'plot_4', (9,4), 100, 'vertical', file_type)
except Exception, e:
print 'Density', e.__str__()
#os.system('convert -resize 150x120\! total_plot.png icon.png')
print "===== plotting done ===== "
data = data_1 + data_2 +data_3 + (data_4 if 'n_e' in locals() else [])
paralel_multiplot(data, "" , 'icon', (4,3), 40)
def main():
if sys.argv[1] == "analysis":
analysis()
if sys.argv[1] == "plots":
plot_data('png')
saveconst('status', 0)
if sys.argv[1] == "postanalysis":
plot_data('png')
saveconst('status', 0)
if __name__ == "__main__":
main()[Return]
