%% Shot numbers
shotNum_vac = 39271;
shotNum = 39272;


%% Path for save

path_gif = 'D:\Czech PhD\GOLEM\matlab\magnetic_coils\';

path_png = 'D:\Czech PhD\GOLEM\matlab\magnetic_coils\';

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

% shot #11364
% data.coils = coils; %plasma data #11364
% data.t_coils = linspace(0, 39999*1e-6,40000)';
% data_vac.coils = coils_vac; %vacuum data #11365
% data_vac.t_coils = linspace(0, 39999*1e-6,40000)';
% data.t_plasma_start = 12;
% data.t_plasma_end = 20;
%% Coils parameters

% 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;
% for i = 1:16
%   subplot(4,4,i);
%   hold on;
%   plot(data.t_coils*1000, data.Bcoils(:,i));
%   %plot(data.t_basic, data.Bt.*Aphi(i));
%   plot(data.t_coils*1000, data_vac.Bcoils(:,i));
%   plot(data.t_coils*1000, data.Bnorm(:,i), 'k', 'LineWidth', 1);
%   hold off;
%   title(i);
%   grid on;
% end


%% crosscorrelation
step = 128;
t1 = find(data.t_coils*1000 > data.t_plasma_start, 1, 'first');
t2 = find(data.t_coils*1000 < data.t_plasma_end, 1, 'last');

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

%% polarplot
lag_num = 160;

figure;
for k = 1:length(t)
    clf;
    y = zeros(1,size(data.Bnorm,2)+1);
    radius = 0;
    hold on;
    
    % normalize, averege and prepare data from different coils for plotting
    for i = 1:size(data.Bnorm,2)
      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)];
      s = (max(x) - min(x));
      x = x./s;
      radius = abs(min(x))*1 + radius;
      x = x + abs(min(x))*1;
      
      y = y + x./mean(x);
    end   
    
    % spline approximation
    [x1,fx, phasespline]=phaseapprox((0:16)/8*pi, y/size(data.Bnorm,2), 0.9999);
    
    polar((0:size(data.Bnorm,2))/8*pi, y/size(data.Bnorm,2), 'or'); %plot points
    polar(x1, fx, 'b');  %plot approximation
    polar(x1, ones(size(fx)).*mean(y)/16, 'k--');  %plot medium level
    text(1.2, 1.5, [num2str(t(k) + data.t_plasma_start) ' ms']) %plot time
    xlim([-2 2]);
    ylim([-2 2]);
    grid on;
    hold off;
    frame(k) = getframe(gcf); % save frames for gif
   % pause(0.2);
   
   % Save pictures with 1 ms step
    if mod(k, 16) == 1
        saveas(gca, [path_png 'graph_' num2str(fix(k/16) + 1) '.png'])
    end
end

close all;

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