%% 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'])