function init_precalculate_P14_EH(lat,lon,max_ele,min_ele)

% 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. 
% N.B This functionality is currently commented out (D. Small). Above
% refers to data in Balco et al. 2023.
%
% Greg Balco
% Berkeley Geochronology Center
% Feb, 2022
% 
% Originally written to accompany a paper entitled 'REVERSIBLE ICE SHEET THINNING IN THE
% AMUNDSEN SEA EMBAYMENT DURING THE LATE HOLOCENE,' by Balco, Brown,
% Nichols, Venturelli, and others. 
%

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

% clear all;

load EH_consts;

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

rlat = lat;
rlon = lon;

% Set up calculation mesh

P.v_surfelv = linspace(min_ele,max_ele,51); % Rows are surface elv in metres
P.v_zg = [0.1 logspace(0,5,99)]; % Cols are mass depth below surface in g cm

% 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));


% 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 EH_calculated_P14


