%**  create plots of oxygen isotope records from sea-floor sediment cores
%    compare with computed insolation 
 
%** data sets
%  insolation_1Ma.mat  
%  INSOL_Jun    1001 x 8 array: June insolation at 8 latitudes, 
%                               1 ka interval
%                               calculated by Berger, A., 1992
%               col1: age, ka,  origin (0) is 1950 A.D.
%               col2: insolation at 90 N W m^-2
%               col3: insolation at 60 N W m^-2
%               col4: insolation at 30 N W m^-2
%               col5: insolation at 0    W m^-2
%               col6: insolation at 30 S W m^-2
%               col7: insolation at 60 S W m^-2
%               col8: insolation at 90 S W m^-2
%
%  INSOL_Dec    1001 x 8 array: December insolation, as above
%
%  insolation_5Ma.mat
%  INSOL_5Ma    5001 x 9        orbital parameters and insolation
%               col1: age, ka,  origin (0) is 1950 A.D. (radiocarbon age)
%               col2: eccentricity, ECC
%               col3: longitude of perihelion from moving vernal equinox 
%                      degree and decimals, OMEGA
%               col4: obliquity, degree and decimals
%               col5: climatic precession, ECC . SIN(OMEGA) 
%               col6: insolation 65N for July, W m^-2   (mid-month)
%               col7: insolation 65S for January, W m^-2  
%               col8: insolation 15N for July, W m^-2 
%               col9: insolation 15S for January, W m^-2 
%
%  Berger, A., 1992, Orbital Variations and Insolation Database. IGBP 
%  PAGES/World Data Center-A for Paleoclimatology Data Contribution Series 
%  # 92-007.  NOAA/NGDC Paleoclimatology Program, Boulder CO, USA.
%
%  Berger A. and M.F. Loutre, 1991, Insolation values for the climate of 
%  the last 10 million of years. Quaternary Sciences Review, 10 (4), 
%  297-317.
%
%  PliPleHol_s_delO18.mat
%  PliPleHol_s_delO18 2116 x 3  stacked benthic oxygen isotopes and errors
%                               age ka, delta O18 PDB per mil, error per mil
%
%  Lisiecki, L. E., and M. E. Raymo (2004), A Plio-Pleistocene stack of 
%  57 globally distributed benthic d18O records, Paleoceanography, 
%  doi:10.1029/2004PA001071.
%

clear
load insolation_1Ma.mat
load insolation_5Ma.mat
load PliPleHol_s_delO18.mat
 

figure(1)
clf
subplot(4,1,1)
axis([0 5300 2 4.5]), hold on
plot(PliPleHol_s_delO18(:,1), PliPleHol_s_delO18(:,2),'b-')
ylabel('\delta O^{18} per mil')
title('Plio-Pleistocene stacked global benthic delta O18')
set(gca,'XTick',[0:500:5300])
set(gca,'XGrid','on')
set(gca,'YTick',[-2:0.5:5])
set(gca,'YGrid','on')
 
subplot(4,1,2)
axis([0 2000 2 4.5]), hold on
plot(PliPleHol_s_delO18(1:1251,1), PliPleHol_s_delO18(1:1251,2),'r-')
ylabel('\delta O^{18} per mil')
set(gca,'XTick',[0:500:2000])
set(gca,'XGrid','on')
set(gca,'YTick',[-2:0.5:5])
set(gca,'YGrid','on')
 
subplot(4,1,3)
axis([2000 4000 2 4.5]), hold on
plot(PliPleHol_s_delO18(1251:1851,1), PliPleHol_s_delO18(1251:1851,2),'r-')
ylabel('\delta O^{18} per mil')
set(gca,'XTick',[2000:500:4000])
set(gca,'XGrid','on')
set(gca,'YTick',[-2:0.5:5])
set(gca,'YGrid','on')
 
subplot(4,1,4)
axis([4000 6000 2 4.5]), hold on
plot(PliPleHol_s_delO18(1851:2116,1), PliPleHol_s_delO18(1851:2116,2),'r-')
ylabel('\delta O^{18} per mil')
xlabel('thousands of years ago')
set(gca,'XTick',[4000:500:6000])
set(gca,'XGrid','on')
set(gca,'YTick',[-2:0.5:5])
set(gca,'YGrid','on')


%*****
% orbital parameters
%*****


figure(2)
clf
subplot(3,1,1)
axis([0 2000 22 24.5]), hold on
plot(INSOL_5Ma(1:2001,1), INSOL_5Ma(1:2001,4))
ylabel('obliquity (degrees)')
title('orbital paramters over last 2 Ma')

subplot(3,1,2)
axis([0 2000 0 0.06]), hold on
plot(INSOL_5Ma(1:2001,1), INSOL_5Ma(1:2001,2))
ylabel('eccentricity')

subplot(3,1,3)
axis([0 2000 -0.05 0.05]), hold on
plot(INSOL_5Ma(1:2001,1), INSOL_5Ma(1:2001,5))
ylabel({'climateic precession'; '(\epsilon sin(\omega))'})
xlabel('thousands of years ago')


 

%*****
% insolation
%*****

figure(3)
clf
subplot(2,1,1)
[AXa,H1a,H2a]=plotyy(INSOL_Jun(1:501,1), INSOL_Jun(1:501,3), INSOL_Dec(1:501,1), INSOL_Dec(1:501,7), 'plot');
set(get(AXa(1),'Ylabel'),'String','W m^-^2')
set(get(AXa(2),'Ylabel'),'String','W m^-^2')
set(AXa,'XGrid','on')
set(AXa,'XTick',[0:100:500])
set(AXa(1),'YLim',[425 575])
set(AXa(1),'YTick',[400:50:600])
set(AXa(2),'YLim',[425 575])
set(AXa(2),'YTick',[400:50:600])
set(gca,'YGrid','on')
title('summer insolation at 60 S and 60 N')

subplot(2,1,2)
[AXa,H1a,H2a]=plotyy(INSOL_Dec(1:501,1), INSOL_Dec(1:501,3), INSOL_Jun(1:501,1), INSOL_Jun(1:501,7), 'plot');
set(get(AXa(1),'Ylabel'),'String','W m^-^2')
set(get(AXa(2),'Ylabel'),'String','W m^-^2')
set(AXa,'XGrid','on')
set(AXa,'XTick',[0:100:500])
set(AXa(1),'YLim',[10 40])
set(AXa(1),'YTick',[0:10:50])
set(AXa(2),'YLim',[10 40])
set(AXa(2),'YTick',[0:10:50])
set(gca,'YGrid','on')
title('winter insolation at 60 S and 60 N')
xlabel('thousands of years ago')



figure(4)
clf
subplot(2,1,1)
[AX,H1,H2]=plotyy(INSOL_Jun(1:501,1), INSOL_Jun(1:501,3), PliPleHol_s_delO18(1:501,1), PliPleHol_s_delO18(1:501,2),'plot');
set(get(AX(1),'Ylabel'),'String','W m^-^2')
set(get(AX(2),'Ylabel'),'String','\delta O^{18} per mil')
set(AX,'XGrid','on')
set(AX,'XTick',[0:100:500])
set(AX(1),'YLim',[400 600])
set(AX(1),'YTick',[400:50:600])
set(AX(2),'YLim',[2.4 4.5])
set(AX(2),'YTick',[2.5:0.5:4.5])
title('June insolation at 60 N and global benthic del O18')
 
subplot(2,1,2)
[AXa,H1a,H2a]=plotyy(INSOL_Dec(1:501,1), INSOL_Dec(1:501,7),  PliPleHol_s_delO18(1:501,1),  PliPleHol_s_delO18(1:501,2), 'plot');
set(get(AXa(1),'Ylabel'),'String','W m^-^2')
set(get(AXa(2),'Ylabel'),'String','\delta O^{18} per mil')
set(AXa,'XGrid','on')
set(AXa,'XTick',[0:100:500])
set(AXa(1),'YLim',[400 600])
set(AXa(1),'YTick',[400:50:600])
set(AXa(2),'YLim',[2.4 4.5])
set(AXa(2),'YTick',[2.5:0.5:4.5])
title('December insolation at 60 S and global benthic \delta O^{18}')
xlabel('thousands of years ago')
 
 
figure(5)
clf
subplot(2,1,1)
[AX,H1,H2]=plotyy(INSOL_Dec(1:501,1), INSOL_Dec(1:501,3),  PliPleHol_s_delO18(1:501,1),  PliPleHol_s_delO18(1:501,2),'plot');
set(get(AX(1),'Ylabel'),'String','W m^-^2')
set(get(AX(2),'Ylabel'),'String','\delta O^{18} per mil')
set(AX,'XGrid','on')
set(AX,'XTick',[0:100:500])
set(AX(1),'YLim',[0 50])
set(AX(1),'YTick',[0:10:50])
set(AX(2),'YLim',[2.4 4.5])
set(AX(2),'YTick',[2.5:0.5:4.5])
title('December insolation at 60 N (note scale) and global benthic \delta O^{18}')
 
subplot(2,1,2)
[AXa,H1a,H2a]=plotyy(INSOL_Jun(1:501,1), INSOL_Jun(1:501,7),  PliPleHol_s_delO18(1:501,1),  PliPleHol_s_delO18(1:501,2), 'plot');
set(get(AXa(1),'Ylabel'),'String','W m^-^2')
set(get(AXa(2),'Ylabel'),'String','\delta O^{18} per mil')
set(AXa,'XGrid','on')
set(AXa,'XTick',[0:100:500])
set(AXa(1),'YLim',[0 50])
set(AXa(1),'YTick',[0:10:50])
set(AXa(2),'YLim',[2.4 4.5])
set(AXa(2),'YTick',[[2.5:0.5:4.5]])
title('June insolation at 60 S and global benthic \delta O^{18}')
xlabel('thousands of years ago')
 
 
figure(6)
clf
subplot(2,1,1)
[AX,H1,H2]=plotyy(PliPleHol_s_delO18(1:801,1), PliPleHol_s_delO18(1:801,2), INSOL_5Ma(1:1001,1), INSOL_5Ma(1:1001,7), 'plot');
set(get(AX(2),'Ylabel'),'String','W m^-^2')
set(get(AX(1),'Ylabel'),'String','\delta O^{18} per mil')
set(AX(1),'XGrid','on')
set(AX(1),'XLim',[0 1000])
set(AX(2),'XLim',[0 1000])
set(AX(1),'XTick',[0:100:1000])
set(AX(2),'YLim',[300 500])
set(AX(2),'YTick',[300:50:500])
set(AX(1),'YLim',[2.5 5])
set(AX(1),'YTick',[2:0.5:5])
title('January insolation at 65 S and global benthic \delta O^{18}')
 
subplot(2,1,2)
[AX,H1,H2]=plotyy(PliPleHol_s_delO18(801:1251,1), PliPleHol_s_delO18(801:1251,2), INSOL_5Ma(1001:2001,1), INSOL_5Ma(1001:2001,7), 'plot');
set(get(AX(2),'Ylabel'),'String','W m^-^2')
set(get(AX(1),'Ylabel'),'String','\delta O^{18} per mil')
set(AX(2),'XGrid','on')
set(AX(1),'XLim',[1000 2000])
set(AX(2),'XLim',[1000 2000])
set(AX(1),'XTick',[1000:100:2000])
set(AX(1),'YLim',[2.5 5])
set(AX(1),'YTick',[2:0.5:5])
set(AX(2),'YLim',[300 500])
set(AX(2),'YTick',[300:50:500])
%title('January insolation at 65 S and global benthic \delta O^{18}')
xlabel('thousands of years ago')
 
 
 
figure(7)
clf
subplot(2,1,1)
[AX,H1,H2]=plotyy(PliPleHol_s_delO18(1:801,1), PliPleHol_s_delO18(1:801,2), INSOL_5Ma(1:1001,1), INSOL_5Ma(1:1001,6), 'plot');
set(get(AX(2),'Ylabel'),'String','W m^-^2')
set(get(AX(1),'Ylabel'),'String','\delta O^{18} per mil')
set(AX(1),'XGrid','on')
set(AX(1),'XLim',[0 1000])
set(AX(2),'XLim',[0 1000])
set(AX(1),'XTick',[0:100:1000])
set(AX(2),'YLim',[300 500])
set(AX(2),'YTick',[300:50:500])
set(AX(1),'YLim',[2.5 5])
set(AX(1),'YTick',[2:0.5:5])
title('July insolation at 65 N and global benthic \delta O^{18}')
 
subplot(2,1,2)
[AX,H1,H2]=plotyy(PliPleHol_s_delO18(801:1251,1), PliPleHol_s_delO18(801:1251,2), INSOL_5Ma(1001:2001,1), INSOL_5Ma(1001:2001,6), 'plot');
set(get(AX(2),'Ylabel'),'String','W m^-^2')
set(get(AX(1),'Ylabel'),'String','\delta O^{18} per mil')
set(AX(2),'XGrid','on')
set(AX(1),'XLim',[1000 2000])
set(AX(2),'XLim',[1000 2000])
set(AX(1),'XTick',[1000:100:2000])
set(AX(1),'YLim',[2.5 5])
set(AX(1),'YTick',[2:0.5:5])
set(AX(2),'YLim',[300 500])
set(AX(2),'YTick',[300:50:500])
%title('July insolation at 65 N and global benthic \delta O^{18}')
xlabel('thousands of years ago')