This is a quick plasma start and end detection from the peaks in the derivative of the $U_{loop}$ signal. This procedure is run as soon as possible in order to provide the information on plasma existance to other routines.
import numpy as np
from scipy import signal # for peak detection
import pandas as pd
import requests
import subprocess #workarround for db operation
import os
shot_no = 39318
def update_db_current_shot_alchemy(field_name, value):
try:
# local UNIX socket connection
engine = sqlalchemy.create_engine('postgresql://golem@/golem_database?host=/var/run/postgresql')
#engine = sqlalchemy.create_engine('postgresql://golem:rabijosille@192.168.2.116/golem_database')
except:
return
engine.execute(f"""UPDATE shots SET "{field_name}"={value} WHERE shot_no IN(SELECT max(shot_no) FROM shots)""")
def update_db_current_shot(field_name, value):
os.system('export PGPASSWORD="rabijosille";psql -c "UPDATE shots SET '+field_name+'='+str(value)+'WHERE shot_no IN(SELECT max(shot_no) FROM shots)" -q -U golem golem_database')
os.makedirs('Results', exist_ok=True)
def save_scalar(phys_quant, value, format_str='%.3f'):
with open("Results/"+phys_quant, 'w') as f:
f.write(format_str % value)
update_db_current_shot(phys_quant,value)
t_CD = requests.get(
f'http://golem.fjfi.cvut.cz/shots/{shot_no}/Production/Parameters/Tcd')
t_CD = float(t_CD.content) * 1e-3 # from us to ms
t_CD
1.9000000000000001
U_loop = pd.read_csv(f'http://golem.fjfi.cvut.cz/shots/{shot_no}/Diagnostics/BasicDiagnostics/U_Loop.csv',
names=['Time', 'U_loop'], index_col='Time', squeeze=True)
The processing is done only beyond the CD system activation time (to prevent thyristor noise from e.g. the $B_t$ system)
U_loop.index *= 1e3 # s to ms
U_loop = U_loop.loc[t_CD:]
dt = np.diff(U_loop.index[[0,-1]]).item() / U_loop.index.size
dt
0.0009999532263797943
def to_points(time, must_be_odd=False):
n = int(np.rint(time / dt))
if n % 2 == 0 and must_be_odd:
n += 1
return n
Estimate $\frac{d U_{loop}}{dt}$ by a second order Savitzky-Golay 1. derivative filter.
Dt = 0.1 # window length in ms
dU = pd.Series(signal.savgol_filter(U_loop.loc[t_CD:], to_points(Dt, must_be_odd=True),
2, 1, delta=dt), index=U_loop.index, name='dU_dt')
def argrelmax_time(sig, Dt, res_index, threshold_value=None,
comparator=np.greater):
idxs = signal.argrelextrema(sig.values, comparator, order=to_points(Dt))
sig_sel = sig.iloc[idxs]
if threshold_value is not None:
sig_sel = sig_sel[comparator(sig_sel, threshold_value)]
return sig_sel.index[res_index]
The plasma starts when the $U_{loop}$ decreases significantly for the first (index 0) time.
try:
t_plasma_start = argrelmax_time(dU, Dt/2, 0, threshold_value=-3,
comparator=np.less)
except IndexError:
t_plasma_start = -1.0 # not found, nan convention
t_plasma_start
3.9481939200000005
The plasma ends when the $U_{loop}$ increases significantly for the last (index -1) time.
try:
t_plasma_end = argrelmax_time(dU.loc[t_plasma_start:], Dt/2, -1,
threshold_value=10,
comparator=np.greater)
except IndexError:
t_plasma_end = -1.0 # not found, nan convention
t_plasma_end
16.2251939
is_plasma = int(t_plasma_start > 0 and t_plasma_end > 0)
if is_plasma:
t_plasma_duration = t_plasma_end - t_plasma_start
else:
t_plasma_duration = -1.0 # convention instead of nan
print(f'plasma: {bool(is_plasma)}, from {t_plasma_start:.1f} to {t_plasma_end:.1f} ms (lifetime {t_plasma_duration:.1f} ms)')
plasma: True, from 3.9 to 16.2 ms (lifetime 12.3 ms)
save_scalar("is_plasma", is_plasma)
save_scalar("t_plasma_start", t_plasma_start)
save_scalar("t_plasma_end", t_plasma_end)
save_scalar("t_plasma_duration", t_plasma_duration)