%% Tokamak GOLEM - electron energy confinement time $\tau_E$ 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 $R_{ch}$ 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.
%
% $U_l$ from ID:uloop needs only multiplication
U_l = Ul(:,2)*UloopCalibration;
%%
% Chamber current $I_{ch}$ 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/I_{ch}$$
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 $\tau_E$ 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
%
% $U_l$ from ID:uloop needs only multiplication
U_l = Ul(:,2)*UloopCalibration;
%%
% Plasma current $I_p$ needs offset correction, integration calibration
% and chamber current $I_{ch}$ 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 $R_p^{-2/3}$ problems @ T_e calculation
R_p = U_l./I_p;
R_p(R_p<0 & isnan(R_p)) = 0;
%%
% Plasma heating power $P_{OH}$
P_OH = U_l.*I_p;
%%
% Central electron temperature $T_e$ (the GOLEM specific case)
T_e = 0.9.*R_p.^(-2/3);
%%
% Energy content in the plasma $W_p$
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'])
