function P = calculate_production_profile(in,consts,plotFlag)

% This code generates a Cl-36 production profile. 
% Uses "St" scaling. 
%
% In is data structure with various fields
%   in.xrf and in.xrftarget contain XRF analyses of WR and target, 
%       respectively. They are vectors like:
% [sio2 tio2 al2o3 fe2o3 feo mno mgo cao na2o k2o p2o5 so4 s h2o co2]
%       units are weight percent oxide
%       Order of oxides is nonnegotiable
%   in.trace and in.tracetarget contain ppm concentrations of trace
%   elements for WR and target, respectively. Elemental order is
% [li b cl cr sm gd u th]
%       units are ppm wt. Order is nonnegotiable. 
% 
% If you have weight percent elements instead, you can enter
%   a single vector containing fractions of major and trace elements 
%   together as in.wpe and in.wpetarget.
%   For this, use official 24 elements in atomic number order:
%   [H Li B C N O Na Mg Al Si P S Cl Ar K Ca Ti Cr Mn Fe Sm Gd Th U]
%   Note that this is fractions (e.g., 0.38) not percent (e.g., 38).  
%
% Also:
%   in.zbot max depth (g/cm2) at which to compute production
%       (optional, defaults to 20000)
%   in.lat, in.lon lat/long of site
%       (optional, defaults to 72.5796 N, 38.4592 W (Greenland Summit)
%   in.pressure site pressure (hPa)
%       (optional, defaults to 1013.25)
%   in.Lsp spallogenic effective attenuation
%       (optional, defaults to 140 g/cm2)
%
% Returns a data structure with:
%   out.zz z coordinate (g/cm2)
%   out.Pxx Production rate from lots of different pathways - lots of 'xx'
%   out.P36 Total production rate in atoms/g/yr
%   out.params Another data structure with calculated parameters
% 
% Part of calculations in Balco et al.,
% 'Chlorine-36/beryllium-10 burial dating of alluvial fan sediments
% associated with the Mission Creek strand of the San Andreas fault system,
% California, USA
%
% Copyright Greg Balco
% Berkeley Geochronology Center
% Feb 2018 - March 2019
% Provided for review of the accompanying paper. 
% Not licensed for other reuse or redistribution. 
% No guarantee of accuracy or fitness for any purpose. 


if nargin < 3
    plotFlag = 0;
end

% load constants

if nargin < 2
    load Cl36_consts;
end

% Get parameters

d = generate_Cl36_params(in,consts);

% Define depth
if ~isfield(in,'zbot')
    in.zbot = 20000;
end
zz = [0 logspace(0,log10(in.zbot),200)];
P.zz = zz;

% Other parameters

if ~isfield(in,'lat') || ~isfield(in,'lon')
    in.lon = -38.4592; in.lat = 72.5796;
end

if ~isfield(in,'pressure')
    in.pressure = 1013.25;
end

% 1. Spallogenic

if ~isfield(in,'Lsp')
    Lsp = 140;
else
    Lsp = in.Lsp;
end

% Obtain spallogenic scaling factor and spallogenic production on Ca, K
% If in.normalized is 1, then this is just considered 1 and send upstream
if isfield(in,'normalized_sp')
    if in.normalized_sp == 1
        SFSP = 1;
        P.P36sp_Ca = exp(-zz./Lsp).*d.W_k_target(16);
        P.P36sp_K = exp(-zz./Lsp).*d.W_k_target(15);
    end
else
    % Not set, actually calculate using St
    SFSP = stone2000(in.lat,in.pressure,1);
    P.P36sp_Ca = SFSP.*consts.Psp_0_Ca.*exp(-zz./Lsp).*d.W_k_target(16);
    P.P36sp_K = SFSP.*consts.Psp_0_K.*exp(-zz./Lsp).*d.W_k_target(15);
end

% Spallogenic production on Ca and K
% This needs to be fixed for proper scaling. 


% 2. Neutron capture
% Follows Alfimov

% 2.1. Epithermal
k_epi_rock = SFSP.*consts.P_f_0.*d.R_eth_rock./(d.Sigma_eth_rock - (d.D_eth_rock./(Lsp.^2)));
k_epi_air = SFSP.*consts.P_f_0./(consts.Sigma_eth_air - (consts.D_eth_air./(Lsp.^2)));

A = ((consts.D_eth_air.*k_epi_air) - (d.D_eth_rock.*k_epi_rock))./Lsp;
B = 1./((consts.D_eth_air./consts.L_eth_air) + (d.D_eth_rock./d.L_eth_rock));

k_esc_epi_rock = B.*(A + ( (k_epi_air - k_epi_rock).*consts.D_eth_air./consts.L_eth_air));
k_esc_epi_air = B.*(A + ( (k_epi_rock - k_epi_air).*d.D_eth_rock./d.L_eth_rock));

phi_eth = k_epi_rock.*(exp(-zz./Lsp)) + k_esc_epi_rock.*exp(-zz./d.L_eth_rock);

P.P36sp_eth_Cl = phi_eth.*d.k_35_eth_target;

% 2.2.Thermal

k_th_rock = (d.pE_th_rock./d.Lambda_eth_rock).*k_epi_rock./(d.Sigma_th_rock - (d.D_th_rock./(Lsp.^2)));
k_corr_rock =    (d.pE_th_rock./d.Lambda_eth_rock).*k_esc_epi_rock./(d.Sigma_th_rock - (d.D_th_rock./((d.L_eth_rock).^2)));
k_th_air =  (consts.pE_th_air./consts.Lambda_eth_air).*k_epi_air./(consts.Sigma_th_air - (consts.D_th_air)./(Lsp.^2));
k_corr_air = (consts.pE_th_air./consts.Lambda_eth_air).*k_esc_epi_air./(consts.Sigma_th_air - (consts.D_th_air)./((consts.L_eth_air).^2));

E = consts.D_th_air.*((k_th_air./Lsp) - (k_corr_air./consts.L_eth_air)) - d.D_th_rock.*((k_th_rock./Lsp) + (k_corr_rock./d.L_eth_rock));
F = k_th_rock - k_th_air + k_corr_rock - k_corr_air;
G = 1./((consts.D_th_air./consts.L_th_air)+(d.D_th_rock./d.L_th_rock));

k_esc_th_rock = (E - F.*(consts.D_th_air./consts.L_th_air)).*G;
k_esc_th_air = (E + F.*(d.D_th_rock./d.L_th_rock)).*G;

phi_th = k_th_rock.*exp(-zz./Lsp) + k_esc_th_rock.*exp(-zz./d.L_th_rock) + k_corr_rock.*exp(-zz./d.L_eth_rock);

P.P36sp_th_Cl = phi_th.*d.k_35_th_target;

% 3. Muon stuff

% Allocate
P.P36_neg_Ca = zeros(size(zz));
P.P36_neg_K = zeros(size(zz));
P.P36_fast_Ca = zeros(size(zz));
P.P36_fast_K = zeros(size(zz));
P.Pn_neg = zeros(size(zz));
P.Pn_fast = zeros(size(zz));

% Muon interaction cross-sections
mcCa.Natoms = d.N_k_target(16);
mcK.Natoms = d.N_k_target(15);
mcCa.k_neg = d.Y_Ca_target; 
mcK.k_neg = d.Y_K_target;
mcCa.sigma0 = 7.3e-30; % Heisinger
mcK.sigma0 = 9.4e-30; % Marrero

% Precalculation is needed here. I guess you
% have to precalculate fast separately, and then precalculate capture 
% except for the Y-value, and then multiply that on at the end. 
for a = 1:length(zz)
    mCa = P_mu_total_alpha1(zz(a),in.pressure,mcCa,'yes');
    mK = P_mu_total_alpha1(zz(a),in.pressure,mcK,'yes');
    P.P36_neg_Ca(a) = mCa.P_neg;
    P.P36_neg_K(a) = mK.P_neg;
    P.P36_fast_Ca(a) = mCa.P_fast;
    P.P36_fast_K(a) = mK.P_fast;
    P.Pn_neg(a) = mCa.R.*d.Y_n_mu; % This should be neutrons/cm2/yr
    P.Pn_fast(a) = mCa.phi.*mCa.Ebar.*4.8e-6; % Also should be neutrons/cm2/yr
end
    
P.Pn_mu = P.Pn_neg + P.Pn_fast; % Units of P.Pn_fast correct? 

% Now compute neutron capture production from muon-produced neutrons
% Ignoring near-surface diffusion

phi_eth_mu = P.Pn_mu .* d.Lambda_eth_rock; 
phi_th_mu = d.Lambda_th_rock .* d.pE_th_rock .* P.Pn_mu;

P.P36mu_eth_Cl = phi_eth_mu.* d.k_35_eth_target;
P.P36mu_th_Cl = phi_th_mu.*d.k_35_th_target;

% 4. Finally, calculate entire production profile

P.P36 = P.P36sp_Ca + P.P36sp_K + P.P36sp_eth_Cl + P.P36sp_th_Cl + P.P36_neg_Ca + P.P36_fast_Ca + P.P36_neg_K + P.P36_fast_K + P.P36mu_eth_Cl + P.P36mu_th_Cl;

% Assign params structure to output

P.params = d;

% 5. Plotting

if plotFlag == 1
    Cacolor = [0.5 0.5 0.5];
    Kcolor = [0 0 1];
    Clcolor = [0.3 0.7 0.3];
    
    % Generate figure like Fig. 5 in paper
    figure('pos',[440   376   800   400]);
    
    % Log axes on left side
    subplot(1,2,1);
    % Spallation is solid line, neutrons dashed, muons dotted
    % Calcium is blue
    plot(P.P36sp_Ca,zz,'-','color',Cacolor,'linewidth',1.5); hold on;
    plot(P.P36_neg_Ca+P.P36_fast_Ca,zz,'o','color',Cacolor','markersize',2,'color',Cacolor);
    plot(P.P36_neg_Ca+P.P36_fast_Ca,zz,'-','color',Cacolor);
    % Potassium is lighter blue
    plot(P.P36sp_K,zz,'-','color',Kcolor,'linewidth',1.5); hold on;
    plot(P.P36_neg_K+P.P36_fast_K,zz,'o','markerfacecolor',Kcolor','markersize',2,'color',Kcolor);
    plot(P.P36_neg_K+P.P36_fast_K,zz,'-','color',Kcolor);
    % Chlorine is green
    plot(P.P36sp_eth_Cl+P.P36sp_th_Cl+P.P36mu_eth_Cl+P.P36mu_th_Cl,zz,'-','color',Clcolor,'linewidth',1.5);
    %plot(P.P36mu_eth_Cl+P.P36mu_th_Cl,zz,'o','color',Clcolor,'markersize',2,'color',Clcolor);
    %plot(P.P36mu_eth_Cl+P.P36mu_th_Cl,zz,'-','color',Clcolor);
    % Total is black
    plot(P.P36,zz,'k','linewidth',2.5);
    
    set(gca,'ydir','reverse','ylim',[0 500],'xscale','log');
    set(gca,'fontsize',14);
    xlabel('^{36}Cl production rate (atoms g^{-1} yr^{-1})');
    ylabel('Mass depth (g cm^{-2})');
    grid on
    
    % Linear at right
    
    % Linear axes on left side
    subplot(1,2,2);
    % Spallation is solid line, neutrons dashed, muons dotted
    % Calcium is blue
    plot(P.P36sp_Ca,zz,'-','color',Cacolor,'linewidth',1.5); hold on;
    plot(P.P36_neg_Ca+P.P36_fast_Ca,zz,'o','color',Cacolor','markersize',2,'color',Cacolor);
    plot(P.P36_neg_Ca+P.P36_fast_Ca,zz,'-','color',Cacolor);
    % Potassium is lighter blue
    plot(P.P36sp_K,zz,'-','color',Kcolor,'linewidth',1.5); hold on;
    plot(P.P36_neg_K+P.P36_fast_K,zz,'o','markerfacecolor',Kcolor','markersize',2,'color',Kcolor);
    plot(P.P36_neg_K+P.P36_fast_K,zz,'-','color',Kcolor);
    % Chlorine is green
    plot(P.P36sp_eth_Cl+P.P36sp_th_Cl+P.P36mu_eth_Cl+P.P36mu_th_Cl,zz,'-','color',Clcolor,'linewidth',1.5);
    %plot(P.P36mu_eth_Cl+P.P36mu_th_Cl,zz,'o','color',Clcolor,'markersize',2,'color',Clcolor);
    %plot(P.P36mu_eth_Cl+P.P36mu_th_Cl,zz,'-','color',Clcolor);
    % Total is black
    plot(P.P36,zz,'k','linewidth',2.5);
    
    set(gca,'ydir','reverse','ylim',[0 500]);
    set(gca,'fontsize',14);
    xlabel('^{36}Cl production rate (atoms g^{-1} yr^{-1})');
    grid on
   
    
    
end
    
% Return

