% This script runs all calculations for random search algorithm to find
% elevation change histories that match Kay Peak Ridge core data. 
%
% 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;

%% First, run all initialization scripts. 

% This runs all initialize scripts necessary to do Kay Peak forward model
% and random search calculations. If you change the data or parameters in
% the initialization script, you have to run it again to make the new data
% available to other functions. 

% 1. Be-10 and Al-26 data - 
% writes files called 'data_Be10.mat' and 'data_Al26.mat'
%
% Be-10 data source: paste Be-10 concentration data from Imperial lab
% spreadsheet/supplemental table into 'init_Be10_data.m'. 
%
% Al-26 data source: paste Al-26 concentration data from Imperial lab
% spreadsheet/supplemental table into 'init_Al26_data.m'. 


disp('Generating Be-10 concentrations data file');
init_Be10_data()

disp('Generating Al-26 concentrations data file');
init_Al26_data()

% -------

% 2. C-14 data - writes files:
%   'data_cores_C14.mat' -- C-14 concentrations from core samples
%          -- also includes depth in core
%   'data_CRONUSA_C14.mat' -- C-14 concentrations from contemporaneous CRONUS-A
%   'data_surface_C14.mat' -- C-14 concentrations from ridge samples
%
% Data source: paste C-14 data from Tulane lab spreadsheet/supplemental 
% table  into 'init_C14_data.m'. 
%
% Note: blank corrections for C-14 data are in this script. 

disp('Generating C-14 concentrations data files');
init_C14_data()

% -------

% 3. Ice/firn density approximation parameters
%
% Writes 'icedensity_params.mat' with parameters for equation approximating
% firn/ice mass depth - depth relation
%
% Data source: init_data/KP_mass_depths_202011.mat, which contains mass
% depths derived from integrating field firn density measurements.

disp('Generating ice/firn mass depth fitting parameters')
init_ice_density_fitting()

% -------

% 4. Production rate precalculation
%
% First, generates constants file from 'init_make_consts_KP.m'. Then
% precalculates Be-10 and C-14 production rates on a grid in surface
% elevation and mass depth below the surface. 
%
% Uses various code from the online exposure age calculators. 
%
% Results in two file writes: 'KP_consts.mat' and 'KP_calculated_P.mat'
% 
% Data source: parameters in 'init_make_consts_KP.mat'

disp('Precalculating production rates')
init_precalculate_P10_P14()

% -------

% 5. Calculate exposure ages for surface samples at Kay Peak ridge.
%
% Data sources:
%
%   -- ICE-D data for previously published Be-10 measurments. Not obtained
%   dynamically, text input format is included in
%   'init_get_surface_exposure_ages.m'. 
%
%  -- ICE-D sample and site data for all previously published and 
%  newly reported. C-14 samples. Likewise, not obtained dynamically, 
%  included in m-file. Note: does not use C-14 concentrations from 
%  database - gets them from file written by init_C14_data.m. 
%
% Note: requires internet connection to get exposure ages from online
% calculator. 
% 
% Writes 'data_exposure_ages.mat' containing exposure-age data.


disp('Calculating surface exposure ages');
init_get_surface_exposure_ages()

%% Second, run random search wrapper for all cores
%
% This saves data in 'current_results_HX.mat' data files. Note: there 
% are some settings in 'random_search_wrapper.m' and 'generate_random_et.m'
% that may affect the results.
%
% This can take a while. 

random_search_wrapper(1)
random_search_wrapper(4)
random_search_wrapper(5)

%% Third, make main figure (Fig. 3 in paper)

plot_main_figure

%% Fourth, make some of the figures in the supplementary methods
%
% Note: this doesn't make ALL the figures. A couple of them require
% changing the code in a few places to perform the experiments shown in 
% the figures (e.g., fitting to only Be-10 or C-14 data separately), so
% it is not possible to make the figures from the default results 
% generated above. 

% First plot all the random search results for the supplement, for all
% three cores. This actually takes a while because it has to load the
% result files, which are fairly large.

close all;
plot_supplement_randomsearch_results(1); drawnow;
plot_supplement_randomsearch_results(4); drawnow;
plot_supplement_randomsearch_results(5); drawnow;


%% 
close all;
plot_age_thickness_tradeoff
