import numpy as np
import matplotlib.pyplot as plt
import requests
from scipy import integrate, signal, interpolate
import pandas as pd
import os
def remove_old(name):
if os.path.exists(name):
os.remove(name)
names = ['analysis.html', 'scan.html', 'icon.html', 'icon-fig.png','IexStab.png']
for name in names:
remove_old(name)
ds = np.DataSource()
# shot_no=35991
shot_no=0
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 = ['Time', 'mc1', 'mc5', 'mc9', 'mc13', '5', 'RogQuadr', '7', '8', '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)
url='http://golem.fjfi.cvut.cz/shots/%i/Diagnostics/LimiterMirnovCoils/plasma_position.csv'%shot_no
tab=ds.open(url)
df=pd.read_csv(tab)
end=df['Time'].iat[-1]
start=df['Time'].iat[0]
df=df.set_index('Time')
try:
url_U = "http://golem.fjfi.cvut.cz/shots/%i/Devices/Oscilloscopes/RigolMSO5104-a/ch2" %shot_no
url_I = "http://golem.fjfi.cvut.cz/shots/%i/Devices/Oscilloscopes/RigolMSO5104-a/ch4" %shot_no
U_exStab = read_signal(shot_no, 'U_fg', url_U)
dt = 2 #time scale setting on the oscilloscope (...ms per window)
t_osc = pd.Series(np.linspace(0,dt*10, len(U_exStab))).rename('Time')
U_exStab = pd.Series(U_exStab.index[:], index = t_osc)
kI=1/0.05
RogCoil = read_signal(shot_no, 'I_fg', url_I)
I_exStab = pd.Series(RogCoil.index[:]*kI, index = t_osc)
exStab=True
except OSError:
exStab=False
print('External Stabilization:', exStab)
if exStab==True:
IexStab_cut = I_exStab.loc[start:end]
UexStab_cut = U_exStab.loc[start:end]
df_processed = pd.concat([df, IexStab_cut.rename('I_exStab'), UexStab_cut.rename('U_exStab')], axis='columns')
#replace NaN value
k_I = 0
k_U = 0
for i in range(len(df_processed['I_exStab'])):
if not np.isnan(df_processed['I_exStab'].iat[i]):
k_I = df_processed['I_exStab'].iat[i]
k_U = df_processed['U_exStab'].iat[i]
else:
df_processed['I_exStab'].iat[i]=k_I
df_processed['U_exStab'].iat[i]=k_U
else:
df_processed = df
# df_processed
External Stabilization: True
import holoviews as hv
hv.extension('bokeh')
import hvplot.pandas
hline = hv.HLine(0)
hline.opts(color='k', line_dash='dashed', alpha = 0.4, line_width=1.0)
if exStab==True:
ax_param=dict(ylim=(-85, 85),height=250,width=600, grid = True)
axI_param=dict(title = '',ylabel='I_exStab [A]', xlabel='Time [ms]', yaxis='left', height=250, width=600,
color='r', ylim=(-max(abs(df_processed['I_exStab']))-1.5,max(abs(df_processed['I_exStab']))+1.5),grid=True)
plot = df_processed['r'].hvplot(title = '', ylabel='r [mm]/ I[A]', xlabel='', **ax_param) *\
df_processed['I_exStab'].hvplot(**axI_param)+\
df_processed['z'].hvplot(title = '',ylabel='z [mm] / I[A]', xlabel='Time [ms]', **ax_param) *\
df_processed['I_exStab'].hvplot(**axI_param)+\
df_processed['a'].hvplot(title = '',ylabel='a[mm]/ I[A]', xlabel='Time [ms]', **ax_param) *\
df_processed['I_exStab'].hvplot(**axI_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
if exStab==True:
fig, ax = plt.subplots(dpi=100)
ax = I_exStab.plot(grid=True,label = '$I_{exStab}$', c = 'r')
ax2 = ax.twinx()
ax2 = U_exStab.plot(label = '$U_{fg}$')
ax.legend(loc='upper left')
ax.set_ylabel('I [A]')
ax.set_ylim= (-max(abs(I_exStab))-0.5, max(abs(I_exStab))+0.5)
ax2.legend(loc='upper right')
ax2.set_ylabel('U [V]')
ax2.set_ylim(-max(abs(U_exStab))-0.5, max(abs(U_exStab))+0.5)
plt.savefig('IexStab')
maxI=max(abs(IexStab_cut))
print('Maximum current:', round(maxI,2), 'A')
with open('Results/Max_IexStab', 'w') as f:
f.write(str(round(maxI,2)))
Maximum current: 9.27 A
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')
if Scan==True:
for shot in ScanDef:
url='http://golem.fjfi.cvut.cz/shots/%i/Diagnostics/LimiterMirnovCoils/plasma_position.csv'%shot
tab=ds.open(url)
df2=pd.read_csv(tab)
end=df2['Time'].iat[-1]
start=df2['Time'].iat[0]
df2=df2.set_index('Time')
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')
ax_param=dict(ylim=(-85, 85),height=250,width=600, grid = True)
hline = hv.HLine(0)
hline.opts(color='k', line_dash='dashed', alpha = 0.4, line_width=1.0)
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)
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=plotScan*hline
plotScan.cols(1)
hvplot.save(plotScan, 'scan.html')
tdur=[]
StabParam=[]
MaxIstab=[]
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_Stab=requests.get('http://golem.fjfi.cvut.cz/shots/%i/Production/Parameters/diagnostics_stabilization'%shot)
StabParam.append(url_Stab.text)
url_maxI=requests.get('http://golem.fjfi.cvut.cz/shots/%i/Operation/Discharge/Stabilization/Results/Max_IexStab'%shot)
MaxIstab.append(url_maxI.text)
data={'Plasma duration': tdur, 'Stabilization parameters': StabParam, 'Max current in stab': MaxIstab, 'Shot_no': ScanDef}
df_results=pd.DataFrame(data)
df_results=df_results.set_index('Shot_no')
df_results
Plasma duration | Stabilization parameters | Max current in stab | |
---|---|---|---|
Shot_no | |||
36015 | 12.071 | 5000,0;6000,0;9000,0;12000,0\n | 3.86 |
36016 | 10.667 | 5000,0;6000,20;9000,20;12000,0\n | 13.74 |
0 | 13.277 | 5000,0;6000,-20;9000,-20;12000,0\n | 9.27 |
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['I_exStab'].plot(grid=True, c='r')
ax4.set(xlim=(start,end), ylim=(-maxI,maxI),xlabel= 'Time [ms]', ylabel = 'IexStab [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')
R=0.063
# L=0.000156
# maxI=max(abs(df_processed['I_exStab']))
u=df_processed['U_exStab'].loc[:]
t=u.index
for i in range(len(u)):
if abs(u.iat[i])>0.9:
StabStart=t[i]
print('Stabilization start:', round(StabStart,2),'ms')
break
t_maxI=df_processed['I_exStab'].idxmax()
eth=maxI/np.e
i_eth=maxI-eth
after_start=df_processed['I_exStab'].loc[StabStart:t_maxI]
dt=(t[1]-t[0])*1e-3
tau=0
for value in after_start.values:
if abs(value) < abs(i_eth):
tau += dt
else:
break
L_exp=R*tau
# R_exp=L/tau
print('MaxI:',round(maxI,3),'A','\nI(tau):',round(i_eth,3),'A','\ntau:',round(tau*1e3,3),'ms',
'\nL:', round(L_exp*1e6,4),'uH','\nR:', R,'ohm')
Stabilization start: 7.11 ms MaxI: 9.265 A I(tau): 5.857 A tau: 0.42 ms L: 26.46 uH R: 0.063 ohm