Tokamak GOLEM - electron energy confinement time 
calculation
Strategic flowchart
Initialization
First we clear matlab workspace and define fundamental constants and parameters
clc; clear; close all
VacuumShotNo = 22475; % Vacuum shot number
PlasmaShotNo=22471; % Plasma shot number
baseURL='http://golem.fjfi.cvut.cz/utils/data/';
RogowskiCalibration=5.3e6; % Calibration for Rogowski coil [A/V(10us)] 03/2010, CASTOR calibration (Mr. Zacek), 10 us sampling, Mr. Zacek states 5.3e6 A/Vs
UloopCalibration=5.5; % Calibration for U_loop measurement
Mu0=4*pi*1e-7; % [H/m] Permeability
kB=1.3806503e-23; % [J/K] Botlzmann constant
eV=1.60217646e-19; % [J/eV] 1 eV equivalent to J
MajorRadius=0.4; % [m]
LimiterPosition=0.085; % [m] max plasma radius 0.05 for ring with Mach probes
MeanPlasmaRadius=0.060; % [m] mean plasma radius without any stabilization
PlasmaVolume=pi*LimiterPosition^2*2*pi*MajorRadius; % [m^3]
Chamber resistence 
calculation from vacuum shot
Plot raw data from NIstandard DAS (Data Acquisition System)
First let's get the data from the DAS in the columns specified by the given identifiers
Ul = str2num(webread([baseURL num2str(VacuumShotNo) '/uloop']));
Ir = str2num(webread([baseURL num2str(VacuumShotNo) '/irog']));
Let's show raw data
figure(1);
h1 = subplot(2,1,1);
plot(Ul(:,1)*1e3,Ul(:,2))
ylabel('channel #1 [V]'); set(gca,'XTickLabels',[]); ylim([0 2]); xlim([0 40])
legend('ID:uloop')
title(['Vacuum shot: #' num2str(VacuumShotNo) ' - raw data from NIstandard DAS'])
h2 = subplot(2,1,2);
p = get(h1, 'pos'); p(2) = p(2) - p(4)*1.1; set(h2,'pos',p) % proper positon of subplot
plot(Ir(:,1)*1e3,Ir(:,2))
ylabel('channel #3 [V]'); xlabel('Time [ms]'); ylim([-.1 0.3]); xlim([0 40])
legend('ID:irog')
Get real physical quantities from raw data
Physical quantities can be obtained by multiplying the raw data by appropriate calibration constants. In the case of the current, the derivative is measured, so the signal has to be numerically integrated.
from ID:uloop needs only multiplication
U_l = Ul(:,2)*UloopCalibration;
Chamber current 
needs offset correction, integration and calibration
offset = mean(Ir(Ir(:,1)<=.004,2)); % Up to 4 ms no signal, good zero specification
I_ch = cumtrapz(Ir(:,1),Ir(:,2)-offset);
I_ch = I_ch*RogowskiCalibration;
And let´s plot it
figure(2);
h1 = subplot(2,1,1);
plot(Ul(:,1)*1e3,U_l)
ylabel('U_l [V]'); set(gca,'XTickLabels',[]); ylim([-4 12]); xlim([0 40])
legend('Loop voltage')
title(['Vacuum shot: #' num2str(VacuumShotNo) ' - real data'])
h2 = subplot(2,1,2);
p = get(h1, 'pos'); p(2) = p(2) - p(4)*1.1; set(h2,'pos',p) % proper positon of subplot
plot(Ir(:,1)*1e3,I_ch)
ylabel('I_{ch} [A]'); xlabel('Time [ms]'); ylim([-1e3 1.5e3]); xlim([0 40])
legend('Chamber current I_{ch}')
Final chamber resistance specification
The chanber resistnace should be calculated in the time span where the physical quantities are reasonably stationary.
t_start = find(Ul(:,1)>=0.014,1,'first');
t_end = find(Ul(:,1)>=0.015,1,'first');
Calculation via Ohm's law 
R_ch = U_l(t_start:t_end)./I_ch(t_start:t_end);
R_mean = mean(R_ch);
R_err = std(R_ch);
disp(['Chamber resistance = (' num2str(R_mean*1e3,'%.2f') ' +/- ' num2str(R_err*1e3,'%.2f') ') mOhm'])
Energy confinement time 
calculation from plasma shot
Plot raw data from NIstandard DAS and interferometer
Get data
Ul = str2num(webread([baseURL num2str(PlasmaShotNo) '/uloop']));
Ir = str2num(webread([baseURL num2str(PlasmaShotNo) '/irog']));
ne = str2num(webread([baseURL num2str(PlasmaShotNo) '/electron_density']));
Interpolate the density on the same time axis as the rest of the quantities
n_e = interp1(ne(:,1),ne(:,2),Ul(:,1),'pchip');
Let´s show raw data from NIstandard
figure(3);
h1 = subplot(2,1,1);
plot(Ul(:,1)*1e3,Ul(:,2))
ylabel('channel #1 [V]'); set(gca,'XTickLabels',[]); ylim([0 2]); xlim([0 40])
legend('ID:uloop')
title(['Plasma shot: #' num2str(PlasmaShotNo) ' - raw data from NIstandard DAS'])
h2 = subplot(2,1,2);
p = get(h1, 'pos'); p(2) = p(2) - p(4)*1.1; set(h2,'pos',p) % proper positon of subplot
plot(Ir(:,1)*1e3,Ir(:,2))
ylabel('channel #3 [V]'); xlabel('Time [ms]'); ylim([-.1 0.3]); xlim([0 40])
legend('ID:irog')
Let´s show raw data from interferometer diagnostics
figure(4)
plot(Ul(:,1)*1e3,n_e)
ylabel('n_e [m^{-3}]'); ylim([0 6e18]); xlim([0 40])
legend('ID:electron\_density')
title(['Plasma shot: #' num2str(PlasmaShotNo) ' - raw data from TektronixDPO'])
Get real physical data
from ID:uloop needs only multiplication
U_l = Ul(:,2)*UloopCalibration;
Plasma current 
needs offset correction, integration calibration and chamber current subtraction (getting from Ohm´s law)
offset = mean(Ir(Ir(:,1)<=.004,2));
I_chANDp = cumtrapz(Ir(:,1),Ir(:,2)-offset);
I_chANDp = I_chANDp*RogowskiCalibration;
I_p = I_chANDp-U_l/R_mean;
And let´s plot it
figure(5);
h1 = subplot(3,1,1);
plot(Ul(:,1)*1e3,U_l)
ylabel('U_l [V]'); set(gca,'XTickLabels',[]); ylim([-4 12]); xlim([0 40])
legend('Loop voltage')
title(['Plasma shot: #' num2str(PlasmaShotNo) ' - real data'])
h2 = subplot(3,1,2);
p = get(h1, 'pos'); p(2) = p(2) - p(4)*1.1; set(h2,'pos',p) % proper positon of subplot
plot(Ir(:,1)*1e3,I_chANDp/1e3)
ylabel('I_{ch} [kA]'); xlabel('Time [ms]'); ylim([-1 5]); xlim([0 40])
legend('Chamber + plasma current I_{ch+p}')
h3 = subplot(3,1,3);
p = get(h2, 'pos'); p(2) = p(2) - p(4)*1.1; set(h3,'pos',p) % proper positon of subplot
plot(Ir(:,1)*1e3,I_p/1e3)
ylabel('I_p [kA]'); xlabel('Time [ms]'); ylim([-1 5]); xlim([0 40])
legend('Plasma current I_p')
Data correlation ...
Select only a signal with plasma
tstart = str2num(webread([baseURL num2str(PlasmaShotNo) '/plasma_start']));
tend = str2num(webread([baseURL num2str(PlasmaShotNo) '/plasma_end']));
i_start = find(Ul(:,1)>=tstart,1,'first');
i_end = find(Ul(:,1)>=tend,1,'first');
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 = Ul(i_start:i_end,1);
And let´s plot it
figure(6);
h1 = subplot(3,1,1);
plot(time*1e3,U_l)
ylabel('U_l [V]'); set(gca,'XTickLabels',[]); ylim([-4 12]); xlim([tstart tend]*1e3)
legend('Loop voltage U_l')
title(['Plasma shot: #' num2str(PlasmaShotNo) ' - real data'])
h2 = subplot(3,1,2);
p = get(h1, 'pos'); p(2) = p(2) - p(4)*1.1; set(h2,'pos',p) % proper positon of subplot
plot(time*1e3,I_p/1e3)
ylabel('I_p [kA]'); set(gca,'XTickLabels',[]); ylim([-1 5]); xlim([tstart tend]*1e3)
legend('Plasma current I_p')
h3 = subplot(3,1,3);
p = get(h2, 'pos'); p(2) = p(2) - p(4)*1.1; set(h3,'pos',p) % proper positon of subplot
plot(time*1e3,n_e)
ylabel('n_e [m^{-3}]'); ylim([0 7e18]); xlabel('Time [ms]'); xlim([tstart tend]*1e3)
legend('Electron density n_e')
Final calculations
Plasma resistivity adopted to avoid 
problems @ T_e calculation
R_p = U_l./I_p;
R_p(R_p<0 & isnan(R_p)) = 0;
Plasma heating power 
P_OH = U_l.*I_p;
Central electron temperature 
(the GOLEM specific case)
T_e = 0.9.*R_p.^(-2/3);
Energy content in the plasma 
W_p = PlasmaVolume.*n_e.*eV.*T_e./3;
Let´s make final calculations in the stationary phase to avoid complex derivation calculations of plasma energy balance
imin = find(time>=0.015,1,'first');
imax = find(time>=0.020,1,'first');
P_OH = P_OH(imin:imax);
W_p = W_p(imin:imax);
time = time(imin:imax);
tau = W_p./P_OH;
Plot
figure(7);
h1 = subplot(3,1,1);
plot(time*1e3,P_OH/1e3)
ylabel('P_{OH} [kW]'); set(gca,'XTickLabels',[]); ylim([20 27]);
legend('Ohmic heating power P_{OH}')
title(['Plasma shot: #' num2str(PlasmaShotNo) ' - energy balance'])
h2 = subplot(3,1,2);
p = get(h1, 'pos'); p(2) = p(2) - p(4)*1.1; set(h2,'pos',p) % proper positon of subplot
plot(time*1e3,W_p)
ylabel('W_p [J]'); set(gca,'XTickLabels',[]); ylim([0.6 1.3]);
legend('Plasma energy W_p')
h3 = subplot(3,1,3);
p = get(h2, 'pos'); p(2) = p(2) - p(4)*1.1; set(h3,'pos',p) % proper positon of subplot
plot(time*1e3,tau*1e6)
ylabel('n_e [m^{-3}]'); ylim([25 50]); xlabel('Time [ms]')
legend('Energy confinement time \tau_E')
grid on
Final calculation
tau_mean = mean(tau);
tau_err = std(tau);
disp(['Energy confinement time = (' num2str(tau_mean*1e6,'%.0f') ' +/- ' num2str(R_err*1e6,'%.0f') ') us'])