%% Script working on "observed data" only
% Check if picking work correctly and gather data into a matrix for each
% night
% 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
F = 1; 


%% Night 22
night22   = load('night22.txt');
TAE22_real = [134; 184; 183; 223; 348; 240; 591; 210; 281; 352];

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

X22_e    = X22(4292:7316); 
dTs22_e  = dTs22(4292:7316);
idTs22_e = - dTs22_e;

[~, loc_Ts22_e]  = findpeaks(dTs22_e, 'MinPeakDistance', 60, 'MinPeakProminence',1); 
[~, iloc_Ts22_e] = findpeaks(idTs22_e, 'MinPeakDistance', 60, 'MinPeakProminence',1); 
iloc_Ts22_e = iloc_Ts22_e(1:end-1);
 
% find "true eruption" multiple eruptions
[~, im_Ts22] = findpeaks(idTs22_e, 'MinPeakDistance', 60, 'MinPeakProminence',8); 

% split between single and multiple eruption
id = find(ismember(iloc_Ts22_e,im_Ts22));
imloc_Ts22  = iloc_Ts22_e(id);
mloc_Ts22   = loc_Ts22_e(id);
isloc_Ts22  = setdiff(iloc_Ts22_e, imloc_Ts22);
sloc_Ts22   = setdiff(loc_Ts22_e, mloc_Ts22);

% Matrix gathering info on eruption cycle
DATA22      = zeros(size(iloc_Ts22_e,1),8);
DATA22(:,2) = diff(loc_Ts22_e);
DATA22(:,1) = 1; DATA22(id, 1) = 2;
DATA22(:,3) = dTs22_e(iloc_Ts22_e) - dTs22_e(loc_Ts22_e(1:end-1));
DATA22(:,4) = X22_e(iloc_Ts22_e) - X22_e(loc_Ts22_e(1:end-1));
DATA22(:,5) = DATA22(:,3)./DATA22(:,4);
DATA22(:,6) = dTs22_e(loc_Ts22_e(2:end)) - dTs22_e(iloc_Ts22_e);
DATA22(:,7) = X22_e(loc_Ts22_e(2:end)) - X22_e(iloc_Ts22_e);
DATA22(:,8) = DATA22(:,6)./DATA22(:,7);

%% Night 23
night23 = load('night23_1.txt');

Ts23 = night23(:,end);
X23 = 0:length(Ts23)-1;
dTs23  = detrend(Ts23);
TAE23_1_real = [291; 240; 227; 183; 182; 247; 396; 201; 248; 233; 118; 330; 203; 521; 215];
TAE23_2_real = [345; 120; 172; 206; 174; 280; 171; 394; 131; 250; 163; 346; 187; 450; 268; 159; ...
 653; 232; 183; 178; 363; 209; 111; 249; 193; 189; 218; 253; 395; 498; 505; 398; 225; 176]; 

% find "true eruption" all eruption
% part one 
dTs23_1  = dTs23(700:4760);
idTs23_1 = -dTs23_1; 
X23_1    = X23(700:4760);
[val_Ts23_1, loc_Ts23_1]   = findpeaks(dTs23_1, 'MinPeakDistance', 60, 'MinPeakProminence',1); 
[ival_Ts23_1, iloc_Ts23_1] = findpeaks(idTs23_1, 'MinPeakDistance', 60, 'MinPeakProminence',1); 

ival_Ts23_1 = ival_Ts23_1(1:end-1);
iloc_Ts23_1 = iloc_Ts23_1(1:end-1);

% find "true eruption" multiple eruptions
[~, im_Ts23_1] = findpeaks(idTs23_1, 'MinPeakDistance', 60, 'MinPeakProminence',10); 
% split between single and multiple eruption
id = find(ismember(iloc_Ts23_1,im_Ts23_1));
imloc_Ts23_1  = iloc_Ts23_1(id);
mloc_Ts23_1   = loc_Ts23_1(id);
isloc_Ts23_1  = setdiff(iloc_Ts23_1, imloc_Ts23_1);
sloc_Ts23_1   = setdiff(loc_Ts23_1, mloc_Ts23_1);


%%  
% Matrix gathering info on eruption cycle - block 1
DATA23_1      = zeros(size(iloc_Ts23_1,1),8);
DATA23_1(:,2) = diff(loc_Ts23_1);
DATA23_1(:,1) = 1; DATA23_1(id, 1) = 3;
DATA23_1(:,3) = dTs23_1(iloc_Ts23_1) - dTs23_1(loc_Ts23_1(1:end-1));
DATA23_1(:,4) = X23_1(iloc_Ts23_1) - X23_1(loc_Ts23_1(1:end-1));
DATA23_1(:,5) = DATA23_1(:,3)./DATA23_1(:,4);
DATA23_1(:,6) = dTs23_1(loc_Ts23_1(2:end)) - dTs23_1(iloc_Ts23_1);
DATA23_1(:,7) = X23_1(loc_Ts23_1(2:end)) - X23_1(iloc_Ts23_1);
DATA23_1(:,8) = DATA23_1(:,6)./DATA23_1(:,7);

%%
% part two
dTs23_2  = dTs23(5790:end);
idTs23_2 = -dTs23_2;
X23_2    = X23(5790:end);
[val_Ts23_2, loc_Ts23_2]   = findpeaks(dTs23_2, 'MinPeakDistance', 60, 'MinPeakProminence',1);
[ival_Ts23_2, iloc_Ts23_2] = findpeaks(idTs23_2, 'MinPeakDistance', 60, 'MinPeakProminence',1);
ival_Ts23_2 = ival_Ts23_2(1:end-1);
iloc_Ts23_2 = iloc_Ts23_2(1:end-1);

% remove aborted eruption
loc_Ts23_2_wae        = loc_Ts23_2; 
loc_Ts23_2_wae(2)     = NaN; 
loc_Ts23_2_wae(end-5) = NaN;
idNaN                 = isnan(loc_Ts23_2_wae);
loc_Ts23_2_wae        = loc_Ts23_2_wae(~idNaN);

val_Ts23_2_wae        = val_Ts23_2; 
val_Ts23_2_wae(2)     = NaN; 
val_Ts23_2_wae(end-5) = NaN;
idNaN                 = isnan(val_Ts23_2_wae);
val_Ts23_2_wae        = val_Ts23_2_wae(~idNaN);

iloc_Ts23_2_wae        = iloc_Ts23_2; 
iloc_Ts23_2_wae(2)     = NaN; 
iloc_Ts23_2_wae(end-5) = NaN;
idNaN                  = isnan(iloc_Ts23_2_wae);
iloc_Ts23_2_wae        = iloc_Ts23_2_wae(~idNaN);

ival_Ts23_2_wae        = ival_Ts23_2; 
ival_Ts23_2_wae(2)     = NaN; 
ival_Ts23_2_wae(end-5) = NaN;
idNaN                  = isnan(ival_Ts23_2_wae);
ival_Ts23_2_wae        = ival_Ts23_2_wae(~idNaN);

% find multiple eruption and split into single and mult
[~, im_Ts23_2] = findpeaks(idTs23_2, 'MinPeakDistance', 60, 'MinPeakProminence',7); 
im_Ts23_2 = [im_Ts23_2(1:2); im_Ts23_2(4:end)]; 
id = find(ismember(iloc_Ts23_2_wae,im_Ts23_2));
imloc_Ts23_2  = iloc_Ts23_2_wae(id);
mloc_Ts23_2   = loc_Ts23_2_wae(id);
isloc_Ts23_2  = setdiff(iloc_Ts23_2_wae, imloc_Ts23_2);
sloc_Ts23_2   = setdiff(loc_Ts23_2_wae, mloc_Ts23_2);

DATA23_2      = zeros(size(iloc_Ts23_2_wae,1),8);
DATA23_2(:,2) = diff(loc_Ts23_2_wae);
DATA23_2(:,1) = 1; DATA23_2(id, 1) = 2; DATA23_2(id(4), 1) = 3;
DATA23_2(:,3) = dTs23_2(iloc_Ts23_2_wae) - dTs23_2(loc_Ts23_2_wae(1:end-1));
DATA23_2(:,4) = X23_2(iloc_Ts23_2_wae) - X23_2(loc_Ts23_2_wae(1:end-1));
DATA23_2(:,5) = DATA23_2(:,3)./DATA23_2(:,4);
DATA23_2(:,6) = dTs23_2(loc_Ts23_2_wae(2:end)) - dTs23_2(iloc_Ts23_2_wae);
DATA23_2(:,7) = X23_2(loc_Ts23_2_wae(2:end)) - X23_2(iloc_Ts23_2_wae);
DATA23_2(:,8) = DATA23_2(:,6)./DATA23_2(:,7);

%% compare eruption together 
DATA = [DATA22; DATA23_1; DATA23_2];
id_single = find(DATA(:,1)==1); 
id_double = find(DATA(:,1)==2);
id_triple = find(DATA(:,1)==3);

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

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

DATA_triple = [mean(DATA(id_triple,2:end))', std(DATA(id_triple,2:end))', ...
    min(DATA(id_triple,2:end))', max(DATA(id_triple,2:end))', median(DATA(id_triple,2:end))'];


%%
% TAE 
TAE = zeros(3,4);
TAE(1,:) = [mean(DATA(id_single,2)), std(DATA(id_single,2)) , ...
            min(DATA(id_single,2)),  max(DATA(id_single,2))]; 
        
TAE(2,:) = [mean(DATA(id_double,2)), std(DATA(id_double,2)) ...
            min(DATA(id_double,2)),  max(DATA(id_double,2))]; 

TAE(3,:) = [mean(DATA(id_triple,2)),  std(DATA(id_triple,2)) ... 
            min(DATA(id_triple,2)),  max(DATA(id_triple,2))]; 

% temperature variation of cooling phase 
COOLING_temp = zeros(3,4);
COOLING_temp(1,:) = [mean(DATA(id_single,3)), std(DATA(id_single,3)) , ...
                    min(DATA(id_single,3)),  max(DATA(id_single,3))];

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

COOLING_temp(3,:) = [mean(DATA(id_triple,3)),  std(DATA(id_triple,3)) ... 
                    min(DATA(id_triple,3)),  max(DATA(id_triple,3))]; 

% Duration cooling phae                
COOLING_time = zeros(3,4);
COOLING_time(1,:) = [mean(DATA(id_single,4)), std(DATA(id_single,4)) , ...
                    min(DATA(id_single,4)), max(DATA(id_single,4))];

COOLING_time(2,:) = [mean(DATA(id_double,4)), std(DATA(id_double,4)) ...
                    min(DATA(id_double,4)), max(DATA(id_double,4))];
                  
COOLING_time(3,:) = [mean(DATA(id_triple,4)), std(DATA(id_triple,4)) ... 
                    min(DATA(id_triple,4)), max(DATA(id_triple,4))]; 

% Cooling rate                
COOLING_rate = zeros(3,4);
COOLING_rate(1,:) = [mean(DATA(id_single,5)), std(DATA(id_single,5)) , ...
                    min(DATA(id_single,5)), max(DATA(id_single,5))];

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

COOLING_rate(3,:) = [mean(DATA(id_triple,5)), std(DATA(id_triple,5)) ... 
                    min(DATA(id_triple,5)), max(DATA(id_triple,5))]; 

% temperature variation of warming phase 
WARMING_temp = zeros(3,4);
WARMING_temp(1,:) = [mean(DATA(id_single,6)), std(DATA(id_single,6)) , ...
                    min(DATA(id_single,6)), max(DATA(id_single,6))];

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

WARMING_temp(3,:) = [mean(DATA(id_triple,6)), std(DATA(id_triple,6)) ... 
                    min(DATA(id_triple,6)), max(DATA(id_triple,6))]; 

                
WARMING_time = zeros(3,4);                
WARMING_time(1,:) = [mean(DATA(id_single,7)), std(DATA(id_single,7)) , ...
                    min(DATA(id_single,7)), max(DATA(id_single,7))];

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

WARMING_time(3,:) = [mean(DATA(id_triple,7)), std(DATA(id_triple,7)) ... 
                    min(DATA(id_triple,7)), max(DATA(id_triple,7))]; 

WARMING_rate = zeros(3,4);                
WARMING_rate(1,:) = [mean(DATA(id_single,8)), std(DATA(id_single,8)) , ...
                    min(DATA(id_single,8)), max(DATA(id_single,8))];

WARMING_rate(2,:) = [mean(DATA(id_double,8)), std(DATA(id_double,8)) ...
                    min(DATA(id_double,8)), max(DATA(id_double,8))];

WARMING_rate(3,:) = [mean(DATA(id_triple,8)), std(DATA(id_triple,8)) ... 
                    min(DATA(id_triple,8)), max(DATA(id_triple,8))]; 
                
%%
figure()
% TAE
subplot(2,2,1)
errorbar(TAE(:,1)./60, TAE(:,2)./60,'ok')
xlabel('Eruption order')
ylabel('TAE (min)')
% dT 
subplot(2,2,2)
errorbar(abs(COOLING_temp(:,1)), abs(COOLING_temp(:,2)),'bv')
hold on 
errorbar(WARMING_temp(:,1), WARMING_temp(:,2),'r^')
xlabel('Eruption order')
ylabel('\Delta T (ºC)')
% dt
subplot(2,2,3)
errorbar(COOLING_time(:,1)./60, COOLING_time(:,2)./60,'bv')
hold on 
errorbar(WARMING_time(:,1)./60, WARMING_time(:,2)./60,'r^')
xlabel('Eruption order')
ylabel('\Delta t (min)')
% rate
subplot(2,2,4)
errorbar(abs(COOLING_rate(:,1)).*60, abs(COOLING_rate(:,2)).*60,'bv')
hold on 
errorbar(WARMING_rate(:,1).*60, WARMING_rate(:,2).*60,'r^')
xlabel('Eruption order')
ylabel(' \Delta T/ \Delta t (ºC)')



%% Eruption selection
figure()
subplot(2,1,1)
plot(X22/60, dTs22)
hold on
plot(X22_e(sloc_Ts22)/60, dTs22_e(sloc_Ts22),'gv')
plot(X22_e(isloc_Ts22)/60, dTs22_e(isloc_Ts22),'k^')
plot(X22_e(mloc_Ts22)/60, dTs22_e(mloc_Ts22),'rv')
plot(X22_e(imloc_Ts22)/60, dTs22_e(imloc_Ts22),'b^')
xlim([0 170])
xlabel('Time(min)')
ylabel('Temperature variation (ºC)')
title('Eruption Night 22')
subplot(2,1,2)
plot(X23/60, dTs23)
hold on
plot(X23_1(mloc_Ts23_1)/60, dTs23_1(mloc_Ts23_1), 'rv')
plot(X23_1(imloc_Ts23_1)/60, dTs23_1(imloc_Ts23_1), 'b^')
plot(X23_1(sloc_Ts23_1)/60, dTs23_1(sloc_Ts23_1), 'gv')
plot(X23_1(isloc_Ts23_1)/60, dTs23_1(isloc_Ts23_1), 'k^')
plot(X23_2(mloc_Ts23_2)/60, dTs23_2(mloc_Ts23_2), 'rv')
plot(X23_2(imloc_Ts23_2)/60, dTs23_2(imloc_Ts23_2), 'b^')
plot(X23_2(sloc_Ts23_2)/60, dTs23_2(sloc_Ts23_2), 'gv')
plot(X23_2(isloc_Ts23_2)/60, dTs23_2(isloc_Ts23_2), 'k^')
xlim([0 255])
xlabel('Time(min)')
ylabel('Temperature variation (ºC)')
title('Eruption Night 23')

%% Post-processing 

