function P = calculate_production_profile_web(in,consts,plotFlag)
%
% To fix:
%   
%   What happens to LOI? Assume all LOI is H2O? Or do you enter that
%   separately? 
% 
% This code generates a Cl-36 production profile, most parts of which are
% normalized for later scaling adjustments. This version is for online
% use and relies on lots of precalculation. 
%
% 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)
%   in.LOI percent loss-on-ignition
%       (optional, defaults to zero)
%   in.formation_age formation age of host rock (Ma)
%       (optional, defaults to infinity)
%
% Returns a data structure with:
%   out.zz z coordinate (g/cm2)
%   out.Pxx Production rate from lots of different pathways - lots of 'xx'
%       Note some of these are normalized, scaling applied downstream
%   out.params Another data structure with calculated parameters
% 
% Greg Balco
% Berkeley Geochronology Center
% Feb 2018
% Not licensed for any reuse or distribution. 
%

if nargin < 3
    plotFlag = 0;
end

% load constants

if nargin < 2
    load Cl36_consts;
end

% Get parameters

d = generate_Cl36_params(in,consts);

P.version_profile = 'calculate_production_profile_web - 0.1-dev';
P.version_params = d.version;

% Initialize P

P.flags = d.flags;

% 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,
% Fe, Ti
% Note this is normalized to surface = 1; scaling applied downstream

P.P36sp_Ca = exp(-zz./Lsp).*d.W_k_target(16);
P.P36sp_K = exp(-zz./Lsp).*d.W_k_target(15);
P.P36sp_Ti = exp(-zz./Lsp).*d.W_k_target(17);
P.P36sp_Fe = exp(-zz./Lsp).*d.W_k_target(20);


% 2. Neutron capture
% Follows Alfimov
% Again, this is normalized, scaling and P_f_0 are set downstream

% 2.1. Epithermal
k_epi_rock = d.R_eth_rock./(d.Sigma_eth_rock - (d.D_eth_rock./(Lsp.^2)));
k_epi_air = 1./(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;

% Sum neutrons
P.P36n_Cl = P.P36sp_eth_Cl + P.P36sp_th_Cl;

% 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 = consts.sigma0_Ca;
mcK.sigma0 = consts.sigma0_K;

% Muon surface flux approximations
P.P_neg_0 = consts.mu.P0_neg.*exp((in.pressure-1013.25)./consts.mu.Leff_neg); 
P.L_neg = 1./(consts.mu.a_neg + in.pressure.*consts.mu.b_neg);
P.P_fast_0 = consts.mu.P0_fast.*exp((in.pressure-1013.25)./consts.mu.Leff_fast); 
P.L_fast = 1./(consts.mu.a_fast + in.pressure.*consts.mu.b_fast);

% Turn into depth profile
P.P36_neg_Ca = P.P_neg_0.*mcCa.k_neg.*exp(-zz./P.L_neg);
P.P36_fast_Ca = P.P_fast_0.*mcCa.Natoms.*mcCa.sigma0.*exp(-zz./P.L_fast);
P.P36_neg_K = P.P_neg_0.*mcK.k_neg.*exp(-zz./P.L_neg);
P.P36_fast_K = P.P_fast_0.*mcK.Natoms.*mcK.sigma0.*exp(-zz./P.L_fast);
P.Pn_neg = P.P_neg_0.*d.Y_n_mu.*exp(-zz./P.L_neg); % This should be neutrons/cm2/yr
P.Pn_fast = P.P_fast_0.*4.8e-6.*exp(-zz./P.L_fast); % This should be neutrons/cm2/yr


% The below is previous code that calls P_mu_total_alpha1. Also plotting
% to compare results.
% 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');
%     P2.P36_neg_Ca(a) = mCa.P_neg;
%     P2.P36_neg_K(a) = mK.P_neg;
%     P2.P36_fast_Ca(a) = mCa.P_fast;
%     P2.P36_fast_K(a) = mK.P_fast;
%     P2.Pn_neg(a) = mCa.R.*d.Y_n_mu; % This should be neutrons/cm2/yr
%     P2.Pn_fast(a) = mCa.phi.*mCa.Ebar.*4.8e-6; % Also should be neutrons/cm2/yr
% end
%     
% if 1
%     figure; 
%     subplot(3,1,1);
%     plot(zz,P.P36_fast_Ca,'r',zz,P2.P36_fast_Ca,'r:',zz,P.P36_neg_Ca,'b',zz,P2.P36_neg_Ca,'b:'); title('Ca: red, fast, blue, negative');
%     set(gca,'xlim',[0 20]);
%     subplot(3,1,2);
%     plot(zz,P.P36_fast_K,'r',zz,P2.P36_fast_K,'r:',zz,P.P36_neg_K,'b',zz,P2.P36_neg_K,'b:'); title('K: red, fast, blue, negative');
%     set(gca,'xlim',[0 20]);
%     subplot(3,1,3);
%     plot(zz,P.Pn_fast,'r',zz,P2.Pn_fast,'r:',zz,P.Pn_neg,'b',zz,P2.Pn_neg,'b:'); title('n: red, fast, blue, negative');
%     set(gca,'xlim',[0 20]);
% 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;

% Now compute total muon-derived production, this is
% scaling-method-independent and in real units

P.P36mu = P.P36_neg_Ca + P.P36_fast_Ca + P.P36_neg_K + P.P36_fast_K + P.P36mu_eth_Cl + P.P36mu_th_Cl;

% 4. Nucleogenic production, if rock age = formation age
P.P36nuc = 0; % Normal case, not counted as production during exposure, zero, note not vector
if d.N_36_r_target < 0 % flagged in generate_Cl36_params as -1 atom
    % This case, nucleogenic production gets added to production
    P.P36nuc = P.P36nuc + d.P_36_r_target;
end

% 5. Finally, calculate entire production profile
% Note this is meaningless when some pathways are normalized to unit
% scaling...probably should be removed....
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 + P.P36nuc;

% Assign params structure to output

P.params = d;

% 5. Plotting

if plotFlag == 1
    if 0
        figure;
        plot(P.P36sp_Ca,zz,'r');hold on;
        plot(P.P36sp_K,zz,'r');
        set(gca,'ydir','reverse');
        set(gca,'ylim',[0 500]);

        grid on
        plot(P.P36sp_eth_Cl,zz,'g');
        plot(P.P36sp_th_Cl,zz,'g');

        plot(P.P36_neg_Ca,zz,'b');
        plot(P.P36_neg_K,zz,'b');

        plot(P.P36_fast_Ca,zz,'c');
        plot(P.P36_fast_K,zz,'c');

        plot(P.P36mu_eth_Cl,zz,'m');
        plot(P.P36mu_th_Cl,zz,'m');

        plot(P.P36,zz,'k');
    end
    if 1
        figure;
        plot(P.P36sp_Ca+P.P36sp_K,zz,'r.'); hold on;
        plot(P.P36sp_eth_Cl+P.P36sp_th_Cl,zz,'g.');
        plot(P.P36_neg_Ca+P.P36_neg_K,zz,'b.');
        plot(P.P36_fast_Ca+P.P36_fast_K,zz,'c.');
        plot(P.P36mu_eth_Cl+P.P36mu_th_Cl,zz,'m.');
        % plot(P.P36,zz,'k',P.P36,zz,'k.');  
        set(gca,'ydir','reverse');
        set(gca,'ylim',[0 600]);
    end
end
    
% Return

