function random_search_wrapper(corenum)

% This is the wrapper for the random search algorithm.
%
% Procedure is as follows: 
%  -- At top of this script, select core (1,4, or 5) and number of trials
%  -- Run this. This saves a file with the results. 
%  -- 'generate_model_data_figure.m' and some other plotting scripts
%       can then plot these results as needed.
%  
% Note: some core-specific data, specifically present bedrock and surface
% elevations, are defined in this script and also in some of the plotting 
% scripts, including, for example, but not limited to, 
% 'plot_null_hypotheses.m'. Must change these in all locations. See notes
% on this below. 
%
% Most data are loaded from files generated by initialization scripts. 
%
% 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. 

% Select core 

if nargin < 1 
    corenum = 1;
end

% Number of iterations
numits = 2000;

% Assemble input data. All these things have to be passed in a structure. 

% d.data10 - data structure for N10 
% d.data14 - data structure for N14 
% d.elv_now - meters present ice surface elevation
% d.t1 - output of generate_te_history
% d.l14 - C-14 decay const
% d.l10 - Be-10 decay const
% d.Dmass, d.Lmass - ice density approx parameters
% d.rhoice - density of glacier ice 

% Data loaded from preexisting files

load KP_consts
d.l14 = consts.l14;
d.l10 = consts.l10;
d.l26 = consts.l26;
d.rhoice = consts.rhoice;

load KP_calculated_P
d.P = P;

load icedensity_params
d.Lmass = Lmass;

% Other constants. Note: N10inh is also defined in plot_null_hypotheses.m. 

d.N10inh = 0; d.N26inh = 0;
d.plot_z_cm = 5:5:135; % This is used for plotting

% Core-specific values
%
% Note: some of these values are redundantly defined in 
% plot_null_hypotheses.m. Must change in both places. Bad coding!

dH1 = d;
dH1.Dmass = Dmasses.H1;
dH1.site_elv_now = 84.2;
dH1.surface_elv_now = 120;
dH1.dirty_ice_thickness_m = 0.67;

dH4 = d;
dH4.Dmass = Dmasses.H4;
dH4.site_elv_now = 80;
dH4.surface_elv_now = 119.4;
dH4.dirty_ice_thickness_m = 0.96;

dH5 = d;
dH5.Dmass = Dmasses.H5;
dH5.site_elv_now = 77.8;
dH5.surface_elv_now = 118.7;
dH5.dirty_ice_thickness_m = 0.37; 

% Load data. Note: uses confidence interval form. 

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];

% Note: if H1, need to add axes for Al-26 data

% 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 0.5 0.5],'markersize',9);
set(gca,'ydir','reverse','ylim',[0 140],'xlim',[0 9000],'xtick',tic10,'xticklabel',{}); 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)]);


if corenum == 1
   
    % Plot Al-26 data for H1
    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',[0.5 1 0.5],'markersize',9);
    set(gca,'ydir','reverse','ylim',[0 140],'xlim',[0 60000],'xtick',[0 20000 40000 60000],'xticklabel',{}); grid on;
    ylabel('Depth in core (cm)'); xlabel('[^{26}Al] (atoms g^{-1})');
    set(gca,'ytick',ticz,'fontsize',fs);

end
    
    
    

% 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',[0.5 0.5 1],'markersize',9);
set(gca,'ydir','reverse','ylim',[0 140],'xlim',[0 9000],'xtick',tic10,'xticklabel',{}); grid on;
ylabel('Depth in core (cm)');xlabel('[^{14}C] (atoms g^{-1})');
set(gca,'ydir','reverse','ylim',[0 140],'xlim',[0 60000],'xtick',tic14,'xticklabel',{}); grid on;
set(gca,'ytick',ticz,'fontsize',fs);

% 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 15000],'ylim',[60 360]);

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


%% Generate random search results

results.plot_z_cm = d.plot_z_cm;

for a = 1:numits
    
    d = thisd; % start fresh
    
    % Generate random time-elevation history
    % Note: dirty ice thickness needs to be an added mass thickness. 
    d.t1 = generate_random_Et(thisd.surface_elv_now,thisd.site_elv_now,thisd.dirty_ice_thickness_m); 
    
    % Plot time/elevation on RH plots
    if corenum == 1
        subplot(1,4,4)
    else
        subplot(1,3,3);
    end
    
    current_t_plot = plot(d.t1.times,d.t1.surface_elvs,'r'); hold on; drawnow;
    
    % Run forward model
    current_model_result = trial_forward_integration(d,1);
    
    if current_model_result.M < 150 % random value to highlight decent results
        thiscol = 'g';  
    else
        thiscol = [0.8 0.8 0.8];
    end
    
    set(current_t_plot,'color',thiscol);
    if corenum == 1
        subplot(1,4,1);
    else
        subplot(1,3,1);
    end
    plot(current_model_result.plot_N10,d.plot_z_cm,'-','color',thiscol);
    
    if corenum == 1
        subplot(1,4,3);
    else
        subplot(1,3,2);
    end
    plot(current_model_result.plot_N14,d.plot_z_cm,'-','color',thiscol);
    
    
    if corenum == 1
        subplot(1,4,2)
        plot(current_model_result.plot_N26,d.plot_z_cm,'-','color',thiscol);
    end
    
    drawnow
    
    results.Ms(a) = current_model_result.M;
    results.plot_N10s(a,:) = current_model_result.plot_N10;
    results.plot_N14s(a,:) = current_model_result.plot_N14;
    if corenum == 1
        results.plot_N26s(a,:) = current_model_result.plot_N26;
    end
    results.times(a,:) = d.t1.times;
    results.surface_elvs(a,:) = d.t1.surface_elvs;
     
    disp(int2str(a));
end
    

% Save results

eval(['save current_results_H' int2str(corenum) ' results']);

