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 = 37452
shot_no = 37542
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
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: 94.67 A Maximum current: 51.88 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')
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')
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),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')
True
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')
Scan=False
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:
df2_rename=df2[j].rename(j+'_%i'%shot)
df_processed = pd.concat([df_processed, df2_rename], axis='columns')
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)
Stab2=Stab2.set_index('Time')
for name in range(2):
Stab2_rename=Stab2[names[name]].rename(names[name]+'_%i'%shot)
df_processed = pd.concat([df_processed, Stab2_rename], axis='columns')
ax_param=dict(ylim=(-85, 85),height=250,width=600, grid = True)
axStab=dict(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[namees[0]+'_%i'%shot].hvplot(title = '',ylabel='I_vert[A]', xlabel='Time [ms]',**axStab)
plotScan=plotScan*hline
plotScan.cols(3)
hvplot.save(plotScan, 'scan.html')
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())
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')
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=100)
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')
Vertical camera data is missing