Discharge/DischargeDatabase/Examples/22813/includes/analysis/Preionization/1113MWPreionization.ON/main.py

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

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

from matplotlib.pyplot import *

from numpy import *
from pygolem_lite.config import *
from pygolem_lite.modules import *
from pygolem_lite import Shot

resonant_Bt = 0.0875
R = MajorRadius
a = MinorRadius

def where_first_or_nan(time, values, harmonic=1):
    possible = time[values >= resonant_Bt*harmonic]
    if len(possible) == 0:
        return nan
    else:
        return possible[0]


def make_plots():
    shot = Shot()
    preionization = shot["preionization"]
    # if preionization not in [4, 5, 6, 7]:
    #     saveconst("MWactive", 0)
    #     return
    # else:
    #     saveconst("MWactive", 1)
    t, mw_power = shot["nistandard", "mw_energy"]
    t, Iplasma = shot["plasma_current"]
    t, Uloop = shot["loop_voltage"]
    t, Bt = shot["toroidal_field"] 
    t, haphoto = shot["haphoto"]

    t *= 1e3

    Bt_LFS = R / (R + a) * Bt
    Bt_HFS = R / (R - a) * Bt

    loc = where(mw_power > mw_power.mean())[0]
    loc = where(mw_power > mw_power.mean())[0]
    mw_start = t[loc[0]]
    os.system(''.join(('echo ',str(mw_start),'>mw_start')))
    mw_end = t[loc[-1]]
    os.system(''.join(('echo ',str(mw_end),'>mw_end')))
                    

    t_resonnant_Bt_HFS = where_first_or_nan(t, Bt_HFS)
    os.system(''.join(('echo ',str(t_resonnant_Bt_HFS),'>t_resonnant_Bt_HFS')))
    t_resonnant_Bt = where_first_or_nan(t, Bt)
    os.system(''.join(('echo ',str(t_resonnant_Bt),'>t_resonnant_Bt')))
    t_resonnant_Bt_LFS = where_first_or_nan(t, Bt_LFS)
    os.system(''.join(('echo ',str(t_resonnant_Bt_LFS),'>t_resonnant_Bt_LFS')))
                        
    
    t_resonnant_2_Bt_HFS = where_first_or_nan(t, Bt_HFS, 2)
    t_resonnant_2_Bt = where_first_or_nan(t, Bt, 2)
    t_resonnant_2_Bt_LFS = where_first_or_nan(t, Bt_LFS, 2)

    def mark_regions():
        axvline(mw_start, c="g", ls="--")
        axvline(mw_end, c="g", ls="--")    
        axvline(t_resonnant_Bt, c="r", ls="--")
        axvline(t_resonnant_Bt_LFS, c="r", ls="--")
        axvline(t_resonnant_Bt_HFS, c="r", ls="--")

        axvline(t_resonnant_2_Bt, c="r", ls=":")
        axvline(t_resonnant_2_Bt_LFS, c="r", ls=":")
        axvline(t_resonnant_2_Bt_HFS, c="r", ls=":")

        
    figure(0)
    axes1 = subplot(511)
    axes1.plot(t, Bt_HFS, "m-", label=r"$B_t$ HFS")    
    axes1.plot(t, Bt, "b-", label=r"$B_t$ center")
    axes1.plot(t, Bt_LFS, "c-", label=r"$B_t$ LFS")
    ylabel(r"$B_t$ [T]")
    ylim(0, Bt_HFS.max()*1.1)
    axes1.axhline(resonant_Bt, c="k", label=("%.4f T" % resonant_Bt))
    axes1.axhline(resonant_Bt * 2, linestyle='--', c="k", label=("%.4f T" % (resonant_Bt * 2)))
    mark_regions()
    l = legend(loc=1, fancybox=True, prop={"size": 10})
    #l.get_frame().set_alpha(0.75)
    axes2 = subplot(512, sharex=axes1)
    axes2.plot(t, Uloop)
    ylabel(r"$U_{loop}$ [V]")
    u_max = shot["loop_voltage_max"]
    ylim(0, u_max+ 0.5)
    mark_regions()
    axes3 = subplot(513, sharex=axes1)
    axes3.plot(t, Iplasma)
    ylabel(r"$I_{plasma}$ [A]")
    ylim(ymin=0)
    mark_regions()
    axes4 = subplot(514, sharex=axes1)
    axes4.plot(t, mw_power)
    ylabel(r"$P_{MW}$ [a.u.]")
    mark_regions()
    axes5 = subplot(515, sharex=axes1)
    axes5.plot(t, haphoto)
    ylabel(r"$H_{\alpha}$ intensity [a.u.]")
    xlabel(r"$t$ [ms]")
    mark_regions()
    savefig("okenko.png")


if __name__ == "__main__":
    make_plots()
    saveconst('status', 0)