%% symmetry oscillation 


% Peak to peak analysis 
% Select all main peaks (likely eruptions) and characterise: 
%  - temperature increase rate 
%  - temperature decrease rate
%  - time interval between peaks
% 
% Plotting: 
% Figures with peak selection  
% Figure with bars plots for T increase and decrease 

clear, close all
F = 1; 

%% Night 20 Strokkur 
night20 = xlsread('Strokkur_20_06_all.xlsx','Sheet1');


%% 
Ts20   = night20(16100:25900, end); 
Td20   = night20(16100:25900, end-1);
X20    = 0:length(Ts20)-1;
dTs20  = detrend(Ts20); 
idTs20 = -dTs20;
dTd20  = detrend(Td20); 
idTd20 = -dTd20;
% Strokkur shallow ~9m 
[~, loc_Ts20]   = findpeaks(dTs20, 'MinPeakDistance', 90, 'MinPeakProminence',1);   % get maxima
[~, iloc_Ts20]  = findpeaks(idTs20, 'MinPeakDistance', 60, 'MinPeakProminence',1.25); % get minima
% quick and dirty fix 
id = find(X20(iloc_Ts20)==5307); 
iloc_Ts20 = setdiff(iloc_Ts20, iloc_Ts20(id));
iloc_Ts20 = (iloc_Ts20(1:end-1));


% strokkur deep 
[~, loc_Td20]   = findpeaks(dTd20, 'MinPeakDistance', 60, 'MinPeakProminence',1, 'MinPeakHeight',-2);   % get maxima
[~, iloc_Td20]  = findpeaks(idTd20, 'MinPeakDistance', 60, 'MinPeakProminence',1); % get minima
iloc_Td20 = iloc_Td20(1:end-1);


% Matrix gathering info on eruption cycle
DATA20s      = zeros(size(iloc_Ts20,1),8);
DATA20s(:,2) = diff(loc_Ts20);
DATA20s(:,1) = 1; 
DATA20s(:,3) = abs(dTs20(iloc_Ts20) - dTs20(loc_Ts20(1:end-1)));
DATA20s(:,4) = X20(iloc_Ts20) - X20(loc_Ts20(1:end-1));
DATA20s(:,5) = DATA20s(:,3)./DATA20s(:,4).*60;
DATA20s(:,6) = dTs20(loc_Ts20(2:end)) - dTs20(iloc_Ts20);
DATA20s(:,7) = X20(loc_Ts20(2:end)) - X20(iloc_Ts20);
DATA20s(:,8) = DATA20s(:,6)./DATA20s(:,7).*60;

DATA20d      = zeros(size(iloc_Td20,1),8);
DATA20d(:,2) = diff(loc_Td20);
DATA20d(:,1) = 1; 
DATA20d(:,3) = abs(dTd20(iloc_Td20) - dTd20(loc_Td20(1:end-1)));
DATA20d(:,4) = X20(iloc_Td20) - X20(loc_Td20(1:end-1));
DATA20d(:,5) = DATA20d(:,3)./DATA20d(:,4).*60;
DATA20d(:,6) = dTd20(loc_Td20(2:end)) - dTd20(iloc_Td20);
DATA20d(:,7) = X20(loc_Td20(2:end)) - X20(iloc_Td20);
DATA20d(:,8) = DATA20d(:,6)./DATA20d(:,7).*60;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Night 21 Strokkur and Geysir
night21   = load('night21.txt');

Tg21   = night21(:,end-1);
Ts21   = night21(:,end);
X21    = 0:length(Ts21)-1;
dTg21  = detrend(Tg21);
dTs21  = detrend(Ts21);
idTs21 = -dTs21; 
idTg21 = -dTg21; 

% Strokkur
[~, loc_Ts21]   = findpeaks(dTs21, 'MinPeakDistance', 60, 'MinPeakProminence',1);   % get maxima
[~, iloc_Ts21]  = findpeaks(idTs21, 'MinPeakDistance', 60, 'MinPeakProminence',1);  % get minima
iloc_Ts21 = iloc_Ts21(1:end-1);

DATA21S      = zeros(size(iloc_Ts21,1),8);
DATA21S(:,2) = diff(loc_Ts21);
DATA21S(:,1) = 1; 
DATA21S(:,3) = abs(dTs21(iloc_Ts21) - dTs21(loc_Ts21(1:end-1)));
DATA21S(:,4) = X21(iloc_Ts21) - X21(loc_Ts21(1:end-1));
DATA21S(:,5) = DATA21S(:,3)./DATA21S(:,4).*60;
DATA21S(:,6) = dTs21(loc_Ts21(2:end)) - dTs21(iloc_Ts21);
DATA21S(:,7) = X21(loc_Ts21(2:end)) - X21(iloc_Ts21);
DATA21S(:,8) = DATA21S(:,6)./DATA21S(:,7).*60;

% Geysir
[~, loc_Tg21]   = findpeaks(dTg21,'minPeakDistance',60, 'minPeakProminence', 0.5);  % get maxima
[~, iloc_Tg21]  = findpeaks(idTg21,'minPeakDistance',60, 'minPeakProminence', 0.5); % get minima 

DATA21G      = zeros(size(iloc_Tg21,1),8);
DATA21G(:,2) = diff(loc_Tg21);
DATA21G(:,1) = 1; 
DATA21G(:,3) = abs(dTg21(iloc_Tg21) - dTg21(loc_Tg21(1:end-1)));
DATA21G(:,4) = X21(iloc_Tg21) - X21(loc_Tg21(1:end-1));
DATA21G(:,5) = DATA21G(:,3)./DATA21G(:,4).*60;
DATA21G(:,6) = dTg21(loc_Tg21(2:end)) - dTg21(iloc_Tg21);
DATA21G(:,7) = X21(loc_Tg21(2:end)) - X21(iloc_Tg21);
DATA21G(:,8) = DATA21G(:,6)./DATA21G(:,7).*60;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% load data for night 22 
night22   = load('night22.txt');

Tg22   = night22(:,end-1);
Ts22   = night22(:,end);
X22    = 0:length(Ts22)-1;
dTg22  = detrend(Tg22);
dTs22  = detrend(Ts22);
idTs22 = -dTs22; 
idTg22 = -dTg22; 

[~, loc_Ts22]   = findpeaks(dTs22, 'MinPeakDistance', 60, 'MinPeakProminence',1);   % get maxima
[~, iloc_Ts22]  = findpeaks(idTs22, 'MinPeakDistance', 60, 'MinPeakProminence',1);  % get minima

DATA22S      = zeros(size(iloc_Ts22,1),8);
DATA22S(:,2) = diff(loc_Ts22);
DATA22S(:,1) = 1; 
DATA22S(:,3) = abs(dTs22(iloc_Ts22) - dTs22(loc_Ts22(1:end-1)));
DATA22S(:,4) = X22(iloc_Ts22) - X22(loc_Ts22(1:end-1));
DATA22S(:,5) = DATA22S(:,3)./DATA22S(:,4).*60;
DATA22S(:,6) = dTs22(loc_Ts22(2:end)) - dTs22(iloc_Ts22);
DATA22S(:,7) = X22(loc_Ts22(2:end)) - X22(iloc_Ts22);
DATA22S(:,8) = DATA22S(:,6)./DATA22S(:,7).*60;

[~, loc_Tg22]   = findpeaks(dTg22,'minPeakDistance',60, 'minPeakProminence', 0.5);  
[~, iloc_Tg22]  = findpeaks(idTg22,'minPeakDistance',60, 'minPeakProminence', 0.5);  

DATA22G      = zeros(size(iloc_Tg22,1),8);
DATA22G(:,2) = diff(loc_Tg22);
DATA22G(:,1) = 1; 
DATA22G(:,3) = abs(dTg22(iloc_Tg22) - dTg22(loc_Tg22(1:end-1)));
DATA22G(:,4) = X22(iloc_Tg22) - X22(loc_Tg22(1:end-1));
DATA22G(:,5) = DATA22G(:,3)./DATA22G(:,4).*60;
DATA22G(:,6) = dTg22(loc_Tg22(2:end)) - dTg22(iloc_Tg22);
DATA22G(:,7) = X22(loc_Tg22(2:end)) - X22(iloc_Tg22);
DATA22G(:,8) = DATA22G(:,6)./DATA22G(:,7).*60;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% load data for night 23 
night23   = load('night23_1.txt');

Tg23 = night23(:,end-1);
Ts23 = night23(:,end);
X23 = 0:length(Ts23)-1;
dTg23  = detrend(Tg23);
dTs23  = detrend(Ts23);
idTs23 = -dTs23; 
idTg23 = -dTg23; 

[~, loc_Ts23]   = findpeaks(dTs23, 'MinPeakDistance', 60, 'MinPeakProminence',1);   % get maxima
[~, iloc_Ts23]   = findpeaks(idTs23, 'MinPeakDistance', 60, 'MinPeakProminence',1); % get minima
iloc_Ts23 = iloc_Ts23(1:end-1);

DATA23S      = zeros(size(iloc_Ts23,1),8);
DATA23S(:,2) = diff(loc_Ts23);
DATA23S(:,1) = 1; 
DATA23S(:,3) = abs(dTs23(iloc_Ts23) - dTs23(loc_Ts23(1:end-1)));
DATA23S(:,4) = X23(iloc_Ts23) - X23(loc_Ts23(1:end-1));
DATA23S(:,5) = DATA23S(:,3)./DATA23S(:,4).*60;
DATA23S(:,6) = dTs23(loc_Ts23(2:end)) - dTs23(iloc_Ts23);
DATA23S(:,7) = X23(loc_Ts23(2:end)) - X23(iloc_Ts23);
DATA23S(:,8) = DATA23S(:,6)./DATA23S(:,7).*60;

[~, loc_Tg23]    = findpeaks(dTg23,'minPeakDistance',60, 'minPeakProminence', 0.5);  
[~, iloc_Tg23]   = findpeaks(idTg23,'minPeakDistance',60, 'minPeakProminence', 0.5);  

iloc_Tg23 = iloc_Tg23(2:end-1);

DATA23G      = zeros(size(iloc_Tg23,1),8);
DATA23G(:,2) = diff(loc_Tg23);
DATA23G(:,1) = 1; 
DATA23G(:,3) = abs(dTg23(iloc_Tg23) - dTg23(loc_Tg23(1:end-1)));
DATA23G(:,4) = X23(iloc_Tg23) - X23(loc_Tg23(1:end-1));
DATA23G(:,5) = DATA23G(:,3)./DATA23G(:,4).*60;
DATA23G(:,6) = dTg23(loc_Tg23(2:end)) - dTg23(iloc_Tg23);
DATA23G(:,7) = X23(loc_Tg23(2:end)) - X23(iloc_Tg23);
DATA23G(:,8) = DATA23G(:,6)./DATA23G(:,7).*60;

%% Statistics data 
% Strokkur 
DATA20s_stats = [mean(DATA20s(:,2:end))', std(DATA20s(:,2:end))', min(DATA20s(:,2:end))'...
    max(DATA20s(:,2:end))', median(DATA20s(:,2:end))'];

DATA20d_stats = [mean(DATA20d(:,2:end))', std(DATA20d(:,2:end))', min(DATA20d(:,2:end))'...
    max(DATA20d(:,2:end))', median(DATA20d(:,2:end))'];

DATA21S_stats = [mean(DATA21S(:,2:end))', std(DATA21S(:,2:end))', min(DATA21S(:,2:end))'...
    max(DATA21S(:,2:end))', median(DATA21S(:,2:end))'];

DATA22S_stats = [mean(DATA22S(:,2:end))', std(DATA22S(:,2:end))', min(DATA22S(:,2:end))'...
    max(DATA22S(:,2:end))', median(DATA22S(:,2:end))'];

DATA23S_stats = [mean(DATA23S(:,2:end))', std(DATA23S(:,2:end))', min(DATA23S(:,2:end))'...
    max(DATA23S(:,2:end))', median(DATA23S(:,2:end))'];

% Great Geysir
DATA21G_stats = [mean(DATA21G(:,2:end))', std(DATA21G(:,2:end))', min(DATA21G(:,2:end))'...
    max(DATA21G(:,2:end))', median(DATA21G(:,2:end))'];

DATA22G_stats = [mean(DATA22G(:,2:end))', std(DATA22G(:,2:end))', min(DATA22G(:,2:end))'...
    max(DATA22G(:,2:end))', median(DATA22G(:,2:end))'];

DATA23G_stats = [mean(DATA23G(:,2:end))', std(DATA23G(:,2:end))', min(DATA23G(:,2:end))'...
    max(DATA23G(:,2:end))', median(DATA23G(:,2:end))'];

%% Reorganise data for plotting
% mean, std

% temperature 
COOLING_temp_S = zeros(5, 2); 
COOLING_temp_S(1,:) = [mean(DATA21S(:,3)), std(DATA21S(:,3))];
COOLING_temp_S(2,:) = [mean(DATA20s(:,3)), std(DATA20s(:,3))];
COOLING_temp_S(3,:) = [mean(DATA22S(:,3)), std(DATA22S(:,3))];
COOLING_temp_S(4,:) = [mean(DATA23S(:,3)), std(DATA23S(:,3))];
COOLING_temp_S(5,:) = [mean(DATA20d(:,3)), std(DATA20d(:,3))];

COOLING_temp_G = zeros(3, 2); 
COOLING_temp_G(1,:) = [mean(DATA22G(:,3)), std(DATA22G(:,3))];
COOLING_temp_G(2,:) = [mean(DATA23G(:,3)), std(DATA23G(:,3))];
COOLING_temp_G(3,:) = [mean(DATA21G(:,3)), std(DATA21G(:,3))];

WARMING_temp_S = zeros(5, 2); 
WARMING_temp_S(1,:) = [mean(DATA21S(:,6)), std(DATA21S(:,6))];
WARMING_temp_S(2,:) = [mean(DATA20s(:,6)), std(DATA20s(:,6))];
WARMING_temp_S(3,:) = [mean(DATA22S(:,6)), std(DATA22S(:,6))];
WARMING_temp_S(4,:) = [mean(DATA23S(:,6)), std(DATA23S(:,6))];
WARMING_temp_S(5,:) = [mean(DATA20d(:,6)), std(DATA20d(:,6))];

WARMING_temp_G = zeros(3, 2);
WARMING_temp_G(1,:) = [mean(DATA22G(:,6)), std(DATA22G(:,6))];
WARMING_temp_G(2,:) = [mean(DATA23G(:,6)), std(DATA23G(:,6))];
WARMING_temp_G(3,:) = [mean(DATA21G(:,6)), std(DATA21G(:,6))];

% rate 
COOLING_rate_S = zeros(5, 2); 
COOLING_rate_S(1,:) = [mean(DATA21S(:,5)), std(DATA21S(:,5))];
COOLING_rate_S(2,:) = [mean(DATA20s(:,5)), std(DATA20s(:,5))];
COOLING_rate_S(3,:) = [mean(DATA22S(:,5)), std(DATA22S(:,5))];
COOLING_rate_S(4,:) = [mean(DATA23S(:,5)), std(DATA23S(:,5))];
COOLING_rate_S(5,:) = [mean(DATA20d(:,5)), std(DATA20d(:,5))];

COOLING_rate_G = zeros(3, 2); 
COOLING_rate_G(1,:) = [mean(DATA22G(:,5)), std(DATA22G(:,5))];
COOLING_rate_G(2,:) = [mean(DATA23G(:,5)), std(DATA23G(:,5))];
COOLING_rate_G(3,:) = [mean(DATA21G(:,5)), std(DATA21G(:,5))];

WARMING_rate_S = zeros(5, 2); 
WARMING_rate_S(1,:) = [mean(DATA21S(:,8)), std(DATA21S(:,8))];
WARMING_rate_S(2,:) = [mean(DATA20s(:,8)), std(DATA20s(:,8))];
WARMING_rate_S(3,:) = [mean(DATA22S(:,8)), std(DATA22S(:,8))];
WARMING_rate_S(4,:) = [mean(DATA23S(:,8)), std(DATA23S(:,8))];
WARMING_rate_S(5,:) = [mean(DATA20d(:,8)), std(DATA20d(:,8))];

WARMING_rate_G = zeros(3, 2);
WARMING_rate_G(1,:) = [mean(DATA22G(:,8)), std(DATA22G(:,8))];
WARMING_rate_G(2,:) = [mean(DATA23G(:,8)), std(DATA23G(:,8))];
WARMING_rate_G(3,:) = [mean(DATA21G(:,8)), std(DATA21G(:,8))];


%% Picking Geysir figure
figure()
subplot(3,1,1)
plot(X21/60, dTg21)
hold on 
plot(X21(loc_Tg21)/60, dTg21(loc_Tg21),'gv')
plot(X21(iloc_Tg21)/60, dTg21(iloc_Tg21),'k^')
xlim([0 192])
ylim([-2 2])
title('Geysir N21')
xlabel('Time(min)')
ylabel('Temperature variation (ºC)')
subplot(3,1,2)
plot(X22/60, dTg22)
hold on 
plot(X22(loc_Tg22)/60, dTg22(loc_Tg22),'gv')
plot(X22(iloc_Tg22)/60, dTg22(iloc_Tg22),'k^')
xlim([0 170])
ylim([-2 2])
title('Geysir N22')
xlabel('Time(min)')
ylabel('Temperature variation (ºC)')
subplot(3,1,3)
plot(X23/60, dTg23)
hold on 
plot(X23(loc_Tg23)/60, dTg23(loc_Tg23),'gv')
plot(X23(iloc_Tg23)/60, dTg23(iloc_Tg23),'k^')
xlim([0 253])
ylim([-2 2])
title('Geysir N23')
xlabel('Time(min)')
ylabel('Temperature variation (ºC)')


%% Picking Strokkur 

figure()
subplot(5,1,1)
plot(X21/60, dTs21)
xlim([0 192])
ylim([-12 8])
hold on
plot(X21(loc_Ts21)/60, dTs21(loc_Ts21), 'gv')
plot(X21(iloc_Ts21)/60, dTs21(iloc_Ts21), 'k^')
xlabel('Time(min)')
ylabel('Temperature variation (ºC)')
title('S21')

subplot(5,1,2)
plot(X20/60, dTs20)
xlim([0 163])
ylim([-12 8])
hold on
plot(X20(loc_Ts20)/60, dTs20(loc_Ts20), 'gv')
plot(X20(iloc_Ts20)/60, dTs20(iloc_Ts20), 'k^')
xlabel('Time(min)')
ylabel('Temperature variation (ºC)')
title('S20a')

subplot(5,1,3)
plot(X22/60, dTs22)
xlim([0 170])
ylim([-12 8])
hold on
plot(X22(loc_Ts22)/60, dTs22(loc_Ts22), 'gv')
plot(X22(iloc_Ts22)/60, dTs22(iloc_Ts22), 'k^')
xlabel('Time(min)')
ylabel('Temperature variation (ºC)')
title('S22')

subplot(5,1,4)
plot(X23/60, dTs23)
xlim([0 253])
ylim([-12 8])
hold on
plot(X23(loc_Ts23)/60, dTs23(loc_Ts23), 'gv')
plot(X23(iloc_Ts23)/60, dTs23(iloc_Ts23), 'k^')
xlabel('Time(min)')
ylabel('Temperature variation (ºC)')
title('S23')

subplot(5,1,5)
plot(X20/60, dTd20)
xlim([0 163])
ylim([-12 8])
hold on
plot(X20(loc_Td20)/60, dTd20(loc_Td20), 'gv')
plot(X20(iloc_Td20)/60, dTd20(iloc_Td20), 'k^')
xlabel('Time(min)')
ylabel('Temperature variation (ºC)')
title('S20b')


%%

figure()
subplot(1,2,1)
errorbar(COOLING_temp_G(:,1), COOLING_temp_G(:,2),'bv')
hold on 
errorbar(WARMING_temp_G(:,1), WARMING_temp_G(:,2),'r^')
subplot(1,2,2)
errorbar(COOLING_temp_S(:,1), COOLING_temp_S(:,2),'bv')
hold on 
errorbar(WARMING_temp_S(:,1), WARMING_temp_S(:,2),'r^')

figure()
subplot(1,2,1)
errorbar(COOLING_rate_G(:,1), COOLING_rate_G(:,2),'bv')
hold on 
errorbar(WARMING_rate_G(:,1), WARMING_rate_G(:,2),'r^')
ylim([0 1])
subplot(1,2,2)
errorbar(COOLING_rate_S(:,1), COOLING_rate_S(:,2),'bv')
hold on 
errorbar(WARMING_rate_S(:,1), WARMING_rate_S(:,2),'r^')
ylim([0 10])


%% 
COOLING_temp_all = [COOLING_temp_G; COOLING_temp_S];
WARMING_temp_all = [WARMING_temp_G; WARMING_temp_S];
COOLING_rate_all = [COOLING_rate_G; COOLING_rate_S];
WARMING_rate_all = [WARMING_rate_G; WARMING_rate_S];


figure()
subplot(1,2,1)
errorbar(COOLING_temp_all(:,1), COOLING_temp_all(:,2),'bv')
hold on 
errorbar(WARMING_temp_all(:,1), WARMING_temp_all(:,2),'r^')
xlim([0 10])
ylim([0 11])
subplot(1,2,2)
errorbar(COOLING_rate_all(:,1), COOLING_rate_all(:,2),'bv')
hold on 
errorbar(WARMING_rate_all(:,1), WARMING_rate_all(:,2),'r^')
xlim([0 10])
ylim([0 9])


