function out = init_make_consts()

% This function creates and saves a structure with relevant constants and
% external data for exposure age calculations. 
%
% 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. 
%


consts.version = datestr(now,'yyyy-mm-dd');
consts.prepdate = fix(clock);

% Density of ice

consts.rhoice = 0.917;

% Decay constants

% Be-10
consts.l10 = -log(0.5)./1.387e6; % Chmeleff/Korschinek value
dldt = -log(0.5).*(1.387e6^-2);
consts.dell10 = sqrt((dldt.*0.012e6)^2); 

% Al-26 -- value compatible with Nishiizumi standards
% lambda = 9.83e-7 --> t(1/2) = 0.705e6 yr
% See Nishiizumi (2004) for details.
consts.l26 = 9.83e-7;
consts.dell26 = 2.5e-8;

% C-14 decay constant
consts.l14 = -log(0.5)./5730; 
consts.dell14 = consts.l14.*0.005; % ????

% Production rates for use with Lal-Stone scheme. 

consts.P10_ref_St = 4.086; % From recalibrating with CRONUS primary 20161204
consts.delP10_ref_St = 4.086.*0.079; % bootstrap from secondary cal data set

consts.P26_ref_St = 28.535; % From recalibrating with CRONUS primary 20161204
consts.delP26_ref_St = 28.535.*0.104;

consts.P14_ref_St = 12.1; % From Tulane CRONUS-A data 
consts.delP14_ref_St = 12.1 * 0.05; % ?? 

% LSDn -- these are nondimensional correction factors, not dimensional
% production rates. At the moment I don't think these are used. 

% Be-10 is calibrated
consts.P10_ref_LSDn = 0.849; % 20161204 Be-10 calibration from CRONUS primary
consts.delP10_ref_LSDn = 0.849.*0.059;

% C-14 calibrated from Tulane CRONUS-A data 
consts.P14_ref_LSDn = 0.725;
consts.delP14_ref_LSDn = 0.725.*0.05; % ???  

% Al-26 is calibrated
consts.P26_ref_LSD = 0.819; % 20161204 Al-26 calibration from CRONUS primary
consts.delP26_ref_LSDn = 0.819 .* 0.086; % Bootstrap

% Muon interaction parameters

% Constants. See paper in QG (2017).  

% Be-10, model 1A
consts.mc10.Natoms = 2.006e22;
consts.mc10.k_neg = 0.00191 .* 0.704 .* 0.1828; % From BCO fit
consts.mc10.sigma0 = 0.280e-30; % From BCO fit

% Al-26, model 1A
consts.mc26.Natoms = 1.003e22;
consts.mc26.k_neg = 0.0133 .* 0.296 .* 0.6559; % From BCO fit
consts.mc26.sigma0 = 3.89e-30; % From BCO fit

% C-14, model 1A
consts.mc14.Natoms = 2.006e22;
consts.mc14.k_neg = 0.116 .* 0.704 .* 0.1828; % From Leymon High fit
consts.mc14.sigma0 = 0.45e-27./190; % From Heisinger, alpha = 1

% Lsp

consts.Lsp = 140;


save EH_consts consts

disp(['Constants version ' consts.version]);
disp('Saved'); 


