% This generates figure showing all data from Kay Peak ridge and subglacial 
% bedrock cores, with results of random search optimization. This is Figure
% 3 in the paper. 
%
% Greg Balco
% Berkeley Geochronology Center
% Feb, 2022
% 
% Accompanies a paper entitled 'REVERSIBLE ICE SHEET THINNING IN THE
% AMUNDSEN SEA EMBAYMENT DURING THE LATE HOLOCENE,' by Balco, Brown,
% Nichols, Venturelli, and others. 
%
% This code is provided solely for review of the accompanying paper. It
% is not licensed for other use or distribution. 

clear all; %close all;

%% Set up axes

figure('pos',[800 300 600 700],'tag','modeldatafigure');

raxw = 0.36;
laxw = 0.215;
lm = 0.08;
rm = 0.03;
thin = 0.01;
axhsp = (1 - raxw - laxw -laxw - lm - rm - thin);

axht = 0.285;
bm = 0.08;

fs = 12; % fontsize global

% colors global
col14 = [0.6 0.6 0.9];
%col14 = [166 206 227]./255;
col10 = [0.9 0.6 0.6];
%col10 = [31 120 180]./255;

colnull = [51 160 44]./244;
colnull = 'k';

lcol = [0.4 0.7 0.4];
%lcol = [178 223 138]./255;

% linewidth global
lw68 = 1.5;
lw95 = 0.5;

ax1l = axes('pos',[lm (bm+2*axht+2*thin) laxw axht]); hold on;

ax1m = axes('pos',[lm+laxw+thin (bm+2*axht+2*thin) laxw axht]); hold on;
ax1r = axes('pos',[lm+laxw+thin+laxw+axhsp (bm+2*axht+2*thin) raxw axht]); hold on;

ax4l = axes('pos',[lm (bm+axht+thin) laxw axht]); hold on;
ax4m = axes('pos',[lm+laxw+thin (bm+axht+thin) laxw axht]); hold on;
ax4r = axes('pos',[lm+laxw+thin+laxw+axhsp (bm+axht+thin) raxw axht]); hold on;

ax5l = axes('pos',[lm (bm) laxw axht]); hold on;
ax5m = axes('pos',[lm+laxw+thin (bm) laxw axht]); hold on;
ax5r = axes('pos',[lm+laxw+thin+laxw+axhsp (bm) raxw axht]); hold on; 

%% Findability

set(ax1l,'tag','ax1l');
set(ax1m,'tag','ax1m');
set(ax1r,'tag','ax1r');

set(ax4l,'tag','ax4l');
set(ax4m,'tag','ax4m');
set(ax4r,'tag','ax4r');

set(ax5l,'tag','ax5l');
set(ax5m,'tag','ax5m');
set(ax5r,'tag','ax5r');


%% 

plot_randomsearch_bad_results()

%% Plot data on axes

% 1. Plot core measurements

tic10 = [0 2000 4000 6000 8000];
ticz = [10:20:130];

% Be-10

load data_Be10.mat

axes(ax1l);
% Observed Be-10
for a= 1:length(d10.H1.z)
    % 95%
    xx = d10.H1.Natomsg_bounds(a,[1 5]);
    yy = [1 1].*d10.H1.z(a);
    plot(xx,yy,'-','linewidth',lw95,'color',[0.4 0.4 0.4]); hold on;
    % 68%
    xx = d10.H1.Natomsg_bounds(a,[2 4]);
    yy = [1 1].*d10.H1.z(a);
    plot(xx,yy,'k','linewidth',lw68); hold on;
end
plot(d10.H1.Natomsg_bounds(:,3),d10.H1.z,'ks','markerfacecolor',col10,'markersize',8);
set(gca,'ydir','reverse','ylim',[0 140],'xlim',[0 9000],'xtick',tic10,'xticklabel',{}); grid on;
ylabel('Depth in core (cm)');
set(gca,'ytick',ticz,'fontsize',fs);

axes(ax4l);
% Observed Be-10
for a= 1:length(d10.H4.z)
    % 95%
    xx = d10.H4.Natomsg_bounds(a,[1 5]);
    yy = [1 1].*d10.H4.z(a);
    plot(xx,yy,'-','linewidth',lw95,'color',[0.4 0.4 0.4]); hold on;
    % 68%
    xx = d10.H4.Natomsg_bounds(a,[2 4]);
    yy = [1 1].*d10.H4.z(a);
    plot(xx,yy,'k','linewidth',lw68); hold on;
end
plot(d10.H4.Natomsg_bounds(:,3),d10.H4.z,'ks','markerfacecolor',col10,'markersize',8);
set(gca,'ydir','reverse','ylim',[0 140],'xlim',[0 9000],'xtick',tic10,'xticklabel',{}); grid on;
ylabel('Depth in core (cm)');
set(gca,'ytick',ticz,'fontsize',fs);

axes(ax5l);
% Observed Be-10
for a= 1:length(d10.H5.z)
    % 95%
    xx = d10.H5.Natomsg_bounds(a,[1 5]);
    yy = [1 1].*d10.H5.z(a);
    plot(xx,yy,'-','linewidth',lw95,'color',[0.4 0.4 0.4]); hold on;
    % 68%
    xx = d10.H5.Natomsg_bounds(a,[2 4]);
    yy = [1 1].*d10.H5.z(a);
    plot(xx,yy,'k','linewidth',lw68); hold on;
end
plot(d10.H5.Natomsg_bounds(:,3),d10.H5.z,'ks','markerfacecolor',col10,'markersize',8);
set(gca,'ydir','reverse','ylim',[0 140],'xlim',[0 9000],'xtick',tic10,'xticklabel',{'0','2','4','6','8'}); grid on;
xlabel('[^{10}Be] (katoms g^{-1})');
ylabel('Depth in core (cm)');
set(gca,'ytick',ticz,'fontsize',fs);

%% Now plot C-14 data

load data_C14.mat
tic14 = [0 20000 40000 60000];

axes(ax1m);
for a= 1:length(d14.H1.z)
    % 95%
    xx = d14.H1.Natomsg_bounds(a,[1 5]);
    yy = [1 1].*d14.H1.z(a);
    plot(xx,yy,'-','linewidth',lw95,'color',[0.4 0.4 0.4]); hold on;
    % 68%
    xx = d14.H1.Natomsg_bounds(a,[2 4]);
    yy = [1 1].*d14.H1.z(a);
    plot(xx,yy,'k','linewidth',lw68); hold on;
end
plot(d14.H1.Natomsg_bounds(:,3),d14.H1.z,'ko','markerfacecolor',col14,'markersize',8);
set(gca,'ydir','reverse','ylim',[0 140],'xlim',[0 60000],'xtick',tic14,'xticklabel',{}); grid on;
set(gca,'ytick',ticz,'yticklabel',{},'fontsize',fs);

text(58000,45,'19-KP-H1','fontsize',12,'fontweight','bold','horizontalalignment','right','verticalalignment','middle','background','w');

axes(ax4m);
for a= 1:length(d14.H4.z)
    % 95%
    xx = d14.H4.Natomsg_bounds(a,[1 5]);
    yy = [1 1].*d14.H4.z(a);
    plot(xx,yy,'-','linewidth',lw95,'color',[0.4 0.4 0.4]); hold on;
    % 68%
    xx = d14.H4.Natomsg_bounds(a,[2 4]);
    yy = [1 1].*d14.H4.z(a);
    plot(xx,yy,'k','linewidth',lw68); hold on;
end
plot(d14.H4.Natomsg_bounds(:,3),d14.H4.z,'ko','markerfacecolor',col14,'markersize',8);
set(gca,'ydir','reverse','ylim',[0 140],'xlim',[0 60000],'xtick',tic14,'xticklabel',{}); grid on;
set(gca,'ytick',ticz,'yticklabel',{},'fontsize',fs);
text(58000,10,'19-KP-H4','fontsize',12,'fontweight','bold','horizontalalignment','right','verticalalignment','middle','background','w');


axes(ax5m);
for a= 1:length(d14.H5.z)
    % 95%
    xx = d14.H5.Natomsg_bounds(a,[1 5]);
    yy = [1 1].*d14.H5.z(a);
    plot(xx,yy,'-','linewidth',lw95,'color',[0.4 0.4 0.4]); hold on;
    % 68%
    xx = d14.H5.Natomsg_bounds(a,[2 4]);
    yy = [1 1].*d14.H5.z(a);
    plot(xx,yy,'k','linewidth',lw68); hold on;
end
plot(d14.H5.Natomsg_bounds(:,3),d14.H5.z,'ko','markerfacecolor',col14,'markersize',8);
set(gca,'ydir','reverse','ylim',[0 140],'xlim',[0 60000],'xtick',tic14); grid on;
set(gca,'ytick',ticz,'yticklabel',{},'fontsize',fs);
xlabel('[^{14}C] (katoms g^{-1})');
set(gca,'xticklabel',{'0','20','40','60'});
text(58000, 10,'19-KP-H5','fontsize',12,'fontweight','bold','horizontalalignment','right','verticalalignment','middle','background','w');



%% Now plot exposure-age data

tc = [0.5 0.5 0.7];

load data_exposure_ages.mat

%% H1 R axis

axes(ax1r);

% plot([1000 7500],[140 140],'--','color',tc); hold on;
% text(7600,140,'Lowest exposed rock','fontsize',10,'horizontalalignment','left','verticalalignment','middle','background','w','fontangle','italic','color',tc);

plot([800 9800],[140 140],'-','color',tc); hold on;
text(9900,140,'2020','fontsize',10,'horizontalalignment','left','verticalalignment','middle','background','w','fontangle','italic','color',tc);

plot([800 9800],[171 171],'-','color',tc); hold on;
text(9900,171,'1966','fontsize',10,'horizontalalignment','left','verticalalignment','middle','background','w','fontangle','italic','color',tc);

plot([500 9000],[120 120],'--','color',tc);
text(9100,120,'H1 ice surface','fontsize',10,'horizontalalignment','left','verticalalignment','middle','background','w','fontangle','italic','color',tc);

plot([0 9000],[84.5 84.5],'-','color',tc);
text(9100,84.5,'H1 core top','fontsize',10,'horizontalalignment','left','verticalalignment','middle','background','w','fontangle','italic','color',tc);




for a = 1:length(t14.elv_m)
    % 95%
    xx = t14.t_St_bounds(a,[1 5]);
    yy = [1 1].*t14.elv_m(a);
    plot(xx,yy,'-','linewidth',lw95,'color',[0.4 0.4 0.4]); hold on;
    % 68%
    xx = t14.t_St_bounds(a,[2 4]);
    yy = [1 1].*t14.elv_m(a);
    plot(xx,yy,'k','linewidth',lw68); hold on;
end
plot(t14.t_St_bounds(:,3),t14.elv_m,'ko','markerfacecolor',col14,'markersize',7);


for a = 1:length(t10.elv_m)
    % 95%
    xx = [2 -2].*t10.dtint_St(a) + t10.t_St(a);
    yy = [1 1].*t10.elv_m(a);
    plot(xx,yy,'-','linewidth',lw95,'color',[0.4 0.4 0.4]); hold on;
    % 68%
    xx = [1 -1].*t10.dtint_St(a) + t10.t_St(a);
    yy = [1 1].*t10.elv_m(a);
    plot(xx,yy,'k','linewidth',lw68); hold on;
end
plot(t10.t_St,t10.elv_m,'ks','markerfacecolor',col10,'markersize',7);


set(gca,'xlim',[0 14000],'xtick',[0:2000:14000],'xticklabel',{});
set(gca,'ylim',[60 360],'ytick',[100 200 300],'fontsize',fs);
ylabel('Elevation (m)');



grid on;


%% H4 R axis

axes(ax4r);

plot([800 9800],[140 140],'-','color',tc); hold on;
text(9900,140,'2020','fontsize',10,'horizontalalignment','left','verticalalignment','middle','background','w','fontangle','italic','color',tc);

plot([800 9800],[171 171],'-','color',tc); hold on;
text(9900,171,'1966','fontsize',10,'horizontalalignment','left','verticalalignment','middle','background','w','fontangle','italic','color',tc);

plot([500 9000],[119.4 119.4],'--','color',tc);
text(9100,119.4,'H4 ice surface','fontsize',10,'horizontalalignment','left','verticalalignment','middle','background','w','fontangle','italic','color',tc);

plot([0 9000],[80 80],'-','color',tc);
text(9100,80,'H4 core top','fontsize',10,'horizontalalignment','left','verticalalignment','middle','background','w','fontangle','italic','color',tc);


plot([800 800],[140 190],'-','color',tc);

text(500,190,['Lowest' newline 'exposed' newline 'rock'],'horizontalalignment','left','verticalalignment','bottom','background','w','fontangle','italic','color',tc);



for a = 1:length(t14.elv_m)
    % 95%
    xx = t14.t_St_bounds(a,[1 5]);
    yy = [1 1].*t14.elv_m(a);
    plot(xx,yy,'-','linewidth',lw95,'color',[0.4 0.4 0.4]); hold on;
    % 68%
    xx = t14.t_St_bounds(a,[2 4]);
    yy = [1 1].*t14.elv_m(a);
    plot(xx,yy,'k','linewidth',lw68); hold on;
end
plot(t14.t_St_bounds(:,3),t14.elv_m,'ko','markerfacecolor',col14,'markersize',7);


for a = 1:length(t10.elv_m)
    % 95%
    xx = [2 -2].*t10.dtint_St(a) + t10.t_St(a);
    yy = [1 1].*t10.elv_m(a);
    plot(xx,yy,'-','linewidth',lw95,'color',[0.4 0.4 0.4]); hold on;
    % 68%
    xx = [1 -1].*t10.dtint_St(a) + t10.t_St(a);
    yy = [1 1].*t10.elv_m(a);
    plot(xx,yy,'k','linewidth',lw68); hold on;
end
plot(t10.t_St,t10.elv_m,'ks','markerfacecolor',col10,'markersize',7);


set(gca,'xlim',[0 14000],'xtick',[0:2000:14000],'xticklabel',{});
set(gca,'ylim',[60 360],'ytick',[100 200 300],'fontsize',fs);
ylabel('Elevation (m)');

grid on;

%% H5 R axis

axes(ax5r);

plot([800 9800],[140 140],'-','color',tc); hold on;
text(9900,140,'2020','fontsize',10,'horizontalalignment','left','verticalalignment','middle','background','w','fontangle','italic','color',tc);

plot([800 9800],[171 171],'-','color',tc); hold on;
text(9900,171,'1966','fontsize',10,'horizontalalignment','left','verticalalignment','middle','background','w','fontangle','italic','color',tc);

plot([500 9000],[118.7 118.7],'--','color',tc);
text(9100,118.7,'H5 ice surface','fontsize',10,'horizontalalignment','left','verticalalignment','middle','background','w','fontangle','italic','color',tc);

plot([0 9000],[77.8 77.8],'-','color',tc);
text(9100,77.8,'H5 core top','fontsize',10,'horizontalalignment','left','verticalalignment','middle','background','w','fontangle','italic','color',tc);

for a = 1:length(t14.elv_m)
    % 95%
    xx = t14.t_St_bounds(a,[1 5]);
    yy = [1 1].*t14.elv_m(a);
    plot(xx,yy,'-','linewidth',lw95,'color',[0.4 0.4 0.4]); hold on;
    % 68%
    xx = t14.t_St_bounds(a,[2 4]);
    yy = [1 1].*t14.elv_m(a);
    plot(xx,yy,'k','linewidth',lw68); hold on;
end
plot(t14.t_St_bounds(:,3),t14.elv_m,'ko','markerfacecolor',col14,'markersize',7);


for a = 1:length(t10.elv_m)
    % 95%
    xx = [2 -2].*t10.dtint_St(a) + t10.t_St(a);
    yy = [1 1].*t10.elv_m(a);
    plot(xx,yy,'-','linewidth',lw95,'color',[0.4 0.4 0.4]); hold on;
    % 68%
    xx = [1 -1].*t10.dtint_St(a) + t10.t_St(a);
    yy = [1 1].*t10.elv_m(a);
    plot(xx,yy,'k','linewidth',lw68); hold on;
end
plot(t10.t_St,t10.elv_m,'ks','markerfacecolor',col10,'markersize',7);


set(gca,'xlim',[0 14000],'xtick',[0:2000:14000],'xticklabel',{'0','2','4','6','8','10','12','14'});
set(gca,'ylim',[60 360],'ytick',[100 200 300],'fontsize',fs);
ylabel('Elevation (m)');
xlabel('Age (kyr)');

grid on;

%% 

plot_randomsearch_results(lcol)

%% Globals


set(get(gcf,'children'),'fontname','gill sans','fontangle','normal','fontweight','normal')
set(findobj('type','axes'),'box','on')

%% Annotations

allaxtop = bm + 3 * axht + 2 * thin;

al1 = annotation('line',[0.082 0.35],[0.015 0.015]+allaxtop,'color','r');
al2 = annotation('line',[0.082 0.082],[0.005 0.015]+allaxtop,'color','r');
al3 = annotation('line',[0.313 0.313],[0.005 0.015]+allaxtop,'color','r');
set([al1 al2 al3],'color',colnull);

axes(ax1l);
at1 = text(11500,-8,'Concentrations predicted by null hypotheses','verticalalignment','middle','fontangle','italic','color',colnull,'fontsize',9);

axes(ax1r);
plot([8300 9300],[260 240],'-','color',colnull);
text(8450,235,'Null hypotheses','verticalalignment','top','fontangle','italic','color',colnull,'fontsize',9,'background','w');

plot(800,340,'ko','markerfacecolor',col14,'markersize',7);
text(1300,340,'^{14}C','fontsize',9,'verticalalignment','middle','horizontalalignment','left','background','w');
plot(800,317,'ks','markerfacecolor',col10,'markersize',7);
text(1300,317,'^{10}Be','fontsize',9,'verticalalignment','middle','background','w','horizontalalignment','left');

abcup = 0.035;

 


al4 = annotation('line',[lm+laxw.*0.55 lm+laxw],[abcup abcup]+allaxtop,'color','k');
al5 = annotation('line',[lm lm+laxw.*0.45],[abcup abcup]+allaxtop,'color','k');

am4 = annotation('line',[lm+laxw+thin lm+laxw+thin+laxw.*0.45],[abcup abcup]+allaxtop,'color','k');
am4 = annotation('line',[lm+laxw+thin+laxw.*0.55 lm+laxw+thin+laxw],[abcup abcup]+allaxtop,'color','k');

ar4 = annotation('line',[lm+laxw+thin+laxw+axhsp lm+laxw+thin+laxw+axhsp+raxw.*0.46],[abcup abcup]+allaxtop,'color','k');
ar4 = annotation('line',[lm+laxw+thin+laxw+axhsp+raxw.*0.54 lm+laxw+thin+laxw+axhsp+raxw],[abcup abcup]+allaxtop,'color','k');


% get coord for 'A'
diff_rel = abcup./axht;
diff_plot = diff(get(ax1l,'ylim')).*diff_rel;
mid_plot = sum(get(ax1l,'xlim'))./2;
axes(ax1l); 
atl1 = text(mid_plot,-diff_plot,'A','fontsize',14,'fontweight','bold','color','k','horizontalalignment','center','verticalalignment','middle');

% for 'B'
diff_plot = diff(get(ax1m,'ylim')).*diff_rel;
mid_plot = sum(get(ax1m,'xlim'))./2;
axes(ax1m); 
atm1 = text(mid_plot,-diff_plot,'B','fontsize',14,'fontweight','bold','color','k','horizontalalignment','center','verticalalignment','middle');

% for 'C'
diff_plot = diff(get(ax1r,'ylim')).*diff_rel;
mid_plot = sum(get(ax1r,'xlim'))./2;
axes(ax1r); 
atm1 = text(mid_plot,max(get(ax1r,'ylim'))+diff_plot,'C','fontsize',14,'fontweight','bold','color','k','horizontalalignment','center','verticalalignment','middle');



%% Plot null hypotheses

plot_null_hypotheses(colnull)

