GOLEM Example: Electron Energy Confinement Time $\tau_E$ Calculation¶

# First let's get the data from the DAS in the columns specified by the given identifiers.
IDuloop = open_data(VacuumShotNo,'uloop')
IDirog = open_data(VacuumShotNo,'irog')
time =  np.linspace(IDuloop['t_start'],IDuloop['t_end'], len(IDuloop['data']))

# Let's show raw data:
f,ax = plt.subplots(2,sharex=True)
f.suptitle('Vacuum shot: #'+str(VacuumShotNo)+ '- raw data from NIstandard DAS')

f.subplots_adjust(hspace=0.05)
ax[0].set_ylim(0,2)
ax[0].set_ylabel('channel #1 [V]')
ax[0].plot(time*1e3, IDuloop['data'],label='ID:uloop')
ax[1].set_ylim(-.1,0.3)
ax[1].set_ylabel('channel #3 [V]')
ax[1].plot(time*1e3, IDirog['data'],label='ID:irog')
ax[1].set_xlabel('Time [ms]')
[a.legend(loc='best') for a in ax]
plt.show()

# Uloop from IDuloop needs only multiplication
U_l = IDuloop['data']*UloopCalibration
# Chamber current I_ch needs offset correction, integration and calibration
i_offset = time.searchsorted(.004)
offset = np.mean(IDirog['data'][:i_offset]) # Up to 4 ms no signal, good zero specification
I_ch = cumtrapz(IDirog['data']-offset, time,initial=0)  # offset correction and integration
I_ch *= RogowskiCalibration # calibration

# And let's plot it
f,ax = plt.subplots(2,sharex=True)
f.subplots_adjust(hspace=0.05)
f.suptitle('Vacuum shot: #'+str(VacuumShotNo)+ '- real data')

ax[0].set_ylabel('$U_l$ [V]')
ax[0].plot(time*1e3,U_l,label='Loop voltage $U_l$')
ax[1].set_ylabel('$I_{ch}$ [A]')
ax[1].plot(time*1e3,I_ch,label='Chamber current $I_{ch}$')
ax[1].set_xlabel('Time [ms]')
[a.legend(loc='best') for a in ax]
plt.show()

t_start = time.searchsorted(0.014) 
t_end = time.searchsorted(0.015) 
# R_ch calculation via Ohm's law:
R_ch = (U_l[t_start:t_end]/I_ch[t_start:t_end]) # in Ohms
# final value
R_ch, R_err = np.mean(R_ch), np.std(R_ch)

print("Chamber resistance =  (%.2f +/- %.2f) mOhm"%(R_ch*1e3, R_err*1e3))

# Get data
IDuloop = open_data(PlasmaShotNo,'uloop')
IDirog = open_data(PlasmaShotNo,'irog')
#time =  np.linspace(IDuloop['t_start'],IDuloop['t_end'], len(IDuloop['data']))

IDelectron_density = open_data(PlasmaShotNo,'electron_density')
#interpolate the density on the same time axis as the rest of the quantities
n_e = IDelectron_density['data']
time_ne =  np.linspace(IDelectron_density['t_start'],IDelectron_density['t_end'], len(n_e))
n_e = np.interp(time, time_ne, n_e)

# Let's show raw data from NIstandard:
f,ax = plt.subplots(2,sharex=True)
f.subplots_adjust(hspace=0.1)
f.suptitle('Plasma shot: #'+str(PlasmaShotNo)+ '- raw data from NIstandard DAS')

ax[0].set_ylim(0,2)
ax[0].set_ylabel('channel #1 [V]')
ax[0].plot(time*1e3, IDuloop['data'],label='ID:uloop')
ax[1].set_ylim(0,0.3)
ax[1].set_ylabel('channel #3 [V]')
ax[1].plot(time*1e3, IDirog['data'],label='ID:irog')
ax[1].set_xlabel('Time [ms]')
[a.legend(loc='best') for a in ax]
plt.show()

# Let's show raw data from interferometer diagnostics :
f,ax = plt.subplots()
ax.set_title('Plasma shot: #'+str(PlasmaShotNo)+ '- raw data from TektronixDPO DAS')

ax.set_xlabel('Time [ms]')
ax.set_ylabel('$n_e\ \mathrm{[m^{-3}]}$')
ax.plot(time*1e3, n_e,label='ID:electron_density')
ax.legend() 
ax.set_ylim(0,None)
plt.show()

# Uloop from IDuloop needs only multiplication
U_l = IDuloop['data']*UloopCalibration
# Chamber current I_ch needs offset correction, integration and calibration
offset = np.mean(IDirog['data'][:i_offset]) # Up to 4 ms no signal, good zero specification
I_chANDp = cumtrapz(IDirog['data']-offset, time, initial=0)  # offset correction and integration
I_chANDp *= RogowskiCalibration # calibration
I_p=I_chANDp-U_l/R_ch # getting real plasma current from Ohm's law

# And let's plot it
f,ax = plt.subplots(3,sharex=True )
f.subplots_adjust(hspace=0.05)
f.suptitle(r'Plasma shot: #'+str(PlasmaShotNo)+ '- real data')

ax[0].set_ylabel('$U_l$ [V]')
ax[0].plot(time*1e3, U_l,label='Loop voltage $U_l$')
ax[1].set_ylabel('$I_{ch+p}$ [kA]')
ax[1].plot(time*1e3,I_chANDp/1e3,label='Chamber+plasma current $I_{ch+p}$')
ax[2].set_ylabel('$I_{p}$ [kA]')
ax[2].plot(time*1e3,I_p/1e3,label='Plasma current $I_{p}$')
ax[2].set_xlabel('Time [ms]')
[a.legend(loc='best') for a in ax]
plt.show()

#select only a signal with the plasma 
i_start,i_end = time.searchsorted((IDuloop['plasma_start']+5e-4,IDuloop['plasma_end']-5e-4))

n_e = n_e[i_start:i_end]
I_p = I_p[i_start:i_end]
U_l = U_l[i_start:i_end]
time = time[i_start:i_end]


# And let's plot it
f,ax = plt.subplots(3,sharex=True)
f.subplots_adjust(hspace=0.05)
f.suptitle('Plasma shot: #'+str(PlasmaShotNo)+ '- real data')

ax[0].set_ylabel('$U_l$ [V]')
ax[0].plot(time*1e3, U_l,label='Loop voltage $U_l$')
ax[1].set_ylabel('$I_{p}$ [kA]')
ax[1].plot(time*1e3,I_p/1e3,label='Plasma current $I_{p}$')
ax[2].set_ylabel('$n_{e}$ [$m^{-3}$]')
ax[2].plot(time*1e3,n_e,label='Electron density $n_{e}$')
ax[-1].set_xlabel('Time [ms]')
[a.legend(loc='best') for a in ax]
plt.show()

R_p = np.maximum(0,U_l/I_p) # plasma resistivity adopted to avoid power(xx,-2./3) problems @ T_e calculation
P_OH = U_l*I_p  # Plasma heating power
T_e = 0.9*R_p**(-2./3) # Electron temperature (the GOLEM specific case)
W_p = V_p*n_e*eV*T_e/3 # Energy content in the plasma
# Let's make final calculations in the stationary phase to avoid complex derivation calculations of plasma energy balance
imin,imax = time.searchsorted((0.016, 0.020))
P_OH = P_OH[imin:imax]
W_p = W_p[imin:imax]
time = time[imin:imax]
tau = W_p/P_OH 
#Plot
f,ax = plt.subplots(3, sharex=True)
f.suptitle('Plasma shot: #'+str(PlasmaShotNo)+ ' energy balance')

f.subplots_adjust(hspace=0.001)
ax[0].set_ylabel('$P_{OH}$ [kW]')
ax[0].plot(time*1e3, P_OH/1000,label='Ohming heating power $P_{OH}$')
ax[1].set_ylabel('$W_p$ [J]')
ax[1].plot(time*1e3,W_p,label='Plasma energy volume $W_p$')
ax[2].set_ylabel(r'$\tau_E$ [us]')
ax[2].plot(time*1e3,tau*1e6,label=r'Energy confinement time $\tau_E$')
ax[2].set_ylim(0,80)
ax[2].grid('on')
[a.legend(loc='best') for a in ax]
ax[-1].set_xlabel('Time [ms]')
plt.show()

# Final calculation
tau, tau_err = np.mean(tau), np.std(tau)
print("Energy confinement time =  (%.0f +/- %.0f) us"%(tau*1e6,tau_err*1e6))