% %% Script working on data from 20th 
% Compare characteristics of Eruption between depth 
% 
%
% 1st col: type 1: single - 2: double - 3: triple
% 2nd col: TAE (sec)
% 3rd col: dT cooling (ºC)
% 4th col: dt cooling (sec)
% 5th col: dT/dt cooling (ºC/sec)
% 6th col: dT warming (ºC)
% 7th col: dt warming (sec)
% 8th col: dT/dt warming (ºC/sec)


clear, close all
Fs = 1; 

%% selection criteron adapted for the night 20 - does not work as well as for the other night 
% much more noise
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  
[~, 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);


[~, im_Ts20]  = findpeaks(idTs20, 'MinPeakDistance', 120, 'MinPeakHeight',3.5); 
% quick and dirty fix 
id_rm = [1100; 7410; 7553];
im_Ts20 = setdiff(im_Ts20, id_rm);

id = find(ismember(iloc_Ts20,im_Ts20));
imloc_Ts20  = iloc_Ts20(id);
mloc_Ts20   = loc_Ts20(id);
isloc_Ts20  = setdiff(iloc_Ts20, imloc_Ts20);
sloc_Ts20   = setdiff(loc_Ts20, mloc_Ts20);

% Matrix gathering info on eruption cycle
DATA20s      = zeros(size(iloc_Ts20,1),8);
DATA20s(:,2) = diff(loc_Ts20);
DATA20s(:,1) = 1; DATA20s(id, 1) = 2;
DATA20s(:,3) = 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);
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);

%% Strokkur deep
% 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);

% find multiple eruption 
[~, im_Td20]  = findpeaks(idTd20, 'MinPeakDistance', 120, 'MinPeakHeight',5); 
id = find(ismember(iloc_Td20,im_Td20));
imloc_Td20  = iloc_Td20(id);
mloc_Td20   = loc_Td20(id);
isloc_Td20  = setdiff(iloc_Td20, imloc_Td20);
sloc_Td20   = setdiff(loc_Td20, mloc_Td20);

% Matrix gathering info on eruption cycle
DATA20d      = zeros(size(iloc_Td20,1),8);
DATA20d(:,2) = diff(loc_Td20);
DATA20d(:,1) = 1; DATA20d(id, 1) = 2;
DATA20d(:,3) = 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);
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);

id_single = find(DATA20s(:,1) == 1);
id_double = find(DATA20s(:,1) ==2);

% id_triple = find(DATA20s(:,2)> 408 & abs(DATA20s(:,5)) > 0.076);
% id_triple = find(DATA20s(:,2)> 408 & abs(DATA20s(:,3)) > 11);
% id_triple = find( abs(DATA20s(:,3)) > 12 & abs(DATA20s(:,5)) > 0.076);


%% 
DATA_all_s = [mean(DATA20s(:,2:end))', std(DATA20s(:,2:end))', ...
    min(DATA20s(:,2:end))', max(DATA20s(:,2:end))'];

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

DATA_single_s = [mean(DATA20s(id_single,2:end))', std(DATA20s(id_single,2:end))', ...
    min(DATA20s(id_single,2:end))', max(DATA20s(id_single,2:end))', median(DATA20s(id_single,2:end))'];

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

DATA_double_s = [mean(DATA20s(id_double,2:end))', std(DATA20s(id_double,2:end))', ...
    min(DATA20s(id_double,2:end))', max(DATA20s(id_double,2:end))', median(DATA20s(id_double,2:end))'];

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

 



%%

TAE_all = zeros(2,4);
TAE_all(1,:) = [mean(DATA20s(:,2)), std(DATA20s(:,2)) , ...
            min(DATA20s(:,2)),  max(DATA20s(:,2))]; 
        
TAE_all(2,:) = [mean(DATA20d(:,2)), std(DATA20d(:,2)) , ...
            min(DATA20d(:,2)),  max(DATA20d(:,2))]; 

TAE_single = zeros(2,4);        
TAE_single(1,:) = [mean(DATA20s(id_single,2)), std(DATA20s(id_single,2)) , ...
            min(DATA20s(id_single,2)),  max(DATA20s(id_single,2))]; 
        
TAE_single(2,:) = [mean(DATA20d(id_single,2)), std(DATA20d(id_single,2)) , ...
            min(DATA20d(id_single,2)),  max(DATA20d(id_single,2))]; 
        
% temperature variation of cooling phase 
COOLING_temp_all = zeros(2,4);
COOLING_temp_all(1,:) = [mean(DATA20s(:,3)), std(DATA20s(:,3)) , ...
                    min(DATA20s(:,3)),  max(DATA20s(:,3))];

COOLING_temp_all(2,:) = [mean(DATA20d(:,3)), std(DATA20d(:,3)) ...
                    min(DATA20d(:,3)),  max(DATA20d(:,3))]; 

COOLING_temp_single = zeros(2,4);
COOLING_temp_single(1,:) = [mean(DATA20s(id_single,3)), std(DATA20s(id_single,3)) , ...
                    min(DATA20s(id_single,3)),  max(DATA20s(id_single,3))];

COOLING_temp_single(2,:) = [mean(DATA20d(id_single,3)), std(DATA20d(id_single,3)) ...
                    min(DATA20d(id_single,3)),  max(DATA20d(id_single,3))]; 

COOLING_temp_double = zeros(2,4);
COOLING_temp_double(1,:) = [mean(DATA20s(id_double,3)), std(DATA20s(id_double,3)) , ...
                    min(DATA20s(id_double,3)),  max(DATA20s(id_double,3))];

COOLING_temp_double(2,:) = [mean(DATA20d(id_double,3)), std(DATA20d(id_double,3)) ...
                    min(DATA20d(id_double,3)),  max(DATA20d(id_double,3))]; 

                
% Duration cooling phase                
COOLING_time_all = zeros(2,4);
COOLING_time_all(1,:) = [mean(DATA20s(:,4)), std(DATA20s(:,4)) , ...
                    min(DATA20s(:,4)),  max(DATA20s(:,4))];

COOLING_time_all(2,:) = [mean(DATA20d(:,4)), std(DATA20d(:,4)) ...
                    min(DATA20d(:,4)),  max(DATA20d(:,4))]; 

COOLING_time_single = zeros(2,4);
COOLING_time_single(1,:) = [mean(DATA20s(id_single,4)), std(DATA20s(id_single,4)) , ...
                    min(DATA20s(id_single,4)),  max(DATA20s(id_single,4))];

COOLING_time_single(2,:) = [mean(DATA20d(id_single,4)), std(DATA20d(id_single,4)) ...
                    min(DATA20d(id_single,4)),  max(DATA20d(id_single,4))]; 
                
COOLING_time_double = zeros(2,4);
COOLING_time_double(1,:) = [mean(DATA20s(id_double,4)), std(DATA20s(id_double,4)) , ...
                    min(DATA20s(id_double,4)),  max(DATA20s(id_double,4))];

COOLING_time_double(2,:) = [mean(DATA20d(id_double,4)), std(DATA20d(id_double,4)) ...
                    min(DATA20d(id_double,4)),  max(DATA20d(id_double,4))]; 


% Cooling rate                
COOLING_rate_all = zeros(2,4);
COOLING_rate_all(1,:) = [mean(DATA20s(:,5)), std(DATA20s(:,5)) , ...
                    min(DATA20s(:,5)),  max(DATA20s(:,5))];

COOLING_rate_all(2,:) = [mean(DATA20d(:,5)), std(DATA20d(:,5)) ...
                    min(DATA20d(:,5)),  max(DATA20d(:,5))]; 

COOLING_rate_single = zeros(2,4);
COOLING_rate_single(1,:) = [mean(DATA20s(id_single,5)), std(DATA20s(id_single,5)) , ...
                    min(DATA20s(id_single,5)),  max(DATA20s(id_single,5))];

COOLING_rate_single(2,:) = [mean(DATA20d(id_single,5)), std(DATA20d(id_single,5)) ...
                    min(DATA20d(id_single,5)),  max(DATA20d(id_single,5))]; 
                
COOLING_rate_double = zeros(2,4);
COOLING_rate_double(1,:) = [mean(DATA20s(id_double,5)), std(DATA20s(id_double,5)) , ...
                    min(DATA20s(id_double,5)),  max(DATA20s(id_double,5))];

COOLING_rate_double(2,:) = [mean(DATA20d(id_double,5)), std(DATA20d(id_double,5)) ...
                    min(DATA20d(id_double,5)),  max(DATA20d(id_double,5))]; 
                

% temperature variation of warming phase 
WARMING_temp_all = zeros(2,4);
WARMING_temp_all(1,:) = [mean(DATA20s(:,6)), std(DATA20s(:,6)) , ...
                    min(DATA20s(:,6)),  max(DATA20s(:,6))];

WARMING_temp_all(2,:) = [mean(DATA20d(:,6)), std(DATA20d(:,6)) ...
                    min(DATA20d(:,6)),  max(DATA20d(:,6))]; 

WARMING_temp_single = zeros(2,4);
WARMING_temp_single(1,:) = [mean(DATA20s(id_single,6)), std(DATA20s(id_single,6)) , ...
                    min(DATA20s(id_single,6)),  max(DATA20s(id_single,6))];

WARMING_temp_single(2,:) = [mean(DATA20d(id_single,6)), std(DATA20d(id_single,6)) ...
                    min(DATA20d(id_single,6)),  max(DATA20d(id_single,6))]; 
                
WARMING_temp_double = zeros(2,4);
WARMING_temp_double(1,:) = [mean(DATA20s(id_double,6)), std(DATA20s(id_double,6)) , ...
                    min(DATA20s(id_double,6)),  max(DATA20s(id_double,6))];

WARMING_temp_double(2,:) = [mean(DATA20d(id_double,6)), std(DATA20d(id_double,6)) ...
                    min(DATA20d(id_double,6)),  max(DATA20d(id_double,6))];                 

% Duration warming phase
WARMING_time_all = zeros(2,4);
WARMING_time_all(1,:) = [mean(DATA20s(:,7)), std(DATA20s(:,7)) , ...
                    min(DATA20s(:,7)),  max(DATA20s(:,7))];

WARMING_time_all(2,:) = [mean(DATA20d(:,7)), std(DATA20d(:,7)) ...
                    min(DATA20d(:,7)),  max(DATA20d(:,7))]; 

WARMING_time_single = zeros(2,4);
WARMING_time_single(1,:) = [mean(DATA20s(id_single,7)), std(DATA20s(id_single,7)) , ...
                    min(DATA20s(id_single,7)),  max(DATA20s(id_single,7))];

WARMING_time_single(2,:) = [mean(DATA20d(id_single,7)), std(DATA20d(id_single,7)) ...
                    min(DATA20d(id_single,7)),  max(DATA20d(id_single,7))]; 

WARMING_time_double = zeros(2,4);
WARMING_time_double(1,:) = [mean(DATA20s(id_double,7)), std(DATA20s(id_double,7)) , ...
                    min(DATA20s(id_double,7)),  max(DATA20s(id_double,7))];

WARMING_time_double(2,:) = [mean(DATA20d(id_double,7)), std(DATA20d(id_double,7)) ...
                    min(DATA20d(id_double,7)),  max(DATA20d(id_double,7))];                 

% Rate warming phase                 
WARMING_rate_all = zeros(2,4);
WARMING_rate_all(1,:) = [mean(DATA20s(:,8)), std(DATA20s(:,8)) , ...
                    min(DATA20s(:,8)),  max(DATA20s(:,8))];

WARMING_rate_all(2,:) = [mean(DATA20d(:,8)), std(DATA20d(:,8)) ...
                    min(DATA20d(:,8)),  max(DATA20d(:,8))]; 

WARMING_rate_single = zeros(2,4);
WARMING_rate_single(1,:) = [mean(DATA20s(id_single,8)), std(DATA20s(id_single,8)) , ...
                    min(DATA20s(id_single,8)),  max(DATA20s(id_single,8))];

WARMING_rate_single(2,:) = [mean(DATA20d(id_single,8)), std(DATA20d(id_single,8)) ...
                    min(DATA20d(id_single,8)),  max(DATA20d(id_single,8))];         

WARMING_rate_double = zeros(2,4);
WARMING_rate_double(1,:) = [mean(DATA20s(id_double,8)), std(DATA20s(id_double,8)) , ...
                    min(DATA20s(id_double,8)),  max(DATA20s(id_double,8))];

WARMING_rate_double(2,:) = [mean(DATA20d(id_double,8)), std(DATA20d(id_double,8)) ...
                    min(DATA20d(id_double,8)),  max(DATA20d(id_double,8))];                 
               
%% Eruption characteristics - all - vs detph
figure()
% dT
subplot(1,3,1)
errorbar(abs(COOLING_temp_all(:,1)), abs(COOLING_temp_all(:,2)),'bv')
hold on 
errorbar(WARMING_temp_all(:,1), WARMING_temp_all(:,2),'r^')
xlabel('depth')
ylabel('\Delta T (ºC)')
ylim([0 16])
% dt
subplot(1,3,2)
errorbar(COOLING_time_all(:,1)./60, COOLING_time_all(:,2)./60,'bv')
hold on 
errorbar(WARMING_time_all(:,1)./60, WARMING_time_all(:,2)./60,'r^')
xlabel('depth')
ylabel('\Delta t (min)')
ylim([0 8])
% rate
subplot(1,3,3)
errorbar(abs(COOLING_rate_all(:,1)).*60, abs(COOLING_rate_all(:,2)).*60,'bv')
hold on 
errorbar(WARMING_rate_all(:,1).*60, WARMING_rate_all(:,2).*60,'r^')
xlabel('depth')
ylabel(' \Delta T/ \Delta t (ºC)')
ylim([0 10])
%% Eruption characteristics - split single vs. double eruption - vs detph
figure()
% dT 
subplot(2,3,1)
errorbar(abs(COOLING_temp_single(:,1)), abs(COOLING_temp_single(:,2)),'bv')
hold on 
errorbar(WARMING_temp_single(:,1), WARMING_temp_single(:,2),'r^')
xlabel('depth')
ylabel('\Delta T (ºC)')
ylim([0 16])
subplot(2,3,4)
errorbar(abs(COOLING_temp_double(:,1)), abs(COOLING_temp_double(:,2)),'bv')
hold on
errorbar(WARMING_temp_double(:,1), WARMING_temp_double(:,2),'r^')
xlabel('depth')
ylabel('\Delta T (ºC)')
ylim([0 16])
% dt
subplot(2,3,2)
errorbar(COOLING_time_single(:,1)./60, COOLING_time_single(:,2)./60,'bv')
hold on 
errorbar(WARMING_time_single(:,1)./60, WARMING_time_single(:,2)./60,'r^')
xlabel('depth')
ylabel('\Delta t (min)')
ylim([0 8])
subplot(2,3,5)
errorbar(COOLING_time_double(:,1)./60, COOLING_time_double(:,2)./60,'bv')
hold on
errorbar(WARMING_time_double(:,1)./60, WARMING_time_double(:,2)./60,'r^')
xlabel('depth')
ylabel('\Delta t (min)')
ylim([0 8])
% rate
subplot(2,3,3)
errorbar(abs(COOLING_rate_single(:,1)).*60, abs(COOLING_rate_single(:,2)).*60,'bv')
hold on 
errorbar(WARMING_rate_single(:,1).*60, WARMING_rate_single(:,2).*60,'r^')
xlabel('depth')
ylabel(' \Delta T/\Delta t (ºC)')
ylim([0 10])
subplot(2,3,6)
errorbar(abs(COOLING_rate_double(:,1)).*60, abs(COOLING_rate_double(:,2)).*60,'bv')
hold on
errorbar(WARMING_rate_double(:,1).*60, WARMING_rate_double(:,2).*60,'r^')
xlabel('depth')
ylabel(' \Delta T/ \Deltat (ºC)')
ylim([0 10])


%% Figure all, single, double -- showing depth -- double where analyse based on visual inspection and 
% obtention of at least 2 criteria TAE / dT cool and dT/dt rate 
figure()
% dT
subplot(3,3,1)
errorbar(abs(COOLING_temp_all(:,1)), abs(COOLING_temp_all(:,2)),'bv')
hold on 
errorbar(WARMING_temp_all(:,1), WARMING_temp_all(:,2),'r^')
xlabel('depth')
ylabel('\Delta T (ºC)')
ylim([0 16])
subplot(3,3,4)
errorbar(abs(COOLING_temp_single(:,1)), abs(COOLING_temp_single(:,2)),'bv')
hold on 
errorbar(WARMING_temp_single(:,1), WARMING_temp_single(:,2),'r^')
xlabel('depth')
ylabel('\Delta T (ºC)')
ylim([0 16])
subplot(3,3,7)
errorbar(abs(COOLING_temp_double(:,1)), abs(COOLING_temp_double(:,2)),'bv')
hold on
errorbar(WARMING_temp_double(:,1), WARMING_temp_double(:,2),'r^')
xlabel('depth')
ylabel('\Delta T (ºC)')
ylim([0 16])
% dt
subplot(3,3,2)
errorbar(COOLING_time_all(:,1)./60, COOLING_time_all(:,2)./60,'bv')
hold on 
errorbar(WARMING_time_all(:,1)./60, WARMING_time_all(:,2)./60,'r^')
xlabel('depth')
ylabel('\Delta t (min)')
ylim([0 8])
subplot(3,3,5)
errorbar(COOLING_time_single(:,1)./60, COOLING_time_single(:,2)./60,'bv')
hold on 
errorbar(WARMING_time_single(:,1)./60, WARMING_time_single(:,2)./60,'r^')
xlabel('depth')
ylabel('\Delta t (min)')
ylim([0 8])
subplot(3,3,8)
errorbar(COOLING_time_double(:,1)./60, COOLING_time_double(:,2)./60,'bv')
hold on
errorbar(WARMING_time_double(:,1)./60, WARMING_time_double(:,2)./60,'r^')
xlabel('depth')
ylabel('\Delta t (min)')
ylim([0 8])

% rate
subplot(3,3,3)
errorbar(abs(COOLING_rate_all(:,1)).*60, abs(COOLING_rate_all(:,2)).*60,'bv')
hold on 
errorbar(WARMING_rate_all(:,1).*60, WARMING_rate_all(:,2).*60,'r^')
xlabel('depth')
ylabel(' \Delta T/ \Delta t (ºC)')
ylim([0 10])
subplot(3,3,6)
errorbar(abs(COOLING_rate_single(:,1)).*60, abs(COOLING_rate_single(:,2)).*60,'bv')
hold on 
errorbar(WARMING_rate_single(:,1).*60, WARMING_rate_single(:,2).*60,'r^')
xlabel('depth')
ylabel(' \Delta T/\Delta t (ºC)')
ylim([0 10])
subplot(3,3,9)
errorbar(abs(COOLING_rate_double(:,1)).*60, abs(COOLING_rate_double(:,2)).*60,'bv')
hold on
errorbar(WARMING_rate_double(:,1).*60, WARMING_rate_double(:,2).*60,'r^')
xlabel('depth')
ylabel(' \Delta T/ \Deltat (ºC)')
ylim([0 10])





%% Select eruptions
figure()
plot(X20/60, Ts20)
xlim([0 163])
hold on
plot(X20(sloc_Ts20)/60, Ts20(sloc_Ts20), 'gv')
plot(X20(isloc_Ts20)/60, Ts20(isloc_Ts20), 'k^')
plot(X20(mloc_Ts20)/60, Ts20(mloc_Ts20), 'rv')
plot(X20(imloc_Ts20)/60, Ts20(imloc_Ts20), 'b^')
plot(X20/60, Td20)
plot(X20(sloc_Td20)/60, Td20(sloc_Td20), 'gv')
plot(X20(isloc_Td20)/60, Td20(isloc_Td20), 'k^')
plot(X20(mloc_Td20)/60, Td20(mloc_Td20), 'rv')
plot(X20(imloc_Td20)/60, Td20(imloc_Td20), 'b^')
xlabel('Time(min)')
ylabel('Temperature variation (ºC)')
title('Strokkur Night 20 - deep')
