%** identify terminations in delO18 time series
%
%   this script uses 
%       Matlab function findpeaks.m; must have signal processing toolbox
%
%       S. Waibel's f_localbndry.m
%       f_localbndry(x,h,flag)
%       x: dependent variable
%       h: half window size
%       flag: choose metric; Mahalanobis distance (D^2) or Euclidean (E^2)
%
%   this script is tailored to work with Lisieki and Raymo data set
%       PliPleHol_s_delO18.mat


clear
load PliPleHol_s_delO18.mat
load insolation_5Ma.mat

%*********
% prepare data
%*********

% extract time and time series values
t = PliPleHol_s_delO18(:,1);
z = PliPleHol_s_delO18(:,2);
clear PliPleHol_s_delO18

% segmenting the time series to handle different intervals a priori
% first segment;  1 kyr interval
t1 = t(1:601);
z1 = z(1:601);
% second segment;  2 kyr interval
t2 = t(601:1051);
z2 = z(601:1051);
% interpolate third segment to every 2 kyr 
t3initial = t(1051:1651);
z3initial = z(1051:1651);
t3 = [t(1051):2:t(1651)]';
z3 = interp1(t3initial, z3initial, t3);
% fourth segment;  5 kyr interval
t4 = t(1651:2116);
z4 = z(1651:2116);



%*********
% analysis and plots
%*********

%* FIRST segment 
h1 = 5;
[seg1_all seg1_pos] = f_localbndry(z1,h1,1);
% find peaks and location indices in E^2
[pks1,locs1] = findpeaks(seg1_pos,'MINPEAKHEIGHT',0.14);
t1bndry = t1(locs1);

figure(51)
clf
subplot(3,1,1)
plot(t1,z1)
ylabel('\delta O^{18} per mil')
set(gca,'XTick',t1(locs1))
set(gca,'XGrid','on')
title({'glacial terminations,'; 'benthic \delta O^{18}, and 65 N summer insolation'})

subplot(3,1,2)
plot(INSOL_5Ma(1:601,1),INSOL_5Ma(1:601,6))
hold on
for n=1:length(locs1)
    locs1p(n)=find(INSOL_5Ma(:,1)==t1(locs1(n))) ;
end
plot(t1(locs1),INSOL_5Ma(locs1p,6),'r.')
ylabel('W m^-^2')

subplot(3,1,3)
plot(t1,seg1_all,'r--')
hold on
plot(t1, seg1_pos,'r')
xlabel('thousands of years ago')
ylabel('E^2')   



%* SECOND segment
h2 = 10;
[seg2_all seg2_pos] = f_localbndry(z2,h2,0);
% find peaks and location indices in E^2
[pks2,locs2] = findpeaks(seg2_pos,'MINPEAKHEIGHT',0.26);
t2bndry = t2(locs2);

figure(52)
clf
subplot(3,1,1)
plot(t2,z2)
hold on
ylabel('\delta O^{18} per mil')
set(gca,'XTick',t2(locs2))
set(gca,'XGrid','on')
title({'glacial terminations,'; 'benthic \delta O^{18}, and 65 N summer insolation'})

subplot(3,1,2)
plot(INSOL_5Ma(601:1501,1),INSOL_5Ma(601:1501,6))
hold on
for n=1:length(locs2)
    locs2p(n)=find(INSOL_5Ma(:,1)==t2(locs2(n))) ;
end
plot(t2(locs2),INSOL_5Ma(locs2p,6),'r.')
ylabel('W m^-^2')

subplot(3,1,3)
plot(t2,seg2_all,'r--')
hold on
plot(t2,seg2_pos,'r')
xlabel('thousands of years ago')
ylabel('D^2')   


%* THIRD segment
h3 = 4;
[seg3_all seg3_pos] = f_localbndry(z3,h3,0);
% find peaks and location indices in E^2
[pks3,locs3] = findpeaks(seg3_pos,'MINPEAKHEIGHT',0.15,'minpeakdistance',10);
t3bndry = t3(locs3);

figure(53)
clf
subplot(3,1,1)
plot(t3,z3)
ylabel('\delta O^{18} per mil')
set(gca,'XTick',t3(locs3))
set(gca,'XGrid','on')
title({'glacial terminations,'; 'benthic \delta O^{18}, and 65 N summer insolation'})

subplot(3,1,2)
plot(INSOL_5Ma(1501:3001,1),INSOL_5Ma(1501:3001,6))
hold on
locs3p=zeros(size(locs3));
for n=1:length(locs3)
    locs3p(n)=find(INSOL_5Ma(:,1)==t3(locs3(n))) ;
end
plot(t3(locs3),INSOL_5Ma(locs3p,6),'r.')
ylabel('W m^-^2')

subplot(3,1,3)
plot(t3,seg3_all,'r--')
hold on
plot(t3,seg3_pos,'r')
xlabel('thousands of years ago')
ylabel('D^2')   



%* FOURTH segment
% h4 = 4;
% [seg4_all seg4_pos] = f_localbndry(z4,h4,0);
% % find peaks and location indices in D^2
% [pks4,locs4] = findpeaks(seg4_pos,'MINPEAKHEIGHT',0.03,'minpeakdistance',6);
% t4bndry = t4(locs4);
% 
% figure(54)
% clf
% subplot(3,1,1)
% axis([min(t4) max(t4) floor(min(z4)) ceil(max(z4))])
% hold on
% plot(t4,z4)
% ylabel('\delta O^{18} per mil')
% set(gca,'XTick',t4(locs4))
% set(gca,'XGrid','on')
% 
% subplot(3,1,2)
% axis([min(t4) max(t4) 400 500])
% hold on
% plot(INSOL_5Ma(3001:5001,1),INSOL_5Ma(3001:5001,6))
% N=length(find(t4(locs4)<5001));
% locs4p=zeros(1,N);   % insolation series ends at 5Ma
% for n=1:N
%     locs4p(n)=find(INSOL_5Ma(:,1)==t4(locs4(n))) ;
% end
% plot(t4(locs4(1:N)),INSOL_5Ma(locs4p,6),'r.')
% ylabel('W m^-^2')
% 
% subplot(3,1,3)
% plot(t4,seg4_all,'r--')
% hold on
% plot(t4,seg4_pos,'r')
% grid on
% xlabel('thousands of years ago')
% ylabel('D^2')