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