function plot_supplement_randomsearch_results(corenum)

% This just generates a figure showing the random search results,
% color-coded by M. Figure is used in supplement. 
%
% 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. 

% This plots supplemental figures showing all the results in a saved random
% search simulation. 

if nargin < 1
    corenum = 1; % select 1,4, or 5
end

% Load data. Various reassignments just to allow reusing other code. 


dH1.site_elv_now = 84.2;
dH1.surface_elv_now = 120;
dH1.dirty_ice_thickness_m = 0.67;


dH4.site_elv_now = 80;
dH4.surface_elv_now = 119.4;
dH4.dirty_ice_thickness_m = 0.96;


dH5.site_elv_now = 77.8;
dH5.surface_elv_now = 118.7;
dH5.dirty_ice_thickness_m = 0.37; 


load data_C14;
load data_Be10;
load data_Al26;

dH1.data10 = d10.H1;
dH1.data14 = d14.H1;
dH1.data26 = d26.H1;

dH4.data10 = d10.H4;
dH4.data14 = d14.H4;

dH5.data10 = d10.H5;
dH5.data14 = d14.H5;

% Now need to choose appropriate data set for core. 

eval(['thisd = dH' int2str(corenum) ';']);

%% Create the working figure

figure('pos',[487 538 1100 450]);
tic10 = [0 2000 4000 6000 8000];
ticz = [10:20:130];
fs = 12;
tic14 = [0 20000 40000 60000];

%% Now load random search results

eval(['load current_results_H' int2str(corenum)]);

%% Plot all data, in order

[sortM,sorti] = sort(results.Ms,'descend');

% colormap
%cm = winter(10);
%cm = cividis(10);


% This one is not bad, kind of blue-to-green     
 cm1 = [0.7071    0.7071    1.0000
    0.7454    0.7454    0.9718
    0.7817    0.7817    0.9428
    0.8165    0.8165    0.9129
    0.7222    0.7222    0.7778
    0.7222    0.7778    0.7222
    0.6667    0.7778    0.6667
    0.6000    0.7778    0.6000
    0.5000    0.7778    0.5000
    0.2    0.5    0.2];
     
% This one is light-to-dark green

cm2  = [    0.9686    0.9882    0.9608
    0.9059    0.9638    0.8876
    0.8065    0.9242    0.7808
    0.6810    0.8719    0.6562
    0.5333    0.8052    0.5272
    0.3660    0.7251    0.4192
    0.2157    0.6288    0.3333
    0.1068    0.5190    0.2488
         0    0.4096    0.1651
         0    0.2667    0.1059];
     
 
cm = cm1;
     
     





% Lowest 90% in gray
pct90 = round(length(sorti).*0.9);

for a = 1:pct90
    if corenum == 1
        subplot(1,4,1);
    else
        subplot(1,3,1);
    end
    plot(results.plot_N10s(sorti(a),:),results.plot_z_cm,'color',[0.9 0.9 0.9]); hold on;
    if corenum == 1
        subplot(1,4,3);
    else
        subplot(1,3,2);
    end
    plot(results.plot_N14s(sorti(a),:),results.plot_z_cm,'color',[0.9 0.9 0.9]); hold on;
    
    if corenum == 1
        subplot(1,4,4);
    else
        subplot(1,3,3);
    end
    plot(results.times(sorti(a),:),results.surface_elvs(sorti(a),:),'color',[0.9 0.9 0.9]); hold on;
    
    if corenum == 1
        subplot(1,4,2);
        plot(results.plot_N26s(sorti(a),:),results.plot_z_cm,'color',[0.9 0.9 0.9]); hold on;
    end
end
    
% Now top 10% 

for bin = 1:10
    starti = (bin-1).* (length(results.Ms) - pct90)./10 + 1;
    endi = starti + (length(results.Ms) - pct90)./10 -1;
    for a = (pct90 + starti):(pct90 + endi)
        if corenum == 1
            subplot(1,4,1);
        else
            subplot(1,3,1);
        end
        plot(results.plot_N10s(sorti(a),:),results.plot_z_cm,'color',cm(bin,:));
        
        if corenum == 1
            subplot(1,4,3);
        else
            subplot(1,3,2);
        end
        plot(results.plot_N14s(sorti(a),:),results.plot_z_cm,'color',cm(bin,:));
        
        if corenum == 1
            subplot(1,4,4);
        else
            subplot(1,3,3);
        end
        plot(results.times(sorti(a),:),results.surface_elvs(sorti(a),:),'color',cm(bin,:));
        
        if corenum == 1
            subplot(1,4,2);
            plot(results.plot_N26s(sorti(a),:),results.plot_z_cm,'color',cm(bin,:)); hold on;
        end
    end
end





%% Be-10 data

if corenum == 1
    subplot(1,4,1)
else
    subplot(1,3,1);
end
    
for a= 1:length(thisd.data10.z)
    % 95%
    xx = thisd.data10.Natomsg_bounds(a,[1 5]);
    yy = [1 1].*thisd.data10.z(a);
    plot(xx,yy,'-','linewidth',0.5,'color',[0.4 0.4 0.4]); hold on;
    % 68%
    xx = thisd.data10.Natomsg_bounds(a,[2 4]);
    yy = [1 1].*thisd.data10.z(a);
    plot(xx,yy,'k','linewidth',1); hold on;
end
plot(thisd.data10.Natomsg_bounds(:,3),thisd.data10.z,'ko','markerfacecolor',[1 1 1],'markersize',9,'linewidth',1);
set(gca,'ydir','reverse','ylim',[0 140],'xlim',[0 9000],'xtick',tic10); grid on;
ylabel('Depth in core (cm)'); xlabel('[^{10}Be] (atoms g^{-1})');
set(gca,'ytick',ticz,'fontsize',fs);

title(['19-KP-H' int2str(corenum)]);

%% C-14 data

if corenum == 1
    subplot(1,4,3);
else
    subplot(1,3,2);
end
    
for a= 1:length(thisd.data14.z)
    % 95%
    xx = thisd.data14.Natomsg_bounds(a,[1 5]);
    yy = [1 1].*thisd.data14.z(a);
    plot(xx,yy,'-','linewidth',0.5,'color',[0.4 0.4 0.4]); hold on;
    % 68%
    xx = thisd.data14.Natomsg_bounds(a,[2 4]);
    yy = [1 1].*thisd.data14.z(a);
    plot(xx,yy,'k','linewidth',1); hold on;
end
plot(thisd.data14.Natomsg_bounds(:,3),thisd.data14.z,'ko','markerfacecolor',[1 1 1],'markersize',9,'linewidth',1);
ylabel('Depth in core (cm)');xlabel('[^{14}C] (atoms g^{-1})');
set(gca,'ydir','reverse','ylim',[0 140],'xlim',[0 60000],'xtick',tic14); grid on;
set(gca,'ytick',ticz,'fontsize',fs);

if corenum == 1
    %% Al-26 data if present

    subplot(1,4,2);

    for a= 1:length(thisd.data26.z)
        % 95%
        xx = thisd.data26.Natomsg_bounds(a,[1 5]);
        yy = [1 1].*thisd.data26.z(a);
        plot(xx,yy,'-','linewidth',0.5,'color',[0.4 0.4 0.4]); hold on;
        % 68%
        xx = thisd.data26.Natomsg_bounds(a,[2 4]);
        yy = [1 1].*thisd.data26.z(a);
        plot(xx,yy,'k','linewidth',1); hold on;
    end
    plot(thisd.data26.Natomsg_bounds(:,3),thisd.data26.z,'k^','markerfacecolor',[1 1 1],'markersize',9,'linewidth',1);
    ylabel('Depth in core (cm)');xlabel('[^{26}Al] (atoms g^{-1})');
    set(gca,'ydir','reverse','ylim',[0 140],'xlim',[0 60000],'xtick',tic14); grid on;
    set(gca,'ytick',ticz,'fontsize',fs);
end

%% Time-elevation axis

if corenum == 1
    subplot(1,4,4)
else
    subplot(1,3,3);
end
    
plot([0 15000],[1 1].*thisd.surface_elv_now,'k'); hold on;
plot([0 15000],[1 1].*thisd.site_elv_now,'k--');
set(gca,'xlim',[0 10000],'ylim',[70 140]);

xlabel('Age (yr)'); ylabel('Elevation (m)'); grid on;

p1 = plot([0 15000],[1 1].*thisd.surface_elv_now,'k'); hold on;
p2 = plot([0 15000],[1 1].*thisd.site_elv_now,'k--');

legend([p1 p2],'Present ice surface','Core top','location','southwest');


colorbar
colormap(cm)
caxis([90 100]);




