Tokamak GOLEM Basic diagnostics

Prerequisities: function definitions

In [1]:
from numpy import *
from scipy.constants import pi,mu_0,k,eV
import sqlalchemy   # high-level library for SQL in Python


from matplotlib.pyplot import *
from scipy.optimize import leastsq
import time, datetime, re
from scipy.signal import medfilt, fftconvolve, convolve
#from urllib  import urlopen #python 2.7 
from urllib.request  import urlopen, urlretrieve #python 3.0 
import os
import pandas as pd
import matplotlib.pyplot as plt

baseURL = "http://golem.fjfi.cvut.cz/utils/data/" #remote access
Destination=''
#os.makedirs( Destination, exist_ok=True );

ShotNo = 31733


def PrintAndSave(PhysQuant, Value, Format ):
    print(PhysQuant+" = %.5f" % Value)
    with open(Destination+PhysQuant, 'w') as f:
        f.write("%.3f" % Value)
    #engine = sqlalchemy.create_engine('postgresql://golem:rabijosille@192.168.2.116/golem_database')
    #engine.execute("""UPDATE shots SET pressure=50 WHERE shot_no IN(SELECT max(shot_no) FROM shots)""")  
    
def read_value(ShotNo, identifier): 
    return np.loadtxt(urlopen(baseURL+ str(ShotNo) + '/' +identifier))

def read_signal(ShotNo, identifier): 
    if not os.path.isfile('/tmp/'+str(ShotNo)+'/'+identifier): # Just to avoid downloading already downloaded data
        os.makedirs( '/tmp/'+str(ShotNo), exist_ok=True );
        urlretrieve(baseURL+ str(ShotNo) + '/' +identifier, '/tmp/'+str(ShotNo)+'/'+identifier )
    return pd.read_csv('/tmp/'+str(ShotNo)+'/'+identifier, skiprows=10, names=['Time',identifier])

$U_l$ management

Check the data availability

In [2]:
identifier='LoopVoltageCoil_raw';LoopVoltage = read_signal(ShotNo,identifier)
LoopVoltage.plot(x="Time",y=identifier,legend=False,grid=True).set(xlabel="Time [s]", ylabel="$U_l$ [V]", title="Loop voltage $U_l$ #"+str(ShotNo));plt.show()

Getting $U_l$ mean

In [3]:
UloopMean = mean(abs(LoopVoltage[identifier]));
PrintAndSave("UloopMean", UloopMean, "%.3f")

UloopMean = 8.03863

Getting $U_l$ max

In [4]:
from scipy.stats.mstats import mquantiles
UloopMax = mquantiles(LoopVoltage[identifier],0.99)[0];
PrintAndSave("UloopMax", UloopMax, "%.3f")

UloopMax = 19.04797

$B_t$ management

Check the data availability

In [5]:
identifier='BtCoil_raw';dBt = read_signal(ShotNo,identifier) # Magnetic measurement, we get dBt/dt
dBt.plot(x="Time",y=identifier, legend=False,grid=True).set(xlabel="Time [s]", ylabel="$dU_{B_t}/dt$ [V]", title="BtCoil_raw signal #"+str(ShotNo));plt.show()

Integration (it is a magnetic diagnostics) & calibration

In [6]:
from scipy import integrate

K_BtCoil=read_value(ShotNo,'K_BtCoil') # Get BtCoil calibration factor
print('BtCoil calibration factor K_BtCoil='+str( K_BtCoil)+' T/(Vs)')
Bt = pd.DataFrame(columns = ['Time', 'Bt'])
Bt['Time'] = dBt['Time']
Bt['Bt'] = integrate.cumtrapz(dBt['BtCoil_raw'], dBt['Time'], initial=0)*K_BtCoil
Bt.plot(x="Time",y="Bt", legend=False,grid=True).set(xlabel="Time [s]", ylabel="$B_t$ [T]", title="Toroidal magnetic field $B_t$ #"+str(ShotNo));
plt.show() 

BtCoil calibration factor K_BtCoil=70.42 T/(Vs)

Getting $B_t$ mean

In [7]:
BtMean = mean(abs(Bt['Bt']));
PrintAndSave("BtMean", BtMean, "%.3f")

BtMean = 0.24637

Getting $B_t$ max

In [8]:
BtMax = mquantiles(Bt['Bt'],0.99)[0];
PrintAndSave("BtMax", BtMax, "%.3f")

BtMax = 0.40969

Plasma + Chamber current $I_{p+ch}$ management

Check the data availability

In [9]:
identifier='RogowskiCoil_raw';dIpch = read_signal(ShotNo,identifier) # Magnetic measurement, we get dI_{p+ch}/dt
dIpch.plot(x="Time",y=identifier, legend=False,grid=True).set(xlabel="Time [s]", ylabel="$dU_{I_{p+ch}}/dt$ [V]", title="RogowskiCoil_raw signal #"+str(ShotNo));plt.show()

Integration (it is a magnetic diagnostics) & calibration

In [10]:
K_RogowskiCoil=read_value(ShotNo,'K_RogowskiCoil') # Get RogowskiCoil calibration factor
print('RogowskiCoil calibration factor K_RogowskiCoil='+str( K_RogowskiCoil)+' A/(Vs)')
Ipch = pd.DataFrame(columns = ['Time', 'Ipch'])
Ipch['Time'] = dIpch['Time']
Ipch['Ipch'] = abs(integrate.cumtrapz(dIpch['RogowskiCoil_raw'], dIpch['Time'], initial=0)*K_RogowskiCoil)
Ipch.plot(x="Time",y="Ipch", legend=False,grid=True).set(xlabel="Time [s]", ylabel="$I_{p+ch}$ [A]", title="Total (plasma+chamber) current $I_{p+ch}$ #"+str(ShotNo));
plt.show() 

RogowskiCoil calibration factor K_RogowskiCoil=5300000.0 A/(Vs)

Plasma current $I_{p}$ management

In [11]:
R_chamber=read_value(ShotNo,'R_chamber') # Get Chamber resistivity
print('Chamber resistivity R_chamber='+str( R_chamber)+' Ohm')
Ip = pd.DataFrame(columns = ['Time', 'Ip'])
Ip['Time'] = Ipch['Time']
Ip['Ip'] = Ipch['Ipch']-LoopVoltage['LoopVoltageCoil_raw']/R_chamber
Ip.plot(x="Time",y="Ip", legend=False,grid=True).set(xlabel="Time [s]", ylabel="$I_{p}$ [A]", title="Plasma current $I_{p}$ #"+str(ShotNo));

Chamber resistivity R_chamber=0.0097 Ohm

Getting $I_p$ mean

In [12]:
IpMean = mean(abs(Ip['Ip']));
PrintAndSave("IpMean", IpMean, "%.3f")

IpMean = 394.43697

Getting $I_p$ max

In [13]:
IpMax = mquantiles(Ip['Ip'],0.99)[0];
PrintAndSave("IpMax", IpMax, "%.3f")

IpMax = 1681.12360

Plasma detect

In [14]:
## Odstrcils:
dIp = pd.DataFrame(columns = ['Time', 'signal'])
dIp['Time'] = Ip['Time']
dIp['signal'] = Ip['Ip'].diff()
dIp['signal'][dIp['Time'] < 0.0025] = 0 # Avoid thyristors

dIp.plot(x="Time",y="signal", legend=False,grid=True).set(xlabel="Time [s]", ylabel="", title="#"+str(ShotNo));

PlasmaStart=dIp.loc[dIp['signal'] > 2].iloc[0]['Time']
PrintAndSave("PlasmaStart", PlasmaStart, "%.3f")

PlasmaEnd=dIp.loc[dIp['signal'] < -6].iloc[0]['Time']
PrintAndSave("PlasmaEnd", PlasmaEnd, "%.3f")

PlasmaStart = 0.00409
PlasmaEnd = 0.00959
In [ ]: