function out = KPP(elvs,zs,nuclide,P)

% This returns the production rate for arbitrary surface elevation and mass
% depth below the surface at the Kay Peak core site (or the adjacent
% ridge). 
%
% Syntax:  out = KPP(elvs,zs,nuclide,P)
% 
% Input args: elvs is scalar or vector surface elevation (m)
% zs is scalar or vector mass depth below the surface (g/cm2)
% Note: elvs and zs must have same dimensions 1 x n or n x 1 
% nuclide is 10, 14, or 26; if omitted, calculates all and returns data
% structure instead of single vector
% P is data structure containing production rate information
% Note: if not supplied, loads data file
%
% Output is vector of production rates for desired nuclide of same size as
% elvs and zs, unless 'nuclide' is not passed, in which case it's a data
% structure with out.P10, out.P14, etc. 
%
% 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. 


if nargin < 4
    load KP_calculated_P
end

if nargin < 3
    nuclide = -1;
end

% Remove zeros by setting to 0.1
zs(find(zs == 0)) = zs(find(zs == 0)) + 0.1;

% 2-d interpolation. Note: log interpolation for depth below the surface
% and linear interpolation for surface elevation. You could use log for the
% surface elevation as well, but that is not changing as much so it is is
% less relevant. 

% Also note: log interpolation vector (P.log_v_zg) has already been
% calculated upstream in the preprocessor. 

if nuclide == 10
    out = interp2(P.log_v_zg,P.v_surfelv,P.P10,log(zs),elvs);
elseif nuclide == 14
    out = interp2(P.log_v_zg,P.v_surfelv,P.P14,log(zs),elvs);
elseif nuclide == 26
    out = interp2(P.log_v_zg,P.v_surfelv,P.P26,log(zs),elvs);
else
    out.P10 = interp2(P.log_v_zg,P.v_surfelv,P.P10,log(zs),elvs);
    out.P14 = interp2(P.log_v_zg,P.v_surfelv,P.P14,log(zs),elvs);
    out.P26 = interp2(P.log_v_zg,P.v_surfelv,P.P26,log(zs),elvs);
end


        
