External stabilization

Current shot

Scan

import numpy as np
import matplotlib.pyplot as plt
import requests
from scipy import integrate, signal, interpolate
import pandas as pd
import holoviews as hv
hv.extension('bokeh')
import hvplot.pandas
import os, glob
ds = np.DataSource()

def remove_old(name):
    if os.path.exists(name):
        os.remove(name)
        
names = ['analysis.html', 'scan.html', 'icon.html', 'icon-fig.png','icon-fig_Radial.png','icon-fig_Vertical.png','Camera_vs_Mirnov_Radial.png','Camera_vs_Mirnov__Vertical.png','InnerQuadr.csv']
for name in names:
    remove_old(name)
    
if os.path.exists('Results/'):
    file=glob.glob('Results/*')
    for f in file:
        os.remove(f)

shot_no = 39501
def open_remote(shot_no, identifier, url_template):
    return ds.open(url_template.format(shot_no=shot_no, identifier=identifier))

def read_signal(shot_no, identifier, url, data_type='csv'): 
    file = open_remote(shot_no, identifier, url)
    if data_type == 'lvm':
        channels = channels = ['Time', 'mc1','mc5','mc9','mc13','Saddle','InnerQuadr', 'IexRad','IexVert','9', '10', '11', '12']
        return pd.read_table(file, sep = '\t', names=channels, index_col = 'Time')
    else:
        return pd.read_csv(file, names=['Time', identifier],
                     index_col = 'Time', squeeze=True)
def read_tab(shot_no,url):
    tab=ds.open(url)
    df=pd.read_csv(tab)
    end=df['Time'].iat[-1]
    start=df['Time'].iat[0]
    df=df.set_index('Time')
    return df, start, end
    
def remove_offset(data, window):
    data-=data.loc[:window].mean()
    return data

def smooth(data,win=11): #41
    smooth_data = signal.savgol_filter(data, win, 3)
    return smooth_data

url = f'http://golem.fjfi.cvut.cz/shots/{shot_no}/Diagnostics/LimiterMirnovCoils/plasma_position.csv'
kI=1/0.05 #Rogowski coil (Flux loop) constant


try:
    df,start,end=read_tab(shot_no,url)
    mirnov=True
except OSError:
    mirnov=False
    
try:
    
    #!!! Prohozene I s U
    url_I_vert  = f'http://golem.fjfi.cvut.cz/shots/{shot_no}/Operation/Discharge/Position_Stabilization/DAS_raw_data_dir_vertical/U%5evertical_fg.csv'
    url_U_vert  = f'http://golem.fjfi.cvut.cz/shots/{shot_no}/Operation/Discharge/Position_Stabilization/DAS_raw_data_dir_vertical/U%5evertical_currclamp.csv'
    
    url_I_rad  = f'http://golem.fjfi.cvut.cz/shots/{shot_no}/Operation/Discharge/Position_Stabilization/DAS_raw_data_dir_radial/U%5eradial_fg.csv'
    url_U_rad  = f'http://golem.fjfi.cvut.cz/shots/{shot_no}/Operation/Discharge/Position_Stabilization/DAS_raw_data_dir_radial/U%5eradial_currclamp.csv'
    
    U_exStabVert = read_signal(shot_no, 'U_fg_Vert', url_U_vert)
    
    url_dt=requests.get(f'http://golem.fjfi.cvut.cz/shotdir/{shot_no}/Devices/Oscilloscopes/RigolMSO5204-b/ScopeSetup/XINC')

    dt=float(url_dt.text)*1e5
    print('dt:',dt)
    t_shift=-1
    t_osc = pd.Series(np.linspace(0,dt*10, len(U_exStabVert))).rename('Time')+t_shift
    
    U_exStabVert = pd.Series(U_exStabVert.index[:], index = t_osc)
    RogCoilVert = read_signal(shot_no, 'I_fg_Vert', url_I_vert)
    I_exStabVert = pd.Series(smooth(RogCoilVert.index[:]*kI), index = t_osc) #the data was multiplied by the consntant [V->A]
#     I_exStabVert = pd.Series(smooth(RogCoilVert.index[:]), index = t_osc) #the data is already recalculated
         
    U_exStabRad = read_signal(shot_no, 'U_fg_Rad', url_U_rad)
    
    U_exStabRad = pd.Series(U_exStabRad.index[:], index = t_osc)
    RogCoilRad = read_signal(shot_no, 'I_fg_Rad', url_I_rad)
    I_exStabRad = pd.Series(smooth(RogCoilRad.index[:]*kI), index = t_osc) #the data was multiplied by the consntant [V->A]
#     I_exStabRad = pd.Series(smooth(RogCoilRad.index[:]), index = t_osc) #the data is already recalculated
        
    exStab=True
    
except OSError:
    
    exStab=False

dt: 2.0

Stabilization

if exStab:
    names=['$I_{vert}$','$I_{rad}$','$U_{vert}$','$U_{rad}$','Vertical', 'Radial']
    Stab=pd.DataFrame({'$I_{vert}$': I_exStabVert, '$I_{rad}$': I_exStabRad, 
                       '$U_{vert}$': U_exStabVert, '$U_{rad}$': U_exStabRad})
    maxI=0
    for stab in range(2):
        I=names[stab]
        U=names[stab+2]
        
        fig, ax = plt.subplots(dpi=100)
        ax = Stab[I].plot(grid=True,label = I, c = 'r') #ocilloscope data
        ax2 = ax.twinx()
        ax2 = Stab[U].plot(label = '$U_{fg}$',c='tab:blue')
        ax.set(ylim=(-max(abs(Stab[I]))-0.5, max(abs(Stab[I]))+0.5), ylabel='I [A]',
               title=f'Shot No.{shot_no}, {names[stab+4]} Stabilization' )
        ax2.set(ylim=(-max(abs(Stab[U]))-0.5, max(abs(Stab[U]))+0.5), ylabel='U [V]')
        ax.legend(loc='upper left');ax2.legend(loc='upper right')
        
        if mirnov==True:
            plt.savefig('icon-fig_'+names[stab+4]+'.png')
        else:
            plt.savefig('icon-fig.png')
            
        maxI_stab=max(abs(Stab[I]))
        if maxI<maxI_stab:
            maxI=maxI_stab
        print('Maximum current:', round(maxI_stab,2), 'A')

        with open('Results/MaxI_'+names[stab+4], 'w') as f:
            f.write(str(round(maxI_stab,2)))

Maximum current: 97.55 A
Maximum current: 99.88 A

if exStab:
    df_processed = pd.concat([Stab], axis='columns')
    df_processed

if exStab:
    data_URL = f"http://golem.fjfi.cvut.cz/shots/{shot_no}/Diagnostics/LimiterMirnovCoils/DAS_raw_data_dir/NIdata_6133.lvm"
    dataNI = read_signal(shot_no, ' ', data_URL, 'lvm').replace([np.inf, -np.inf, np.nan], value = 0)
    dataNI.index*=1e3
    Ivert=dataNI['IexVert']*kI
    Irad=dataNI['IexRad']*kI
    InnerQuadr=dataNI['InnerQuadr']
    
    InnerQuadr.to_csv('Results/InnerQuadr.csv')
    df_processed = pd.concat([Ivert, Irad],axis='columns')
    if mirnov:
        df_processed = pd.concat([df, Ivert, Irad],axis='columns')
    names=['IexVert','IexRad']
df_processed.to_csv('Results/ExStab.csv')    

Plasma position - Mirnov coils

print('Position from Mirnov coils:', mirnov)
if mirnov:
    hline = hv.HLine(0)
    hline.opts(color='k', line_dash='dashed', alpha = 0.4, line_width=1.0)
    I_v=names[0]; I_r=names[1]
    if exStab: 
        ax_param=dict(ylim=(-85, 85),xlim=(start, end),height=250,width=600, grid = True)
        axIvert_param=dict(title = '',ylabel='I_exStab [A]', xlabel='Time [ms]', yaxis='left', height=250, width=600, 
             color='b', ylim=(-maxI-1.5, maxI+1.5),grid=True)

        axIrad_param=dict(title = '',ylabel='I_exStab [A]', xlabel='Time [ms]', yaxis='left', height=250, width=600, 
             color='r', ylim=(-maxI-1.5,maxI+1.5),grid=True)

        plot = df_processed['r'].hvplot(title = '', ylabel='r [mm]/ I[A]', xlabel='', **ax_param) *\
                    df_processed[I_v].hvplot(**axIvert_param) *\
                        df_processed[I_r].hvplot(**axIrad_param)+\
        df_processed['z'].hvplot(title = '',ylabel='z [mm] / I[A]',  xlabel='Time [ms]', **ax_param) *\
            df_processed[I_v].hvplot(**axIvert_param) *\
                df_processed[I_r].hvplot(**axIrad_param)+\
        df_processed['a'].hvplot(title = '',ylabel='a[mm]/ I[A]', xlabel='Time [ms]', **ax_param) *\
            df_processed[I_v].hvplot(**axIvert_param) *\
                    df_processed[I_r].hvplot(**axIrad_param)
        plot=plot*hline
        plot.cols(2)
        
    else:
        layout = hv.Layout([df_processed[v].hvplot.line(
            xlabel='', ylabel=l,ylim=(-85,85), xlim=(start,end),legend=False, title='', grid=True, group_label=v)
                            for (v, l) in [('r', ' r [mm]'), ('z', 'z [mm]'), ('a', 'a [mm]')] ])*hline
        plot=layout.cols(1).opts(hv.opts.Curve(width=600, height=200),  hv.opts.Curve('a', xlabel='time [ms]'))

    hvplot.save(plot, 'icon.html')

Position from Mirnov coils: True

New iconfig

if mirnov:
    fig, axs = plt.subplots(3, 1, sharex=True, dpi=200)
    df_processed['r'].plot(grid=True, ax=axs[0])
    df_processed['z'].plot(grid=True, ax=axs[1])
    df_processed['a'].plot(grid=True, ax=axs[2])
    for i in range(3):
        ax4=axs[i].twinx()
        ax4=df_processed[names[0]].plot(grid=True, c='b',alpha=0.3,label=names[0])
        ax4=df_processed[names[1]].plot(grid=True, c='r',alpha=0.3,label=names[1])
        plt.legend()
        ax4.set(xlim=(start,end), ylim=(-2*maxI,maxI*2),xlabel= 'Time [ms]', ylabel = 'I$_{exStab}$ [A]')
        
    axs[2].set(ylim=(-85,85), xlim=(start,end), xlabel= 'Time [ms]', ylabel = '$a$ [mm]')
    axs[1].set(ylim=(-85,85), xlim=(start,end), xlabel= 'Time [ms]', ylabel = '$\Delta$z [mm]')
    axs[0].set(ylim=(-85,85), xlim=(start,end), xlabel= 'Time [ms]', ylabel = '$\Delta$r [mm]', 
               title = 'Horizontal, vertical plasma position and radius #{}'.format(shot_no))
    plt.savefig('icon-fig')

Scan over the session

try:
    ScanDef_url = requests.get("http://golem.fjfi.cvut.cz/shots/%i/Production/Parameters/ScanDefinition" % shot_no)
    ScanDef = ScanDef_url.text
    ScanDef=ScanDef.split(" ")
    for i in range(len(ScanDef)):
        ScanDef[i]=int(ScanDef[i])
    Scan=True
except ValueError:
    Scan=False
    ScanDef=[]
    print('Scan not defined')

#TODO - osa x, zredukovat pocet smycek
if Scan and mirnov:
    for shot in ScanDef:
        url=f'http://golem.fjfi.cvut.cz/shots/{shot}/Diagnostics/LimiterMirnovCoils/plasma_position.csv'
        df2,start,end=read_tab(shot_no, url)
        param = ['r','z', 'a']
        for j in param:
            df_processed = df_processed.join(df2[j],rsuffix=f'_{shot}')
            
        url=f'http://golem.fjfi.cvut.cz/shots/{shot}/Operation/Discharge/Position_Stabilization/Results/ExStab.csv'
        tab=ds.open(url)
        Stab2=pd.read_csv(tab).set_index('Time')  
        
        for name in range(2):
            df_processed = df_processed.join(Stab2[names[name]], rsuffix=f'_{shot}')
     
    ax_param=dict(xlim=(start,end),ylim=(-85, 85),height=250,width=600, grid = True)
    axStab=dict(xlim=(start,end),ylim=(-110,110),height=250,width=600, grid = True)
    hline = hv.HLine(0)
    hline.opts(color='k', line_dash='dashed', alpha = 0.4, line_width=1.0)
    
    display(df_processed)
    plotScan = df_processed['r'].hvplot(title = '', ylabel='r [mm]', xlabel='', **ax_param) +\
                df_processed['z'].hvplot(title = '',ylabel='z [mm]',  xlabel='', **ax_param) +\
                df_processed['a'].hvplot(title = '',ylabel='a[mm]', xlabel='Time [ms]', **ax_param) +\
                df_processed[names[1]].hvplot(title = '',ylabel='I_rad[A]', xlabel='Time [ms]',**axStab)+\
                df_processed[names[0]].hvplot(title = '',ylabel='I_vert[A]', xlabel='Time [ms]',**axStab) 
                
    for shot in ScanDef:
        plotScan=plotScan[0]*df_processed['r_%i'%shot].hvplot(title = '', ylabel='r [mm]', xlabel='', **ax_param)+\
            plotScan[1]*df_processed['z_%i'%shot].hvplot(title = '', ylabel='z [mm]', xlabel='', **ax_param)+\
            plotScan[2]*df_processed['a_%i'%shot].hvplot(title = '', ylabel='a [mm]', xlabel='Time', **ax_param)+\
            plotScan[3]*df_processed[names[1]+'_%i'%shot].hvplot(title = '',ylabel='I_rad[A]', xlabel='Time [ms]',**axStab) +\
            plotScan[4]*df_processed[names[0]+'_%i'%shot].hvplot(title = '',ylabel='I_vert[A]', xlabel='Time [ms]',**axStab)
                    
        plotScan=plotScan*hline
        plotScan.cols(3)
        hvplot.save(plotScan, 'scan.html')

	r	z	a	IexVert	IexRad	r_39492	z_39492	a_39492	IexVert_39492	IexRad_39492	r_39497	z_39497	a_39497	IexVert_39497	IexRad_39497
Time
0.000	NaN	NaN	NaN	-0.04882	-0.04882	NaN	NaN	NaN	0.02442	0.04882	NaN	NaN	NaN	0.14648	0.21972
0.001	NaN	NaN	NaN	0.04882	0.04882	NaN	NaN	NaN	0.04882	0.12208	NaN	NaN	NaN	0.04882	0.14648
0.002	NaN	NaN	NaN	0.00000	-0.04882	NaN	NaN	NaN	-0.02442	0.04882	NaN	NaN	NaN	0.04882	0.19532
0.003	NaN	NaN	NaN	-0.04882	-0.04882	NaN	NaN	NaN	0.09766	0.12208	NaN	NaN	NaN	-0.02442	0.12208
0.004	NaN	NaN	NaN	-0.07324	-0.07324	NaN	NaN	NaN	0.02442	0.04882	NaN	NaN	NaN	0.00000	0.12208
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
39.995	NaN	NaN	NaN	-10.15624	-8.61816	NaN	NaN	NaN	0.09766	0.09766	NaN	NaN	NaN	-10.18066	-3.24708
39.996	NaN	NaN	NaN	-10.20508	-8.56934	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN
39.997	NaN	NaN	NaN	-10.35156	-8.74024	NaN	NaN	NaN	0.02442	0.04882	NaN	NaN	NaN	-10.20508	-3.17382
39.998	NaN	NaN	NaN	-10.35156	-8.74024	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN
39.999	NaN	NaN	NaN	-10.30274	-8.74024	NaN	NaN	NaN	0.07324	0.09766	NaN	NaN	NaN	-10.13184	-3.14942

76566 rows × 15 columns

def scan():
    tdur=[]; StabVertParam=[]; StabRadParam=[]; MaxIvertStab=[]; MaxIradStab=[]
    ScanDef.append(shot_no)
    for shot in ScanDef:
        url_tdur=requests.get('http://golem.fjfi.cvut.cz/shots/%i/Diagnostics/BasicDiagnostics/Results/t_plasma_duration'%shot)
        tdur.append(float(url_tdur.text))

        url_StabVert=requests.get('http://golem.fjfi.cvut.cz/shots/%i/Operation/Discharge/Position_Stabilization/Parameters/vertical_waveform'%shot)
        StabVertParam.append(url_StabVert.text) 
        url_StabRad=requests.get('http://golem.fjfi.cvut.cz/shots/%i/Operation/Discharge/Position_Stabilization/Parameters/radial_waveform'%shot)
        StabRadParam.append(url_StabRad.text) 
        
        url_maxIvert=requests.get('http://golem.fjfi.cvut.cz/shots/%i/Operation/Discharge/Position_Stabilization/Results/MaxI_Vertical'%shot)
        MaxIvertStab.append(url_maxIvert.text)
        url_maxIrad=requests.get('http://golem.fjfi.cvut.cz/shots/%i/Operation/Discharge/Position_Stabilization/Results/MaxI_Radial'%shot)
        MaxIradStab.append(url_maxIrad.text)
    data={'Plasma duration': tdur, 'Vertical Stabilization parameters': StabVertParam, 
            'Radial Stabilization parameters': StabRadParam, 'Max current in VertStab': MaxIvertStab,
              'Max current in RadStab': MaxIradStab, 'Shot_no': ScanDef}
    df_results=pd.DataFrame(data)
    df_results=df_results.set_index('Shot_no')
    df_results.to_csv('Results/scanResults.csv')
    return df_results

if mirnov and Scan:
    display(scan())

	Plasma duration	Vertical Stabilization parameters	Radial Stabilization parameters	Max current in VertStab	Max current in RadStab
Shot_no
39492	9.332	3000,0;11000,-5;22000,0;30000,0\n	2000,0;3000,0;9000,0;12000,-5;24000,0;25000,0\n	6.63	6.98
39497	16.935	3000,0;6000,-5;8000,-20;30000,0\n	2000,0;3000,0;8000,-5;12000,-10;24000,-2;25000...	97.78	52.66
39501	17.056	3000,0;6000,-5;8000,-20;30000,0\n	2000,0;8000,-5;10000,-10;12000,-15;14000,-20;3...	97.55	99.88

Cameras vs Mirnov Coils

def comparison_image(shot_no, Position):
    
    url=f'http://golem.fjfi.cvut.cz/shots/{shot_no}/Diagnostics/LimiterMirnovCoils/plasma_position.csv'
    tab=ds.open(url)
    df=pd.read_csv(tab)
    end=df['Time'].iat[-1]
    start=df['Time'].iat[0]
    df=df.set_index('Time')
    
    url_VertCam = f'http://golem.fjfi.cvut.cz/shots/{shot_no}/Diagnostics/FastCameras/Camera_Vertical/CameraVerticalPosition'
    url_RadCam = f'http://golem.fjfi.cvut.cz/shots/{shot_no}/Diagnostics/FastCameras/Camera_Radial/CameraRadialPosition'
    
    if Position=='Radial':
        symb='r'
        camera_position=read_signal(shot_no, symb, url_RadCam, data_type='csv')
    if Position=='Vertical': 
        symb='z'
        camera_position=read_signal(shot_no, symb, url_VertCam, data_type='csv')

    fig = plt.figure(figsize=(12,4),dpi=70)
    
    FONT = 'DejaVu Sans'

    ax = df[symb].plot(label='Mirnov coils')
    ax = camera_position.plot(label = 'Fast Camera')

    ax.set_ylabel('$\Delta$'+symb+' [mm]',fontname = FONT, fontsize = 12)
    ax.set_xlabel('Time [ms]',fontname = FONT, fontsize = 12)
    ax.set_ylim(-85,85) 

    loclegend='best'
    leg = plt.legend(loc = loclegend, shadow = True, fancybox=False) #with marker

    leg.get_frame().set_linewidth(1)
    leg.get_frame().set_edgecolor('k')
    for text in leg.get_texts():
          plt.setp(text, fontname=FONT, fontsize = 12)
    for line, text in zip(leg.get_lines(), leg.get_texts()):
          text.set_color(line.get_color())

    ax.grid(which = 'major', c = 'gray', linewidth = 0.5, linestyle = 'solid') 
    ax.grid(which = 'minor', c = 'gray', linewidth = 0.3, linestyle = 'dashed') 

    ax.axhline(y=0, color='k', ls='--', lw=1, alpha=0.4)

    fig.savefig(f'Camera_vs_Mirnov_{Position}')
    

for position in ['Radial','Vertical']:
    try:
        comparison_image(shot_no, position)
    except:
        print(f'{position} camera data is missing')