Plasma column position

Adapted version of InnerStabilisation.ipynb by Daniela Kropáčková, Honza Stöckel and Vojtěch Svoboda.

This notebook loads Mirnov coil data and calculates the plasma column position from them.

Import libraries

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 = 35762
#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='')  
http://golem.fjfi.cvut.cz/shots/0//Diagnostics/LimiterMirnovCoils/Parameters/vacuum_shot
Vacuum shot: 35746

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')
http://golem.fjfi.cvut.cz/shots/35762/Diagnostics/BasicDiagnostics/Results/t_plasma_start
http://golem.fjfi.cvut.cz/shots/35762/Diagnostics/BasicDiagnostics/Results/t_plasma_end
Plasma start = 3.759 ms
Plasma end = 20.554 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]:
<matplotlib.lines.Line2D at 0x7fd7403b98d0>

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]:
<matplotlib.lines.Line2D at 0x7fd7366c22d0>

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]:
[(0, 85),
 Text(0, 0.5, '$a$ [mm]'),
 (3.759, 20.554),
 Text(0.5, 0, 'Time [ms]'),
 Text(0.5, 1.0, 'Plasma column radius #35762')]

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')
start = 3.758 ms
end = 20.553 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]:
r z a
Time
3.759 -30.271789 65.071979 13.231318
3.760 -29.679453 63.987338 14.464552
3.761 -29.320489 63.392648 15.155037
3.762 -29.010551 62.569817 16.031935
3.763 -28.682863 61.938728 16.742307
... ... ... ...
20.548 -41.426466 55.788420 15.512592
20.549 -42.336846 57.106920 13.911261
20.550 -42.478190 58.251452 12.905421
20.551 -40.186198 58.555080 13.981496
20.552 -38.703509 59.243473 14.234538

16794 rows × 3 columns

In [11]:
savedata = 'plasma_position.csv'
df_processed.to_csv(savedata)

Markdown("[Plasma position data - $r, z, a$ ](./{})".format(savedata))

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')