Rail probe - simple data rendering¶
The aim of the script is to plot raw signals from individual probes: rail, Langmuir and ball-pen probe.
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.ticker as ticker
import pandas as pd
import io
import os
import requests
import scipy.stats
import scipy.signal
from scipy.signal import butter, filtfilt, lfilter
from scipy.signal import welch
from scipy import integrate
from scipy.optimize import curve_fit
%pylab notebook
Populating the interactive namespace from numpy and matplotlib
shot_no = 36277
adress_lp = "RailProbe/I_LP"
adress_rp = "RailProbe/I_RP"
adress_bpp = "RailProbe/U_fl_BPP"
adress_bias = "RailProbe/Null"
adress_lp_r = "r_lp_vachar"
adress_rp_r = "r_rp_vachar"
save_destination= "Results/"
The aim of the script is to plot raw signals from individual probes: rail, Langmuir and ball-pen probe.
def print_and_save(phys_quant, value, format_str='%.3f'):
print(phys_quant+" = %.5f" % value)
with open(save_destination+phys_quant, 'w') as f:
f.write(format_str % value)
def get_data(shot, channel):
url = "http://golem.fjfi.cvut.cz/shots/"+str(shot)+"/Diagnostics/"+channel+".csv"
s = requests.get(url).content
data = pd.read_csv(io.StringIO(s.decode('utf-8')), sep=',', header=0, names=["t","U"])
return data
def get_param(shot, param):
r = requests.get("http://golem.fjfi.cvut.cz/shots/"+str(shot)+"/Diagnostics/RailProbe/Parameters/"+str(param))
response = r.json()
return float(response)
def get_data_uloop(shot):
calibration_uloop = 5.23*0.1966
data_uloop = get_data(shot,"BasicDiagnostics/U_Loop")
data_uloop["U"] *= calibration_uloop
return data_uloop
def get_data_btor(shot):
calibration_btor = 69.3
data = get_data(shot,"BasicDiagnostics/U_BtCoil")
zerotime = int(0.03*len(data))
data_for_int = data["U"] - np.mean(data["U"][0:zerotime])
data_int = integrate.cumtrapz(data_for_int, data["t"], initial=0)
data_calibrated = calibration_btor*data_int
data_btor = pd.DataFrame({'t': data["t"], 'U': data_calibrated})
return data_btor
def get_data_irog(shot):
calibration_irog = 1.42*3320000/0.8765
resistance_ch = 0.0097/0.915
data_irog = get_data(shot,"BasicDiagnostics/U_RogCoil")
zerotime = int(0.05*len(data_irog))
data_for_int_irog = data_irog["U"] - np.mean(data_irog["U"][0:zerotime])
data_int_irog = integrate.cumtrapz(data_for_int_irog, data_irog["t"], initial=0)
data_irog_calibrated = data_int_irog * calibration_irog
data_uloop = get_data_uloop(shot)
data_irog_calibrated = -1*data_irog_calibrated - data_uloop["U"]/resistance_ch
data_irog = pd.DataFrame({'t': data_irog["t"], 'U': data_irog_calibrated})
return data_irog
def get_plot(shot_no, plot_start, plot_end):
fig, ax = plt.subplots(7, figsize=(8, 14), sharex=True, tight_layout=True, gridspec_kw = {'wspace':0, 'hspace':0})
#irog
ax[0].plot(get_data_irog(shot_no)["t"]*1e3, get_data_irog(shot_no)["U"], "b-",label="Plasma current", linewidth=0.5)
ax[0].set(ylabel="$I_{plasma}$ [A]", xlim=[plot_start,plot_end])
#uloop
ax[1].plot(get_data_uloop(shot_no)["t"]*1e3, get_data_uloop(shot_no)["U"], "b-",label="Loop voltage", linewidth=0.5)
ax[1].set(ylabel="$U_{loop}$ [V]", xlim=[plot_start,plot_end])
#btor
ax[2].plot(get_data_btor(shot_no)["t"]*1e3, get_data_btor(shot_no)["U"], "b-",label="Tor. mag. field", linewidth=0.5)
ax[2].set(ylabel="$B_{T}$ [T]", xlim=[plot_start,plot_end])
#lp
data_lp = get_data(shot_no, adress_lp)
ax[3].plot(data_lp["t"]*1e3, 1e3*data_lp["U"]/get_param(shot_no, adress_lp_r), "b-",label="Langmuir probe", linewidth=0.5)
ax[3].set(ylabel="$I_{lp}$ [mA]", xlim=[plot_start,plot_end])#, ylim=[-500,500])
#rp
data_rp = get_data(shot_no, adress_rp)
ax[4].plot(data_rp["t"]*1e3, 1e3*data_rp["U"]/get_param(shot_no, adress_rp_r), "r-", label="Rail probe", linewidth=0.5)
ax[4].set(ylabel="$I_{rp}$ [mA]", xlim=[plot_start,plot_end])#, ylim=[-500,500])
#bpp
data_bp = get_data(shot_no, adress_bpp)
ax[5].plot(data_bp["t"]*1e3, data_bp["U"], "b-",label="Ball-pen probe", linewidth=0.5)
ax[5].set(ylabel="$U_{bpp}$ [V]", xlim=[plot_start,plot_end])#, ylim=[-0.4,0.15])
#bias
data_bias = get_data(shot_no, adress_bias)
ax[6].plot(data_bias["t"]*1e3, data_bias["U"]*100, "b-", label="Bias voltage", linewidth=0.5)
ax[6].set(xlabel="t [ms]", ylabel="$U_{bias}$ [V]", xlim=[plot_start,plot_end])#, ylim=[-0.05,0.75])
for a in ax:
a.legend(loc=2)
a.grid()
a.xaxis.set_major_locator(ticker.MultipleLocator(2))
plt.show()
plot_start, plot_end = 0, 0.02*1e3
get_plot(shot_no, plot_start, plot_end)
"""
def sin_wave(t,A,f,phi,B):
return A*np.sin(2*np.pi*f*t + phi) + B
def remove_parasitic(data):
popt_lp, pcov_lp = curve_fit(sin_wave, data["t"][0:40000]*1e3, data["U"][0:40000], p0 = [9.43925868e-02, 5, 6.22758426e-01,0])
data["U"] = data["U"] - sin_wave(data["t"]*1e3,*popt_lp)
return data
"""
'\ndef sin_wave(t,A,f,phi,B):\n return A*np.sin(2*np.pi*f*t + phi) + B\n\ndef remove_parasitic(data):\n popt_lp, pcov_lp = curve_fit(sin_wave, data["t"][0:40000]*1e3, data["U"][0:40000], p0 = [9.43925868e-02, 5, 6.22758426e-01,0])\n data["U"] = data["U"] - sin_wave(data["t"]*1e3,*popt_lp)\n return data\n'
def linear(x,a,b):
return a*x + b
def lowpass_filter(signal, cutoff, order):
sample_freq = 1/(signal["t"][1]-signal["t"][0])
nyq = 0.5 * sample_freq
normal_cutoff = cutoff / nyq
b, a = butter(order, normal_cutoff, btype='low', analog=False)
signal["U"] = filtfilt(b, a, signal["U"])
return signal
def remove_parasitic(probe, bias_input):
order = 5
cutoff = 35000
bias = bias_input.copy(deep=True)
bias["U"] = np.gradient(bias["U"])
paras_cur_lowpass = lowpass_filter(bias, cutoff, order)
popt, pcov = curve_fit(linear, paras_cur_lowpass["U"][0:10000], probe["U"][0:10000])
probe["U"] = probe["U"] - linear(paras_cur_lowpass["U"],*popt)
return probe
fig, ax = plt.subplots(3, figsize=(8, 6), sharex=True, tight_layout=True, gridspec_kw = {'wspace':0, 'hspace':0})
data_bias = get_data(shot_no, adress_bias) ### BIAS
data_lp = get_data(shot_no, adress_lp)
data_rp = get_data(shot_no, adress_rp)
data_lp["U"] = 1e3*data_lp["U"]/get_param(shot_no, adress_lp_r)
data_rp["U"] = 1e3*data_rp["U"]/get_param(shot_no, adress_rp_r)
#lp
data_lp_clean = remove_parasitic(data_lp, data_bias)
ax[0].plot(data_lp_clean["t"]*1e3, data_lp_clean["U"], "b-",label="Langmuir probe", linewidth=0.5)
ax[0].set(ylabel="$U_{lp}$ [V]", xlim=[plot_start,plot_end])#, ylim=[-0.005,0.075])
#rp
data_rp_clean = remove_parasitic(data_rp, data_bias)
ax[1].plot(data_rp_clean["t"]*1e3, data_rp_clean["U"], "b-",label="Rail probe", linewidth=0.5)
ax[1].set(ylabel="$U_{rp}$ [V]", xlim=[plot_start,plot_end])#, ylim=[-0.005,0.075])
#bias
ax[2].plot(data_bias["t"]*1e3, data_bias["U"], "b-", label="Bias voltage", linewidth=0.5)
ax[2].set(ylabel="$U_{bias}$ [V]", xlim=[plot_start,plot_end])#, ylim=[-0.05,0.75])
for a in ax:
a.legend(loc=2)
a.grid()
a.xaxis.set_major_locator(ticker.MultipleLocator(2))
plt.show()
def LangmuirFormula_3(U, I_sat, U_fl, T_e):
return I_sat * (1 - np.exp((U - U_fl)/T_e))
def LangmuirFormula_4(U, I_sat, U_fl, T_e, delta):
return I_sat - I_sat*np.exp((U - U_fl)/T_e) + (delta)*(U - U_fl)
def sweeping(probe, bias, left, right):
probe = remove_parasitic(probe, bias)
probe_volt = list(probe.loc[(probe['t'] >= left) & (probe['t'] <= right)]["U"])
bias_volt = list(bias.loc[(bias['t'] >= left) & (bias['t'] <= right)]["U"])
data_table = pd.DataFrame({'U': bias_volt, 'I': probe_volt})
return data_table
left = 0.0120
right = 0.0140
data_lp = get_data(shot_no, adress_lp)
data_lp["U"] = 1e3*data_lp["U"]/get_param(shot_no, adress_lp_r)
data_rp = get_data(shot_no, adress_rp)
data_rp["U"] = 1e3*data_rp["U"]/get_param(shot_no, adress_rp_r)
data_bias = get_data(shot_no, adress_bias)
data_bias["U"] = data_bias["U"]*100
#lp
data_lp_swept = sweeping(data_lp, data_bias, left, right)
data_lp_swept_cut = data_lp_swept.loc[data_lp_swept['I'] >= -300]
popt_lp, pcov_lp = curve_fit(LangmuirFormula_4, data_lp_swept_cut["U"], data_lp_swept_cut["I"], p0 = [0.09,-31,19,0], bounds=([0,-100,0,-1],[100,0,100,0]))
perr_lp = np.sqrt(np.diag(pcov_lp))
#rp
data_rp_swept = sweeping(data_rp, data_bias, left, right)
data_rp_swept_cut = data_rp_swept.loc[data_rp_swept['I'] >= -300]
popt_rp, pcov_rp = curve_fit(LangmuirFormula_4, data_rp_swept_cut["U"], data_rp_swept_cut["I"], p0 = [0.09,-31,19,0], bounds=([0,-100,0,-1],[100,0,100,0]))
perr_rp = np.sqrt(np.diag(pcov_rp))
names = ["I_sat","U_fl","T_e","d"]
for i,name in enumerate(names):
print_and_save(name+"_lp",popt_lp[i])
print_and_save(name+"_rp",popt_rp[i])
I_sat_lp = 29.29672 I_sat_rp = 100.00000 U_fl_lp = -30.88161 U_fl_rp = -26.43021 T_e_lp = 32.52861 T_e_rp = 39.73661 d_lp = -0.00000 d_rp = -0.56342
piv = np.linspace(-150,20,500)
fig, ax = plt.subplots(2, figsize=(8, 8), sharex=True, tight_layout=True, gridspec_kw = {'wspace':0, 'hspace':0.1})
ax[0].errorbar(data_lp_swept["U"], data_lp_swept["I"], fmt="k.", capsize=5, elinewidth=1, markersize=0.5, label="Data", zorder=1)
ax[0].plot(piv, LangmuirFormula_4(piv,*popt_lp), "r-", linewidth=1.5, label="4-parametr fit", zorder=2)#" $T_e$ = "+str(round(popt_lang[2],1))+" eV, $J_{sat}$ = "+str(round(popt_lang[0]/A,1))+" mA, $U_{fl}$ = "+str(round(popt_lang[1],1))+" V")#label=lab,
ax[0].legend(loc=3, title="Langmuir probe")
ax[1].errorbar(data_rp_swept["U"], data_rp_swept["I"], fmt="k.", capsize=5, elinewidth=1, markersize=0.5, label="Data", zorder=1)
ax[1].plot(piv, LangmuirFormula_4(piv,*popt_rp), "r-", linewidth=1.5, label="4-parametr fit", zorder=2)#" $T_e$ = "+str(round(popt_lang[2],1))+" eV, $J_{sat}$ = "+str(round(popt_lang[0]/A,1))+" mA, $U_{fl}$ = "+str(round(popt_lang[1],1))+" V")#label=lab,
ax[1].legend(loc=3, title="Rail probe")
for a in ax:
a.hlines(0,-150,50, linewidth=1)
a.vlines(0,-5,1, linewidth=1)
a.grid()
a.set(xlabel="$U_{probe}$ [V]", ylabel="$I$ [mA]", xlim=[-100,12], ylim=[-600,200])
a.xaxis.set_major_locator(ticker.MultipleLocator(10))
plt.savefig('icon-fig')
plt.show()
<ipython-input-1-ddb5bff661f7>:20: UserWarning: This figure includes Axes that are not compatible with tight_layout, so results might be incorrect.
plt.savefig('icon-fig')
fig, ax = plt.subplots(8, figsize=(8, 16), sharex=True, tight_layout=True, gridspec_kw = {'wspace':0, 'hspace':0})
time = np.arange(0,max(data_lp["t"]),0.001)
index = [0,1,2,3]
names = ["$I_{sat}$","$U_{fl}$","$T_e$","$d$"]
units = ["[mA]","[V]","[eV]","[A/V]"]
lower_bounds = [0,-40,0,-0.002]
higher_bounds = [0.2,0,30,0]
for left, right in zip(time[:-1],time[1:]):
t = 1e3*(right + left)/2
data_lp_swept = sweeping(data_lp, data_bias, left, right)
data_lp_swept_cut = data_lp_swept.loc[data_lp_swept['I'] >= -0.3]
popt_lp, pcov_lp = curve_fit(LangmuirFormula_4, data_lp_swept_cut["U"], data_lp_swept_cut["I"], p0 = [0.09,-31,19,0], bounds=([0,-100,0,-1],[100,0,100,0]))
perr_lp = np.sqrt(np.diag(pcov_lp))
for i, name, unit, lower, higher in zip(index, names, units, lower_bounds, higher_bounds):
ax[i].set(xlabel="$t$ [ms]", ylabel=name+" "+unit, xlim=[0,20], ylim=[lower,higher])
ax[i].grid()
ax[i].errorbar(t, popt_lp[i], yerr=perr_lp[i], fmt="b.", capsize=2, elinewidth=1, markersize=1)
ax[i].text(0.15, 0.9, "Langmuir probe: "+name, horizontalalignment='center', verticalalignment='center', transform=ax[i].transAxes, fontsize=10)
data_rp_swept = sweeping(data_rp, data_bias, left, right)
data_rp_swept_cut = data_rp_swept.loc[data_rp_swept['I'] >= -0.3]
popt_rp, pcov_rp = curve_fit(LangmuirFormula_4, data_rp_swept_cut["U"], data_rp_swept_cut["I"], p0 = [0.09,-31,19,0], bounds=([0,-100,0,-1],[100,0,100,0]))
perr_rp = np.sqrt(np.diag(pcov_rp))
for i, name, unit, lower, higher in zip(index, names, units, lower_bounds, higher_bounds):
ax[i+4].set(xlabel="$t$ [ms]", ylabel=name+" "+unit, xlim=[0,20], ylim=[lower,higher])
ax[i+4].grid()
ax[i+4].errorbar(t, popt_lp[i], yerr=perr_lp[i], fmt="g.", capsize=2, elinewidth=1, markersize=1)
ax[i+4].text(0.15, 0.9, "Rail probe: "+name, horizontalalignment='center', verticalalignment='center', transform=ax[i+4].transAxes, fontsize=10)
plt.show()
/opt/anaconda3/lib/python3.8/site-packages/scipy/optimize/minpack.py:828: OptimizeWarning: Covariance of the parameters could not be estimated
warnings.warn('Covariance of the parameters could not be estimated',
/opt/anaconda3/lib/python3.8/site-packages/scipy/optimize/minpack.py:828: OptimizeWarning: Covariance of the parameters could not be estimated
warnings.warn('Covariance of the parameters could not be estimated',
/opt/anaconda3/lib/python3.8/site-packages/scipy/optimize/minpack.py:828: OptimizeWarning: Covariance of the parameters could not be estimated
warnings.warn('Covariance of the parameters could not be estimated',
/opt/anaconda3/lib/python3.8/site-packages/scipy/optimize/minpack.py:828: OptimizeWarning: Covariance of the parameters could not be estimated
warnings.warn('Covariance of the parameters could not be estimated',