% This wrapper can be used to generate model results presented in Small et al. 
% "A thinner-than-present West Antarctic Ice Sheet in the southern Weddell
% Sea Embayment during the Holocene" submitted to the Cryosphere.

% To run the wrapper for the first time requires that the follwowing files are findable
% 'EH_14_data_v2.xlsx'

% The following functions need to be on the path too
% get_core_data.m - reads in core data (see xlsx file for format)
% init_precalculate_P14_EH.m (from Balco et al., 2023) - Calculates production rates for range of elevations and mass depths
% init_make_consts.m (from Balco et al., 2023) - Makes TCN constants needed to calcualte production rates
% all_model_concs.m - predicts nuclide concentrations at depths of measured samples under synthetic ice elevation histories made by...
% random_3stages.m - this function
% get_fitting_concs.m - does misfit calculation
% fitting_model_conc.m - plots fitting model concentrations and some others (Fig 8a/c)
% plot_samples.m - Adds samples to above plot

% These functions originally written by Greg Balco are also needed
% stone2000.m - calculates scaling factors from Stone (2000)
% antatm.m - converts elevation to atmospheric pressure for Antarctica
% P_mu_total_alpha1.m - does computations for muons

% David Small
% Durham University
% September 2025

%% Start fresh
clear all, close all
addpath(genpath(pwd)) % add folders to path

%% Lets read in 14C data and run init scripts (from Balco et al., 2023)
% N.B. Below only needs to be run once. If core_data.mat, EH_consts.mat and EH_calculated_P14.mat exist then it can be commented out 

file_name = 'EH_14_data_v2.xlsx'; % define file name with core data
[D] = get_core_data(file_name); % creat struct with core data

% Set and add extra data (one location will do for all samples in our case)
D.lat = -79.944; % set latitude of core site
D.lon = -81.429; % set longtitude
D.core_cm = 50; % max length of core, as long as it covers core lengths its fine
D.max_ele = 750; % max ele for claculating PR grid (as long as core top and ice surface elevations are in this range its fine)
D.min_ele = 350; % min ele for calculating PR grid

% Run initialisation scripts (these are from Balco et al., 2023)
% init_make_consts()
% init_precalculate_P14_EH(D.lat,D.lon,D.max_ele,D.min_ele)
 
save('core_data', 'D') % Saves core data as a struct D

%% Now run random ice history model

% Model creates a 3 stage ice history (i.e. simpler than Balco et al.) 
% Stage 1 (S1) = post LGM  ice cover (200 - 50 m)
% Stage 2 = exposure/thinner than present ice
% Stage 3 re-thickening at present thickness. 
% Model calcs PRs during each stage for sample depths within core (i.e with ice cover for S1 and S3, no/thin ice during S2)
% Model then calcs final nuclide concentrations, returned in stuct (results) along with the ice thickness histories that
% produce them.

% results.N14_final is array with [n=numits] rows. 
% Each col is a depth in core (defined in results.depths) 
% Each row is an iteration with different ice history (defined in results.Ts/Zs)
% results.Ts is time vector defining 6 nodes (start/end of each stage)
% results.Zs is accompanying ice thicknesses (3 values but 6 nodes)

load core_data % load if struct already made

% define number ot iterations and core to be analysed
numits = 10000; % 10000 takes ~2 minutes on my machine.
core = 3; % 2 = BH02, 3 = BH03

z_fix = 1; % 1 = set a % of model runs to thin below sample depth. 0 = ice thickness free between 0 and present thickness.
z_pc = 10; % what % of runs to set at zero thickness

[results] = all_model_concs(D, numits, core, z_fix, z_pc); 
% Saves model results as 'model_results_BH##_date_numits_zpc_fix.mat'

% tidy up if you like
% clear all
%% So now we need to look for concentrations that "fit" what was measured

% if you have results already can just start here
load model_results_BH03_27-Jan-2026_10000_10_fix.mat %  change filename if required
load core_data
core = 3; % make sure this matches file being loaded above
z_pc = 10; % again make sure this matches

% define samples of interest for misfit calc
SOI = [1,2,4]; % samples numbered top - bottom (1 = top, etc) - paper used all
thresh = [95]; % what Rchi thresh

% option to save results (maybe dont want to if just exploring data)
to_save = [1]; % 1 = yes

% this fucntion does (very) basic misfit calculation
[results] = get_fitting_concs_RC(SOI, results, core, D, z_pc, thresh, to_save);
% [results] = get_fitting_concs(SOI, results, core, D, n, z_pc, to_save);
% Saves model results as 'fitting_model_results_BH##_date_numits_zpc_fix.mat'

% tidy up if you like
% clear all
%% Can now plot fitting results with depth (Figure 8a and c in paper)

% If you have *fitting* results can just start here
load fitting_model_results_BH03_27-Jan-2026_10000_98_fix_RChi.mat
load core_data
core = 3; % as above make sure this matches file being loaded

% this function plots result (does calcs too, inefficient i know!)
fitting_results = fitting_model_concs(D,core,results); % also plots some non fitting profiles 

if core == 2
    name = 'BH02';
elseif core == 3 
    name = 'BH03';
else 
    name = 'BH06';
end

 % save fitting models plot
 save_name = strcat(name,' fitting models_RChi');
 print(save_name,'-dpng','-r300') % export figure



