clc
clear all
close all

ShotNo=22471;
baseURL='http://golem.fjfi.cvut.cz/utils/data/';

RogowskiCalibration=5.3e6;%   #[A/V(10us)] 03/2010, CASTOR calibration (Mr. F.Zacek), 10 us vzorkovani, p. Zacek uvadi 5.3 e6 A/Vs
UloopCalibration=5.5;
Mu0=4*pi*1e-7;
kB=1.3806503e-23; % [J/K]
eV=1.60217646e-19; %[J/eV]
RoomTemperature=300;%  #[K]

MajorRadius=0.4;
MinorRadius=0.1;
Aspect = MinorRadius/MajorRadius;

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;
%% Time span
identifier='plasma_start';
dataURL=strcat(baseURL,int2str(ShotNo),'/',identifier); 
urlwrite(dataURL,identifier);
data = load(identifier, '\t');
tstart=data;
identifier='plasma_end';
dataURL=strcat(baseURL,int2str(ShotNo),'/',identifier); 
urlwrite(dataURL,identifier);
data = load(identifier, '\t');
tend=data;


%% Loop voltage
identifier='loop_voltage';
dataURL=strcat(baseURL,int2str(ShotNo),'/',identifier); 
urlwrite(dataURL,identifier);
data = load(identifier, '\t');
Ul=data(:,2);
time=data(:,1);
figure
plot(time*1000, Ul, 'r-') ;
xlabel('Time [ms]')
ylabel('U_l [V]')
xlim([tstart tend].*1e3)

%% Plasma current
identifier='plasma_current';
dataURL=strcat(baseURL,int2str(ShotNo),'/',identifier); 
urlwrite(dataURL,identifier);
data = load(identifier, '\t');
Ip=data(:,2);
time=data(:,1);
figure
plot(time*1000, Ip*1e-3, 'r-') ;
xlabel('Time [ms]')
ylabel('I_p [kA]')
xlim([tstart tend].*1e3)

%% Plasma resistivity
Rp=Ul./Ip;
figure
plot(time*1000, Rp*1e3, 'r-') ;
xlabel('Time [ms]')
ylabel('R_p [m\Omega]')
xlim([tstart tend].*1e3)
set(gca,'yscale','log')

%% Central electron temperature
Te=0.9*Rp.^(-2/3);
figure
plot(time*1000, Te, 'r-') ;
xlabel('Time [ms]')
ylabel('T_e(0,t) [eV]')
xlim([tstart tend].*1e3)
set(gca,'yscale','log')

%% Electron density
identifier='electron_density';
dataURL=strcat(baseURL,int2str(ShotNo),'/',identifier); 
urlwrite(dataURL,identifier);
data = load(identifier, '\t');
ne=data(:,2);
tt=data(:,1);
figure
plot(tt*1000, ne, 'r-') ;
xlabel('Time [ms]')
ylabel('n_e [m^{-3}]')
xlim([tstart tend].*1e3)
set(gca,'yscale','log')

%% Total energy content
iii=[];
for kk=1:length(tt)
    iii(kk)=find(time>=tt(kk),1);
end
Wp=ne.*eV.*Te(iii).*PlasmaVolume/3;
figure
plot(tt*1000, Wp, 'r-') ;
xlabel('Time [ms]')
ylabel('W_p [J]')
xlim([tstart tend].*1e3)
set(gca,'yscale','log')

%% Plasma heating power
Poh=Ul.*Ip;

%% Loss power
dWp_dt=zeros(length(tt),1);
dWp_dt(1)=(Wp(2)-Wp(1))/(tt(2)-tt(1));
for kk=2:length(tt)-1
    dWp_dt(kk)=(Wp(kk+1)-Wp(kk-1))/(tt(kk+1)-tt(kk-1));
end
dWp_dt(length(tt))=(Wp(end)-Wp(end-1))/(tt(end)-tt(end-1));
figure
plot(tt*1000, dWp_dt, 'r-') ;
xlabel('Time [ms]')
ylabel('dW_p/dt [W]')
xlim([tstart tend].*1e3)
set(gca,'yscale','log')

Ploss=Poh(iii)-dWp_dt;
figure
plot(time*1000, Poh, 'r-') ;
hold on
plot(tt*1000, Ploss, 'b-') ;
xlabel('Time [ms]')
ylabel('P [W]')
xlim([tstart tend].*1e3)
set(gca,'yscale','log')
legend('P_{OH}','P_{loss}')


%% Energy confinement time
tau_e=Wp./Ploss;
figure
plot(tt*1000, tau_e.*1e6, 'r-') ;
xlabel('Time [ms]')
ylabel('\tau_e [\mus]')
xlim([15 20])
ylim([0 50])
% set(gca,'yscale','log')