Plasma detection from $U_{loop}$

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 = 39316

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.0009999532282507018

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

9.89488848

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

-1.0

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: False, from 9.9 to -1.0 ms (lifetime -1.0 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)