% This makes some diagnostic plots showing time-thickness tradeoffs. Idea
% being that you can fit the data pretty well with a longer period of
% thicker ice or a shorter period of thinner ice. Thus, it isn't
% possible to give precise dates for the lowstand. 
%
% 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; 
cm = jet(50);
nlines = 50;

msize = 5;

load current_results_H1

[sortM,sorti] = sort(results.Ms);
okH1 = sorti(1:nlines);
MH1 = sortM(1:nlines);
sorti2 = fliplr(sorti);

H1.time_lengths = results.times(:,3) - results.times(:,2);
H1.ice_thicknesses = (results.surface_elvs(:,2) + results.surface_elvs(:,3))./2 - 84.2;
H1.Ms = results.Ms;
H1.index = round(50.*results.Ms./400);
H1.index(find(H1.index > 50)) = zeros(size(find(H1.index > 50))) + 50;


figure('pos',[168 500 1290 300]);
subplot(1,3,1);

for a = 1:length(results.Ms)
    plot3(H1.time_lengths(sorti2(a)),H1.ice_thicknesses(sorti2(a)),H1.Ms(sorti2(a)),'o','markersize',msize,'color',cm(H1.index(sorti2(a)),:),'markerfacecolor',cm(H1.index(sorti2(a)),:)); hold on;
end

grid on; xl1 = xlabel(['Lowstand' newline 'duration (yr)']); yl1 = ylabel(['Mean ice thickness' newline 'during lowstand (m)']); zlabel('Misfit value');
set(yl1,'pos',[-3500 10 0],'horizontalalignment','center');
set(xl1,'pos',[4000 -5 0],'horizontalalignment','center');

axis([0 10000 0 20 0 500]);
set(gca,'cameraviewangle',10.8);
title('19-KP-H1');
drawnow;

%% 

load current_results_H4

[sortM,sorti] = sort(results.Ms);
okH4 = sorti(1:nlines);
MH4 = sortM(1:nlines);
sorti2 = fliplr(sorti);

H4.time_lengths = results.times(:,3) - results.times(:,2);
H4.ice_thicknesses = (results.surface_elvs(:,2) + results.surface_elvs(:,3))./2 - 80;
H4.Ms = results.Ms;
H4.index = round(50.*results.Ms./400);
H4.index(find(H4.index > 50)) = zeros(size(find(H4.index > 50))) + 50;


subplot(1,3,2);

for a = 1:length(results.Ms)
    plot3(H4.time_lengths(sorti2(a)),H4.ice_thicknesses(sorti2(a)),H4.Ms(sorti2(a)),'o','markersize',msize,'color',cm(H4.index(sorti2(a)),:),'markerfacecolor',cm(H4.index(sorti2(a)),:)); hold on;
end

grid on; xl4 = xlabel(['Lowstand' newline 'duration (yr)']); yl4 = ylabel(['Mean ice thickness' newline 'during lowstand (m)']); zlabel('Misfit value');
set(yl4,'pos',[-3500 10 0],'horizontalalignment','center');
set(xl4,'pos',[4000 -5 0],'horizontalalignment','center');

axis([0 10000 0 20 0 500]);
title('19-KP-H4');
set(gca,'cameraviewangle',10.8);
drawnow;


%% 

load current_results_H5

[sortM,sorti] = sort(results.Ms);
okH5 = sorti(1:nlines);
MH5 = sortM(1:nlines);
sorti2 = fliplr(sorti);

H5.time_lengths = results.times(:,3) - results.times(:,2);
H5.ice_thicknesses = (results.surface_elvs(:,2) + results.surface_elvs(:,3))./2 - 77.8;
H5.Ms = results.Ms;
H5.index = round(50.*results.Ms./400);
H5.index(find(H5.index > 50)) = zeros(size(find(H5.index > 50))) + 50;


subplot(1,3,3);

for a = 1:length(results.Ms)
    plot3(H5.time_lengths(sorti2(a)),H5.ice_thicknesses(sorti2(a)),H5.Ms(sorti2(a)),'o','markersize',msize','color',cm(H5.index(sorti2(a)),:),'markerfacecolor',cm(H5.index(sorti2(a)),:)); hold on;
end

axis([0 10000 0 20 0 500]);
grid on; xl5 = xlabel(['Lowstand' newline 'duration (yr)']); yl5 = ylabel(['Mean ice thickness' newline 'during lowstand (m)']); zlabel('Misfit value');
set(yl5,'pos',[-3500 10 0],'horizontalalignment','center');
set(xl5,'pos',[4000 -5 0],'horizontalalignment','center');

title('19-KP-H5');
set(gca,'cameraviewangle',10.8);
drawnow;


%% 

figure; 

plot([0 10000],[0.67 0.67],'-','color',[1 0.6 0.6],'linewidth',0.75); hold on;
plot([0 10000],[0.96 0.96],'-','color',[0.6 0.6 1],'linewidth',0.75);
plot([0 10000],[0.37 0.37],'-','color',[0.6 0.6 0.6],'linewidth',0.75);
text(9500,1.1,'Top of debris-rich ice in each core','horizontalalignment','right','verticalalignment','middle','background','w','fontsize',10);

% Choose whether to have a simple plot where there is just one symbol for
% each core, or a more complicated one where the symbol sizes are scaled
% inversely with M. 

if 1
    
    % This block scales the points by how good the fit is. Low value of M,
    % big symbol. Scale is arbitrary. 

    p2 = [];
    p3 = [];
    p1 = [];
    
    
    range1 = max(MH1) - min(MH1);
    range4 = max(MH4) - min(MH4);
    range5 = max(MH5) - min(MH5);
    
    size1 = 5 + 20.*(max(MH1) - MH1)./range1;
    size4 = 5 + 20.*(max(MH4) - MH4)./range4;
    size5 = 5 + 20.*(max(MH5) - MH5)./range5;

    for a = 1 :nlines
        p2(end+1) = plot(H4.time_lengths(okH4(a)),H4.ice_thicknesses(okH4(a)),'ks','markerfacecolor',[0.6 0.6 1],'markersize',size4(a));
        hold on;
    end

    for a = 1 :nlines
        p3(end+1) = plot(H5.time_lengths(okH5(a)),H5.ice_thicknesses(okH5(a)),'k^','markerfacecolor',[0.6 0.6 0.6],'markersize',size5(a));
        hold on;
    end

    for a = 1 :nlines
        p1(end+1) = plot(H1.time_lengths(okH1(a)),H1.ice_thicknesses(okH1(a)),'ko','markerfacecolor',[1 0.6 0.6],'markersize',size1(a));
        hold on;
    end
    
else
    % This just plots them all the same size. 

    p2(1) = plot(H4.time_lengths(okH4),H4.ice_thicknesses(okH4),'ks','markerfacecolor',[0.6 0.6 1],'markersize',9);
    hold on;
    p3(1) = plot(H5.time_lengths(okH5),H5.ice_thicknesses(okH5),'k^','markerfacecolor',[0.6 0.6 0.6],'markersize',9);
    p1(1) = plot(H1.time_lengths(okH1),H1.ice_thicknesses(okH1),'ko','markerfacecolor',[1 0.6 0.6],'markersize',9);
end


grid on
set(gca,'xlim',[0 10000],'ylim',[0 10],'fontsize',12);
xlabel('Duration of lowstand (yr)');
ylabel('Ice thickness at core site during lowstand (m)'); 
legend([p1(1) p2(1) p3(1)],'H1','H4','H5','location','northwest')





