%% Shot numbers

ShotNo = 0;
vacuum_shot = 0;

shotNum_vac = vacuum_shot;
shotNum = ShotNo;

%% Download data

data = data_golem_coils(shotNum);
data_vac = data_golem_coils(shotNum_vac);

% Coil's normalization parameters
Aeff = 1e-4*[-68.9 -140.70 138.9 -140.4 -68.6 134.5 -134.3 142.5 -67.6 -142.8 -140.4 -138 -76.3 -142.2 -139.8 -139.3]';

% Normalization and integration to get B_theta

data.Bcoils = integ(data.coils, 1e-6)./repelem((Aeff)', size(data.t_coils,1), size(data.t_coils,2));
%data.Bcoils(:,12) = data.Bcoils(:,12)*10;
data_vac.Bcoils = integ(data_vac.coils, 1e-6)./repelem((Aeff)', size(data.t_coils,1), size(data.t_coils,2));
% remove external fields
data.Bnorm = data.Bcoils-data_vac.Bcoils;

%% Check offset removal by vacuum shot
figure('visible', 'on');
hold on

for i = 1:16
subplot(4,4,i);
plot(data.t_coils*1000, data.Bcoils(:,i));
plot(data.t_coils*1000, data_vac.Bcoils(:,i));
plot(data.t_coils*1000, data.Bnorm(:,i), 'k');
title(['coil ' num2str(i)]);
xlabel('t, ms');
grid on;
end
hold off

set(gcf, 'Position', [0 0 1000 800]);
saveas(gcf, 'offset.png');

%% safety factor calculation
B_theta =(1.25663706e-6)*data.Ip*1000/(2*pi*0.093);

q = (8.5/40)*data.Bt./B_theta;

% discharge time
t1 = find(data.t_basic > data.t_plasma_start, 1, 'first');
t2 = find(data.t_basic < data.t_plasma_end, 1, 'last');

figure('visible','on');

hold on
plot(data.t_basic(t1+1000:t2),q(t1+1000:t2));
grid on;
xlabel('t, ms');
ylabel('q');
title('Safety factor');
set(gcf, 'Position', [0 0 1000 800]);

hold off
saveas(gcf, 'q.png');

%% Spectrogram
step = 128*2; % window

figure('visible','on');
hold on

%figure;
for i = 1:1
[sp, t, f] = Spectrs(data.coils(t1:t2,i)-data_vac.coils(t1:t2,i), step, 1, step/2, 1);

%subplot(4,4,i);
n = find(f*1000 < 200, 1,'last'); % frequency limit
h = pcolor(t+data.t_coils(t1-1)*1000,f(1:n)*1000, sp(1:n,:));
set(h,'EdgeColor','None');
ylabel('f, kHz');
xlabel('t, ms');

title('Spectrogram');
set(gcf, 'Position', [0 0 1000 800]);

end

hold off
saveas(gcf, 'Spectrogram.png');
%close(fh)

%% crosscorrelation
step = 128*2; %window

[r, lag, t] = corr_cross(diff(data.Bnorm(t1:t2,:)), step, 1, step/2, 1);

%% polarplot gif animation
% plot a croscorrelation coefficients for the base coil as an animation

lag_num = step; % lag = 0

radius = 1;
for i = 1:size(r,2) % set the mean radius of the polar plot
r(lag_num,i,:) = (r(lag_num,i,:) + radius);
end

figure('visible','on');

for k = 1:length(t)
clf;
y = zeros(1,size(data.Bnorm,2)+1);

hold on;

for i = 1:size(data.Bnorm,2)
x2(i,:) = r(lag_num,(i-1)*size(data.Bnorm,2)+1:i*size(data.Bnorm,2),k);
end
% normalize, average and prepare data from different coils for plotting
for i = 1:1 % display only the first coils's correlation
x = [r(lag_num,(i-1)*size(data.Bnorm,2)+1:i*size(data.Bnorm,2),k) r(lag_num,(i-1)*size(data.Bnorm,2)+1,k)];

x1= x;
x1(end) = [];

x1 = x1 - mean(x1);
st(k,i) = std(x1); % activity level

y = x;
x = (0:16)/8*pi;

% spline approximation
[x1,fx, phasespline]=phaseapprox(x, y, 0.9999);

%subplot(4,4,i);
hold on;
polar(x, y, 'or'); %plot points
polar(x1, fx, 'b'); %plot approximation
polar(x1, ones(size(fx)), 'k--'); %plot medium level
text(1.2, 1.5, [num2str(t(k) + data.t_plasma_start) ' ms']) %plot time
[~, t_num] = min(abs((t(k) + data.t_plasma_start)-data.t_basic));
text(-1.2, 1.5, ['q = ' num2str(q(t_num))]); %plot q
xlim([-2 2]);
ylim([-2 2]);
grid on;
set(gcf, 'Position', [0 0 1000 800]);
end

save frame
frame(k) = getframe(gcf);
% pause(0.3);

end

%% Save animation as .gif
for n = 1:size(frame,2)
[A,map] = rgb2ind(frame2im(frame(n)),256);
if n == 1
imwrite(A, map, 'animation.gif','gif','LoopCount',Inf,'DelayTime',0.3);
else
imwrite(A, map, 'animation.gif','gif','WriteMode','append','DelayTime',0.3);
end
end

%% Calculate and plot MHD activity
% MHD activity is estimated as standart deviation of coils signals.
% Points for all coils are ploted. The Black line is a base coil

for k = 1:length(t)
y = zeros(1,size(data.Bnorm,2)+1);

hold on;

for i = 1:size(data.Bnorm,2)
x2(i,:) = r(lag_num,(i-1)*size(data.Bnorm,2)+1:i*size(data.Bnorm,2),k);
end

for i = 1:size(data.Bnorm,2)

x1= r(lag_num,(i-1)*size(data.Bnorm,2)+1:i*size(data.Bnorm,2),k);

x1 = x1 - mean(x1);
st(k,i) = std(x1); % activity level

end
end

figure('visible','on');

hold on;
plot(t + data.t_plasma_start, st, 'ro');
plot(t + data.t_plasma_start, st(:,1), 'ko-', 'LineWidth', 2);
grid on;
title('MHD activity level');
xlabel('t, ms');
set(gcf, 'Position', [0 0 1000 800]);

saveas(gcf, 'MHD_activity.png');

%% polarplot
% plot timeslices of crosscorrelation coefficients for the base coil

time_step = ceil(length(t)/16);
s = 0;

figure('visible','on');

%figure;
for k = 1:time_step:length(t)
y = zeros(1,size(data.Bnorm,2)+1);

hold on;

s = s + 1;
% normalize, average and prepare data from different coils for plotting
for i = 1:1%size(data.Bnorm,2) % display only the first coils's correlation
x = [r(lag_num,(i-1)*size(data.Bnorm,2)+1:i*size(data.Bnorm,2),k) r(lag_num,(i-1)*size(data.Bnorm,2)+1,k)];

y = x;
x = (0:16)/8*pi;

% spline approximation
[x1,fx, phasespline]=phaseapprox(x, y, 0.9999);

subplot(4,4,s);
hold on;
polar(x, y, 'or'); %plot points
polar(x1, fx, 'b'); %plot approximation
polar(x1, ones(size(fx)), 'k--'); %plot medium level
text(1, 1.7, [num2str(t(k) + data.t_plasma_start,3) 'ms']) %plot time
[~, t_num] = min(abs((t(k) + data.t_plasma_start)-data.t_basic));
text(-1.9, 1.7, ['q =', num2str(q(t_num), 2)]); %plot q
xlim([-2 2]);
ylim([-2 2]);
grid on;
set(gcf, 'Position', [0 0 1200 1000]);
end

end

saveas(gcf, 'pictures.png');