In [1]:
import os
#delete files from the previous discharge
if os.path.exists('icon-fig.png'):
os.remove('icon-fig.png')
import numpy as np
import matplotlib.pyplot as plt
import holoviews as hv
hv.extension('bokeh')
import hvplot.pandas
from scipy import integrate, signal, interpolate
import pandas as pd
import requests
from IPython.display import Markdown
Set discharge number¶
To calculate the plasma position, $B_t$ contribution must be subtracted from the Mirnov coil signal. This contribution can be acquired by two ways: from a dedicated vacuum shot with the same toroidal magnetic field, or it can be estimated from the standard diagnostics. If a suitable vacuum shot is available, its number is stored in the variable vacuum_shot. If it isn't available, the vacuum_shot variable shall be set to False and the $B_t$ contribution shall be estimated.
In [2]:
# data_URL = "http://golem.fjfi.cvut.cz/shots/{shot_no}/DASs/LimiterMirnovCoils/{identifier}.csv" #Mirnov coils and quadrupole
BDdata_URL = "http://golem.fjfi.cvut.cz/shots/{shot_no}/Diagnostics/BasicDiagnostics/{identifier}.csv" #BD = basic diagnostic
# data_URL = "http://golem.fjfi.cvut.cz/shots/{shot_no}/DASs/022Daniela0InnerStabilization/{identifier}.csv"
data_URL = "http://golem.fjfi.cvut.cz/shots/{shot_no}/Diagnostics/LimiterMirnovCoils/DAS_raw_data_dir/Nidatap_6133.lvm" #NI Standart (DAS)
from urllib.request import urlopen
def get_file(url, shot, silent=False):
URL = 'http://golem.fjfi.cvut.cz/shots/%i/%s' % (shot, url)
if not silent:
print(URL)
f = urlopen(URL)
try:
return np.loadtxt(f)
except ValueError: # the data couldn't be converted to a row of numbers
print('The data could not be read! Using default values.')
return np.array([35022])
shot_no = 35855
#shot_no = 35021
vacuum_shot = int(get_file('/Diagnostics/LimiterMirnovCoils/Parameters/vacuum_shot', 0))
#vacuum_shot = 35022 #number of the vacuum shot or 'False'
print('Vacuum shot: %i' % vacuum_shot)
ds = np.DataSource(destpath='')
Data access¶
In [3]:
def open_remote(shot_no, identifier, url_template=data_URL):
return ds.open(url_template.format(shot_no=shot_no, identifier=identifier))
def read_signal(shot_no, identifier, url = data_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)
def get_basic_diag_data(shot):
url = ( 'http://golem.fjfi.cvut.cz/shots/{}/Diagnostics/BasicDiagnostics/'
'basig_diagnostics_processed.csv'.format(shot) )
print(url)
return pd.read_csv(url, header=0, index_col='Time')
Mirnov coil signal processing: integration and $B_t$ elimination¶
Load the data of a given Mirnov coil, integrate it and subtract the $B_t$ contribution from them.
In [4]:
def elimination (shot_no, identifier, vacuum_shot = False):
#load data
mc = (read_signal(shot_no, identifier, data_URL, 'lvm'))
mc = mc.set_index('Time')
mc = mc.replace([np.inf, -np.inf, np.nan], value = 0)
mc = mc[identifier]
konst = 1/(3.8e-03)
if vacuum_shot == False:
signal_start = mc.index[0]
length = len(mc)
basic_diag = get_basic_diag_data(shot_no)
Bt = basic_diag['Bt']
if len(Bt)>len(mc):
Bt = Bt.iloc[:length]
if len(Bt)<len(mc):
mc = mc.iloc[:len(Bt)]
if identifier == 'mc1':
k=300
elif identifier == 'mc5':
k= 14
elif identifier == 'mc9':
k = 31
elif identifier == 'mc13':
k = -100
mc_vacuum = Bt/k
else:
mc_vacuum = (read_signal(vacuum_shot, identifier, data_URL, 'lvm'))
mc_vacuum = mc_vacuum.set_index('Time')
mc_vacuum = mc_vacuum[identifier]
mc_vacuum -= mc_vacuum.loc[:0.9e-3].mean() #remove offset
mc_vacuum = mc_vacuum.replace([np.inf, -np.inf], value = 0)
mc_vacuum = pd.Series(integrate.cumtrapz(mc_vacuum, x=mc_vacuum.index, initial=0) * konst,
index=mc_vacuum.index*1000, name= identifier) #integration
mc -= mc.loc[:0.9e-3].mean() #remove offset
mc = pd.Series(integrate.cumtrapz(mc, x=mc.index, initial=0) * konst,
index=mc.index*1000, name= identifier) #integration
#Bt elimination
mc_vacuum = np.array(mc_vacuum)
mc_elim = mc - mc_vacuum
return mc_elim
Find plasma start and end¶
In [5]:
plasma_start = float(get_file('Diagnostics/BasicDiagnostics/Results/t_plasma_start', shot_no))
plasma_end = float(get_file('Diagnostics/BasicDiagnostics/Results/t_plasma_end', shot_no))
print ('Plasma start =', round(plasma_start, 3), 'ms')
print ('Plasma end =', round(plasma_end, 3), 'ms')
Horizontal plasma position $\Delta r$ calculation¶
In [6]:
def horpol(shot_no, vacuum_shot=False):
mc1 = elimination(shot_no, 'mc1', vacuum_shot)
mc9 = elimination (shot_no, 'mc9', vacuum_shot)
b = 93
r = ((mc1-mc9)/(mc1+mc9))*b
r = r.replace([np.nan], value = 0)
return r.loc[plasma_start:plasma_end]
r = horpol(shot_no, vacuum_shot)
ax = r.plot()
ax.set(ylim=(-85,85), xlim=(plasma_start,plasma_end), xlabel= 'Time [ms]', ylabel = '$\Delta$r [mm]',
title = 'Horizontal plasma position #{}'.format(shot_no))
ax.axhline(y=0, color='k', ls='--', lw=1, alpha=0.4)
Out[6]:
Vertical plasma position $\Delta z$ calculation¶
In [7]:
def vertpol(shot_no, vacuum_shot = False):
mc5 = elimination(shot_no, 'mc5', vacuum_shot)
mc13 = elimination (shot_no, 'mc13', vacuum_shot)
b = 93
z = ((mc5-mc13)/(mc5+mc13))*b
z = z.replace([np.nan], value = 0)
return z.loc[plasma_start:plasma_end]
z = vertpol (shot_no, vacuum_shot)
ax = z.plot()
ax.set(ylim=(-85, 85), xlim=(plasma_start,plasma_end), xlabel= 'Time [ms]', ylabel = '$\Delta$z [mm]', title = 'Vertical plasma position #{}'.format(shot_no))
ax.axhline(y=0, color='k', ls='--', lw=1, alpha=0.4)
Out[7]:
Plasma column radius $a$ calculation¶
In [8]:
def plasma_radius(shot_no, vacuum_shot=False):
r = horpol(shot_no, vacuum_shot)
z = vertpol(shot_no, vacuum_shot)
a0 = 85
a = a0 - np.sqrt((r**2)+(z**2))
a = a.replace([np.nan], value = 0)
return a.loc[plasma_start:plasma_end]
a = plasma_radius(shot_no,vacuum_shot)
ax = a.plot()
ax.set(ylim=(0,85), xlim=(plasma_start,plasma_end), xlabel= 'Time [ms]', ylabel = '$a$ [mm]',
title = 'Plasma column radius #{}'.format(shot_no))
Out[8]:
Save the results¶
In [9]:
plasma_time = []
t = 0
for i in a:
if i>85 or i <0:
a = a.replace(i, value = 0)
else:
plasma_time.append(a.index[t])
t+=1
start = plasma_time[0]-1e-03
end = plasma_time[-1]-1e-03
print('start =', round(start, 3), 'ms')
print('end =', round(end, 3), 'ms')
In [10]:
r_cut = r.loc[start:end]
a_cut = a.loc[start:end]
z_cut = z.loc[start:end]
df_processed = pd.concat(
[r_cut.rename('r'), z_cut.rename('z'), a_cut.rename('a')], axis= 'columns')
df_processed
Out[10]:
In [11]:
savedata = 'plasma_position.csv'
df_processed.to_csv(savedata)
Markdown("[Plasma position data - $r, z, a$ ](./{})".format(savedata))
Out[11]:
Plot the results¶
In [12]:
hline = hv.HLine(0)
hline.opts(
color='k',
line_dash='dashed',
alpha = 0.4,
line_width=1.0)
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]'))
plot
Out[12]:
In [13]:
responsestab = requests.get("http://golem.fjfi.cvut.cz/shots/%i/Production/Parameters/PowSup4InnerStab"%shot_no)
stab = responsestab.text
if '-1' in stab or 'Not Found' in stab:
fig, axs = plt.subplots(3, 1, sharex=True, dpi=200)
r.plot(grid=True, ax=axs[0])
z.plot(grid=True, ax=axs[1])
a.plot(grid=True, ax=axs[2])
axs[2].set(ylim=(0,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))
else:
fig, axs = plt.subplots(4, 1, sharex=True, dpi=200)
dataNI = (read_signal(shot_no, 'Rog', data_URL, 'lvm')) #data from NI
dataNI = dataNI.set_index('Time')
dataNI = dataNI.set_index(dataNI.index*1e3)
Irog = dataNI['RogQuadr']
# Irog=abs(Irog)
Irog*=1/5e-3
r.plot(grid=True, ax=axs[0])
z.plot(grid=True, ax=axs[1])
a.plot(grid=True, ax=axs[2])
Irog.plot(grid=True, ax=axs[3])
axs[3].set(xlim=(start,end), ylabel = 'I [A]')
axs[2].set(ylim=(0,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')
