function P = PofZ(zin,m)

% Returns production rate at any depth z (g/cm2) at SLHL.
% m is structure with precalculated data:
%   m.zz = depth vector
%   m.P = corresponding production rate vector
%   m.Lexp = attenuation length for extending beyond max zz
% output is production rate
% vector or matrix inputs for zin are OK
%
% Exponential approximation for muons below 20,000 g/cm2
% 
% The idea is to allow for reasonably fast time integration of muon
% production during erosion.
%
% 
% Greg Balco -- Berkeley Geochronology Center -- May 2016
% Revised to be more general Feb 2018
 
% Unwrap matrix input
z = reshape(zin,1,numel(zin));

inbounds = (z <= max(max(m.zz)));
outofbounds = ~inbounds;

P(inbounds) = interp1(m.zz,m.P,z(inbounds));
P(outofbounds) = m.P(end).*exp(-(z(outofbounds)-max(m.zz))./m.Lexp);

% Rewrap matrix input
P = reshape(P,size(zin,1),size(zin,2));
