import pandas as pd
import numpy as np
Osciloskop, který používáme je také vidět na stránce výboje a data se ukládají na
http://golem.fjfi.cvut.cz/shotdir/XXXX/Devices/Oscilloscopes/RigolMSO5104-a/
konkrétně jsou napění na závit, napětí indukované na mg. cívcem napětí na rogovského pásku a napění na fotodidě na:
shot_no = 45320
Data budeme stahovat pomocí fuknce pandas.read_csv
kromě odkazu na csv soubor s daty přidáme argumenty names = ['t','U'], index_col='t'
,
aby měli sloupce rozumné názvy a sloupec s časem se použíl jako index (což nám zjednoduší práci později)
# postupně načteme data z výboje
U_Loop = pd.read_csv('http://golem.fjfi.cvut.cz/shotdir/'+str(shot_no)+'/Devices/Oscilloscopes/RigolMSO5104-a/U_Loop.csv',
names = ['t','U'], index_col='t')
U_BtCoil = pd.read_csv('http://golem.fjfi.cvut.cz/shotdir/'+str(shot_no)+'/Devices/Oscilloscopes/RigolMSO5104-a/U_BtCoil.csv',
names = ['t','U'], index_col='t')
U_RogCoil = pd.read_csv('http://golem.fjfi.cvut.cz/shotdir/'+str(shot_no)+'/Devices/Oscilloscopes/RigolMSO5104-a/U_RogCoil.csv',
names = ['t','U'], index_col='t')
U_photod = pd.read_csv('http://golem.fjfi.cvut.cz/shotdir/'+str(shot_no)+'/Devices/Oscilloscopes/RigolMSO5104-a/U_photod.csv',
names = ['t','U'], index_col='t')
#takhle vypadá tabulka s daty napětí na jeden závit transformátoru
U_Loop
U | |
---|---|
t | |
0.00000 | 0.000000 |
0.00005 | 0.399716 |
0.00010 | 0.000000 |
0.00015 | -1.199148 |
0.00020 | 0.399716 |
... | ... |
0.04980 | -9.992900 |
0.04985 | -8.793752 |
0.04990 | -9.992900 |
0.04995 | -9.593184 |
NaN | NaN |
1001 rows × 1 columns
# jednoduše ho můžeme vykreslit pomocí funkce plot()
#U_Loop.plot()
U magnetických diagnostik už musíme provést nějaké zpracování dat, protože z obrázku toho moc nepoznáme.
#U_BtCoil.plot()
Podle zákona o elektromagnetické indukci měříme $$U = - \frac{\Delta \Phi}{\Delta t} = - \frac{\Phi_2 - \Phi_1}{t_2 - t_1}$$ kde $\Phi$ je magjetický indukční tok, a $t$.
Tahle rovnice platí pro každou naměřenou hodnotu, postupně můžeme odvodit, že hodnoty indukce spočítáme jako: $$\Phi_n = -\left( U_1 + U_2 + \cdots U_n \right)\cdot \Delta t$$
Tedy potřebujeme získat postupné součty z U_BtCoil
a U_RogCoil
, které pak vynásobíme časovým krokem.
K tomu použijeme fukci cumsum()
# Protože výsledkem funkce cumsum() je další pandas objekt, můžeme funkce řetězit a rovnou vykreslit grafy
#U_BtCoil.cumsum().plot()
#U_RogCoil.cumsum().plot()
Z předchozích grafů vidíme, že musíme opravit znaménko u U_RogCoil
a zbavit se malé chyby, která se ve výpočtu kumuluje (offset)
Z grafů vidíme že až do času 1 ms jsou data okolo 0, ale bohužel ne přesně 0. Vypočítáme si průměrnou hodnotu signálu do času 0.001 s a tu odečteme od uložených dat.
offset_BtCoil = U_BtCoil[:0.001].mean()
offset_RogCoil = U_RogCoil[:0.001].mean()
dt = U_Loop.index[1]-U_Loop.index[0]
Bt_calib = 849765.7342657342
Rog_coil_calib = 4424581.0055865925
Bt = (U_BtCoil-offset_BtCoil).cumsum()*dt*Bt_calib
Itot = (-U_RogCoil+offset_RogCoil).cumsum()*dt*Rog_coil_calib
Bt.plot()
Itot.plot()
#U_RogCoil_no_offset_sum.max()
<AxesSubplot:xlabel='t'>
Nakonec jak vypadají všechny diagnostiky v jednom grafu
R_chamber = 0.0095 #Ohm
e = 1.60217663e-19
n = np.loadtxt("http://golem.fjfi.cvut.cz/shots/"+str(shot_no)+"/Diagnostics/Interferometry/ne_lav_mean")
V = 0.08
U = U_Loop
Ip = Itot-U_Loop/R_chamber
T = 0.9*(U/Ip)**(-2/3)
conf_time = (e*n*T*V)/(3*U*Ip)*1e6 # us
import matplotlib.pyplot as plt
fig, (ax1,ax2,ax3,ax4) = plt.subplots(4,sharex=True, figsize = [6,8])
U_Loop.rename(columns={'U':'U_loop'}).plot(ax=ax1)
Bt.rename(columns={'U':'B_t'}).plot(ax=ax2)
Itot.rename(columns={'U':"I_tot"}).plot(ax=ax3)
Ip.rename(columns={'U':"I_p"}).plot(ax=ax3)
U_photod.rename(columns={'U':"U_photod"}).plot(ax=ax4)
plt.xlabel('t [s]')
plt.savefig('diag.png')
T.rename(columns={'U':"T_e"})[0.004:0.0125].plot()
plt.ylabel('T_e [eV]')
plt.xlabel('t [s]')
plt.savefig('T_e.png')
conf_time[0.004:0.0125].rename(columns={'U':r"$\tau_E$"}).plot()
plt.ylabel(r'$\tau_E$ [$\mu$s]')
plt.xlabel('t [s]')
plt.savefig('icon-fig.png')