
import numpy as np
import matplotlib.pyplot as plt
import requests
from scipy import integrate, signal, interpolate
import pandas as pd
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','IexStab.png']
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 = 37153
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', 'mc_out','IexStab','9', '10', '11', '12']
        return pd.read_table(file, sep = '\t', names=channels)
    else:
        return pd.read_csv(file, names=['Time', identifier],
                     index_col = 'Time', squeeze=True)
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:
    tab=ds.open(url)
    df=pd.read_csv(tab)
    end=df['Time'].iat[-1]
    start=df['Time'].iat[0]
    df=df.set_index('Time')
    mirnov=True
except OSError:
    mirnov=False
    
try:
    
    url_U_vert  = f'http://golem.fjfi.cvut.cz/shots/{shot_no}/Operation/Discharge/Position_Stabilization/DAS_raw_data_dir_vertical/U^vertical_fg.csv'
    url_I_vert  = f'http://golem.fjfi.cvut.cz/shots/{shot_no}/Operation/Discharge/Position_Stabilization/DAS_raw_data_dir_vertical/U^vertical_currclamp.csv'
    
    url_U_rad  = f'http://golem.fjfi.cvut.cz/shots/{shot_no}/Operation/Discharge/Position_Stabilization/DAS_raw_data_dir_radial/U^radial_fg.csv'
    url_I_rad  = f'http://golem.fjfi.cvut.cz/shots/{shot_no}/Operation/Discharge/Position_Stabilization/DAS_raw_data_dir_radial/U^radial_currclamp.csv'
    
    U_exStabVert = read_signal(shot_no, 'U_fg_Vert', url_U_vert)
    
    dt = 2 #time scale setting on the oscilloscope (...ms per window)
    t_osc = pd.Series(np.linspace(0,dt*10, len(U_exStabVert))).rename('Time')
    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)
    
    dt = 2 #time scale setting on the oscilloscope (...ms per window)
    t_osc = pd.Series(np.linspace(0,dt*10, len(U_exStabRad))).rename('Time')
    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

print(f'External Stabilization:{exStab}')

External Stabilization:True


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.legend(loc='upper left')
        ax.set_ylabel('I [A]')
        ax.set_ylim= (-max(abs(Stab[I]))-0.5, max(abs(Stab[I]))+0.5)
        ax2.legend(loc='upper right')
        ax2.set_ylabel('U [V]')
        ax2.set_ylim(-max(abs(Stab[U]))-0.5, max(abs(Stab[U]))+0.5)
        ax.set_title(f'Shot No.{shot_no}, {names[stab+4]} Stabilization')
        if mirnov==True:
            plt.savefig('icon-fig_'+names[stab+4]+'.png')
        else:
            plt.savefig('icon-fig.png')
        #     maxI=max(abs(IexStab_cut)) #use with NI
        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: 6.11 A
Maximum current: 3.28 A


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


if exStab and mirnov:
    Stab_cut=Stab.loc[start:end]
    df_processed = pd.concat([df, Stab_cut],axis='columns')

    #replace NaN value
    k_I = 0
    k_U = 0

    for stab in range(2):
        I=names[stab]
        U=names[stab+2]
        
        for i in range(len(df_processed[I])):
            
            if not np.isnan(df_processed[I].iat[i]):
                
                k_I = df_processed[I].iat[i]
                k_U = df_processed[U].iat[i]
             
            else:
                
                df_processed[I].iat[i]=k_I
                df_processed[U].iat[i]=k_U
    
elif mirnov:
    df_processed = df    
df_processed

df_processed.to_csv('Results/ExStab.csv')


import holoviews as hv
hv.extension('bokeh')
import hvplot.pandas
print(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),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')
# plot

External stabilization¶

Current shot

Scan

Stabilization¶