function init_precalculate_P10_P14()

% This script precalculates production rates for a core site. Once
% production rates at a range of depths have been calculated and stored in
% a data file, you can just look them up when you need them in a forward
% model calculation instead of having to do the whole calculation.
%
% Note: despite the name, this also calculates Al-26 production rates.
% These aren't used in fitting but are used in comparing to a few Al-26
% data in core H1. 
%
% 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. 

% First generate constants, whether needed or not
init_make_consts_KP()

%clear all;

load KP_consts;

% Set location. One location is adequate for all cores and ridge surface
% samples. 

rlat = -75.21534;
rlon = -110.96041;

% Set up calculation mesh

P.v_surfelv = linspace(0,1000,51); % Rows are surface elv
P.v_zg = [0.1 logspace(0,5,99)]; % Cols are depth below surface

% Note: also transform to log to make log-interpolation faster. 

P.log_v_surfelv = log(P.v_surfelv);
P.log_v_zg = log(P.v_zg);

% Preallocate space for production rate grids

P.P14 = zeros(length(P.v_surfelv),length(P.v_zg));
P.P10 = P.P14; 
P.P26 = P.P10;

% Make mesh of calculated production rates

for elv = 1:length(P.v_surfelv)
    currentp = antatm(P.v_surfelv(elv));
    current_spsf = stone2000(rlat,currentp,1);
    
    for depth = 1:length(P.v_zg)
        
        % Production rate for C-14
        current_P14sp = consts.P14_ref_St.*current_spsf.*exp(-P.v_zg(depth)./consts.Lsp);
        current_P14mu = P_mu_total_alpha1(P.v_zg(depth),currentp,consts.mc14);
        P.P14(elv,depth) = current_P14sp + current_P14mu;
        
        % Production rate for Be-10
        current_P10sp = consts.P10_ref_St.*current_spsf.*exp(-P.v_zg(depth)./consts.Lsp);
        current_P10mu = P_mu_total_alpha1(P.v_zg(depth),currentp,consts.mc10);
        P.P10(elv,depth) = current_P10sp + current_P10mu;
        
        % Production rate for Al-26
        current_P26sp = consts.P26_ref_St.*current_spsf.*exp(-P.v_zg(depth)./consts.Lsp);
        current_P26mu = P_mu_total_alpha1(P.v_zg(depth),currentp,consts.mc26);
        P.P26(elv,depth) = current_P26sp + current_P26mu;
        
    end
    disp(int2str(P.v_surfelv(elv)));
end

P.date = datestr(now);

% Save that as a data file. 
save KP_calculated_P P


