function out = generate_Cl36_params(in,consts,dFlag)

% function to calculate neutron-transport properties from xrf and 
% trace element analyses of both whole rock and a target mineral separate.
%
% No uncertainties are computed. 
%
% Syntax: out = generate_Cl36_params(in)
%
% In is data structure with various fields: either you enter XRF results
% and trace elements in ppm (option 1), or you enter weight percents of all
% elements (option 2). 
%  
% Option 1:
%   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. 
%
% Note that at present the only thing that has to be in in.tracetarget is
% the Cl concentration (I think). So that is a bit unnecessary. 
% 
% other inputs that are used are:
% in.LOI loss on ignition percent for whole rock
% in.formation_age host rock age in Ma (for nucleogenic Cl-36 calculation)
% 
% Option 2:
%
% 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)
%   or ppm (e.g., 380000).  
% Irene, note that Co is not in here, but Cr is. Not sure why. 
%
%
% Returns a structure with a whole lot of stuff in it. 
% 
% This version being renovated late 2018. 
% NOTE: nucleogenic Cl-36 calculation needs to be aware of formation age. 
% NOTE: needs to do something intelligent if only K, Ca reported for target
% chemistry. Presumably only applies to muon cross-sections...so assume
% generic feldspar? In any case, needs to be able to pass warnings. 
%

if nargin < 3
    dFlag = 0;
end

% Load consts
if nargin < 2
    load Cl36_consts;
end

out.version = 'generate_Cl36_params - 0.1-dev';

% Other definitions

N_a = 6.02214e23;
out.flags = {};

% Sort various inputs

if isfield(in,'wpe')
    % Case entering element weight fractions directly
    if isfield(in,'xrf') || isfield(in,'trace') || isfield(in,'xrftarget') || isfield(in,'tracetarget');
        error('generate_Cl36_params.m: XRF or trace element data entered redundantly with weight fxns');
    end
    if ~isfield(in,'wpetarget')
        % No separate target entered, target is same as WR
        in.wpetarget = in.wpe;
        out.flags{end+1} = 'No separate target data entered, generate_Cl36_params.m assumes target same as WR.<br>';
    end
    % Assign
    out.W_k_rock = in.wpe;
    out.W_k_target = in.wpetarget;
else
    % Case XRF and trace element vectors were entered
    % Must assemble weight fraction from these
    % unpack xrf data

    rock.sio2 = in.xrf(1); 
    rock.tio2 = in.xrf(2); 
    rock.al2o3 = in.xrf(3);
    rock.fe2o3 = in.xrf(4);
    rock.feo = in.xrf(5);
    rock.mno = in.xrf(6);
    rock.mgo = in.xrf(7);
    rock.cao = in.xrf(8);
    rock.na2o = in.xrf(9); 
    rock.k2o = in.xrf(10); 
    rock.p2o5 = in.xrf(11); 
    rock.so4 = in.xrf(12);
    rock.s = in.xrf(13); 
    rock.h2o = in.xrf(14); 
    rock.co2 = in.xrf(15);

    target.sio2 = in.xrftarget(1); 
    target.tio2 = in.xrftarget(2); 
    target.al2o3 = in.xrftarget(3);
    target.fe2o3 = in.xrftarget(4);
    target.feo = in.xrftarget(5);
    target.mno = in.xrftarget(6);
    target.mgo = in.xrftarget(7);
    target.cao = in.xrftarget(8);
    target.na2o = in.xrftarget(9); 
    target.k2o = in.xrftarget(10); 
    target.p2o5 = in.xrftarget(11); 
    target.so4 = in.xrftarget(12);
    target.s = in.xrftarget(13); 
    target.h2o = in.xrftarget(14); 
    target.co2 = in.xrftarget(15);

    % unpack trace element data

    rock.li = in.trace(1);
    rock.b = in.trace(2);
    rock.cl = in.trace(3);
    rock.cr = in.trace(4);
    rock.sm = in.trace(5);
    rock.gd = trace(6);
    rock.u = in.trace(7);
    rock.th = in.trace(8);

    target.li = in.tracetarget(1);
    target.b = in.tracetarget(2);
    target.cl = in.tracetarget(3);
    target.cr = in.tracetarget(4);
    target.sm = in.tracetarget(5);
    target.gd = in.tracetarget(6);
    target.u = in.tracetarget(7);
    target.th = in.tracetarget(8);

    % convert all to element weight fractions (g/g)

    out.W_k_rock = 0.01 *	[(0.1111 * rock.h2o); % H - 1
                        (1e-4 * rock.li); % Li - 2
                        (1e-4 * rock.b); % B - 3
                        (0.2729 * rock.co2); % C - 4
                        0; % N - 5
                        (0.88888*rock.h2o)+(0.7271*rock.co2)+(0.2581*rock.na2o)+(0.3970*rock.mgo)+(0.4708*rock.al2o3)+(0.5326*rock.sio2)+...
                    (0.5636*rock.p2o5)+(0.6662*rock.so4)+(0.1699*rock.k2o)+(0.2853*rock.cao)+(0.4006*rock.tio2)+(0.2255*rock.mno)+...
                    (0.3006*rock.fe2o3)+(0.2227*rock.feo); % O - 6
                        (0.7419*rock.na2o); % Na - 7
                        (0.6030*rock.mgo); % Mg - 8
                        (0.5292*rock.al2o3); % Al - 9
                        (0.4674*rock.sio2); % Si - 10
                        (0.4364*rock.p2o5); % P - 11
                        (0.3338*rock.so4)+(1.0*rock.s); % S - 12
                        (1e-4*rock.cl); % Cl - 13
                        0; % Ar - 14
                        (0.8301*rock.k2o); % K - 15
                        (0.7147*rock.cao); % Ca - 16
                        (0.5994*rock.tio2); % Ti - 17
                        (1e-4*rock.cr); % Cr - 18
                        (0.7745*rock.mno); % Mn - 19
                        (0.6994*rock.fe2o3)+(0.7773*rock.feo); % Fe - 20
                        (1e-4*rock.sm); % Sm - 21
                        (1e-4*rock.gd); % Gd - 22
                        (1e-4*rock.th); % Th - 23
                        (1e-4*rock.u); % U - 24
                        ];

    out.W_k_target = 0.01 *	[(0.1111 * target.h2o); 
                        (1e-4 * target.li); 
                        (1e-4 * target.b); 
                        (0.2729 * target.co2);
                        0;
                        (0.88888*target.h2o)+(0.7271*target.co2)+(0.2581*target.na2o)+(0.3970*target.mgo)+(0.4708*target.al2o3)+(0.5326*target.sio2)+...
                    (0.5636*target.p2o5)+(0.6662*target.so4)+(0.1699*target.k2o)+(0.2853*target.cao)+(0.4006*target.tio2)+(0.2255*target.mno)+...
                    (0.3006*target.fe2o3)+(0.2227*target.feo);
                        (0.7419*target.na2o);
                        (0.6030*target.mgo);
                        (0.5292*target.al2o3);
                        (0.4674*target.sio2);
                        (0.4364*target.p2o5);
                        (0.3338*target.so4)+(1.0*target.s);
                        (1e-4*target.cl);
                        0;
                        (0.8301*target.k2o);
                        (0.7147*target.cao);
                        (0.5994*target.tio2);
                        (1e-4*target.cr);
                        (0.7745*target.mno);
                        (0.6994*target.fe2o3)+(0.7773*target.feo);
                        (1e-4*target.sm);
                        (1e-4*target.gd);
                        (1e-4*target.th);
                        (1e-4*target.u);
                        ];
end
                    
% Start calculating stuff:

if size(out.W_k_rock,1) == 1
    out.W_k_rock = out.W_k_rock';
end

if size(out.W_k_target,1) == 1
    out.W_k_target = out.W_k_target';
end

% 1. general rock properties

% mole fractions f_k_rock; units of mol element/mol total atoms
% Is this needed? Apparently not. 
%out.f_k_rock =  (out.W_k_rock .* (consts.A_k.^(-1)))/(dot(out.W_k_rock,(consts.A_k.^(-1))));

% atom abundances (at/g rock) N_k_rock
% Alfimov equation 20
out.N_k_rock = N_a .* out.W_k_rock ./ consts.A_k;
out.N_k_target = N_a .* out.W_k_target ./ consts.A_k;

% mean atomic weight of rock A_rock; units are g/mol
% Alfimov equation 22
out.A_rock = dot(consts.A_k,out.N_k_rock)./sum(out.N_k_rock);

% macroscopic scattering cross-section Sigma_sc_rock
% Alfimov equation 21
out.Sigma_sc_rock = dot(consts.sigma_sc_k,out.N_k_rock); 

% 2. Ca, K, Fe, Ti weight percents; needed for spallogenic production

out.Wf_K_target = out.W_k_target(15);
out.Wf_Ca_target = out.W_k_target(16);
out.Wf_Ti_target = out.W_k_target(17);
out.Wf_Fe_target = out.W_k_target(20);	

% 3. Epithermal neutron transport properties

% Epithermal neutron diffusion coefficient D_eth_rock
% Alfimov Eqn. 27
out.D_eth_rock = 1./(3 .* out.Sigma_sc_rock .* (1 - (2./(3.*out.A_rock))));

% Effective resonance integral for absorbtion of epithermal neutrons I_eff_rock
% Alfimov Eqn. 23
out.I_eff_rock = dot(consts.I_abs_k,out.N_k_rock);
out.I_eff_target = dot(consts.I_abs_k,out.N_k_target);

% Macroscopic log decrement energy loss per neutron collision zeta_bar_rock
% Alfimov Eqn. 24
out.zeta_bar_rock = (sum(consts.zeta_k .* consts.sigma_sc_k .* out.N_k_rock))/out.Sigma_sc_rock;

% effective epithermal neutron attenuation length Lambda_eth_rock
% Alfimov Eqn. 25
out.Lambda_eth_rock = 1./(out.zeta_bar_rock .*(out.I_eff_rock + out.Sigma_sc_rock));

% Effective epithermal loss cross-section Sigma_eth_rock
% Alfimov Eqn. 25
out.Sigma_eth_rock = 1./out.Lambda_eth_rock;

% Dimensionless normalisation factor for neutron production rate R_eth_rock
% Alfimov Eqn. 30
out.R_eth_rock = (out.A_rock/consts.A_air).^0.5;

% Epithermal diffusion length L_eth_rock
% Alfimov Eqn. 28
out.L_eth_rock = (out.D_eth_rock / out.Sigma_eth_rock)^0.5;

% 4. Thermal neutron transport properties

out.D_th_rock = out.D_eth_rock; % Alfimov Eqn. 27
out.pE_th_rock = exp(-1*out.I_eff_rock/(out.zeta_bar_rock.*out.Sigma_sc_rock)); 
% Note for pE_th_rock, Alfimov appears to have corrected error in Phillips and Stone eqn. 16. 
% See Alfimov Eqn. 31 and Schimmelpfennig Eqn. 15, which appear to be the same.
out.R_th_rock = out.pE_th_rock./consts.pE_th_air; % Stone and Phillips Eqn 17
out.Sigma_th_rock = dot(consts.sigma_th_k,out.N_k_rock); % Alfimov Eqn. 26
out.Sigma_th_target = dot(consts.sigma_th_k,out.N_k_target); % same, this is a cross-section in cm2/g
out.Lambda_th_rock = 1/out.Sigma_th_rock; % same
out.L_th_rock = (out.D_th_rock./out.Sigma_th_rock)^0.5; % Alfimov eqn.29

% 5. Capture cross-sections for 35Cl --> 36Cl
% Natural abundance of 35Cl/total Cl is 0.7577 
% Total cross-sections for Cl-35 --> Cl-36 by themselves -- (cm2/atom).*(atom/g) = cm2/g
% that are used in the formula below for steady-state neutron fluxes. 
out.k_35_eth_rock = 15e-24.*(out.N_k_rock(13).*0.7577); % 15e-24 is I_a for 35Cl
out.k_35_th_rock = 43e-24.*(out.N_k_rock(13).*0.7577); % 43e-24 is sigma_th for 35Cl
out.k_35_eth_target = 15e-24.*(out.N_k_target(13).*0.7577); %15e-24 is I_a for 35Cl
out.k_35_th_target = 43e-24.*(out.N_k_target(13).*0.7577); %43e-24 is sigma_th for 35Cl

% 6. Neutron production from U,Th

% From Alfimov tables
    
% Decay neutron production
out.P_n_U = (sum(consts.S_k.*out.W_k_rock.*consts.YU_k)./sum(consts.S_k.*out.W_k_rock)).*out.W_k_rock(24).*1e6;
out.P_n_Th = (sum(consts.S_k.*out.W_k_rock.*consts.YTh_k)./sum(consts.S_k.*out.W_k_rock)).*out.W_k_rock(23).*1e6; 

% Spontaneous fission
out.P_n_sf = 0.429.*out.W_k_rock(24).*1e6; % Reflects WR

% Total neutron production
out.P_n_r = out.P_n_U + out.P_n_Th + out.P_n_sf;

% One-stage production from Alfimov

out.alfimov_f_th_target = 43e-24.*out.N_k_target(13).*0.7577./out.Sigma_th_rock; % Note params for 35Cl 
out.alfimov_P_36_r_target = out.alfimov_f_th_target .* out.P_n_r;
out.alfimov_N_36_r_target = out.alfimov_P_36_r_target./ consts.l36;

% Two-stage production from equilibrium neutron fluxes
% based on no-diffusion trivial solution of Phillips equations
% Equilibrium neutron fluxes; reflects whole rock
out.phi_eth_r = out.P_n_r .* out.Lambda_eth_rock; 
out.phi_th_r = out.Lambda_th_rock .* out.pE_th_rock .* out.P_n_r;
% By the way, I think there is an error in Schimmelpfennig (Eqn. 57) here.
% Note: this section is whay the k_th, etc. values don't have to be divided
% by the whole-rock term. 

% Calculate P and N. 
out.P_36_r_eth_target = out.phi_eth_r * out.k_35_eth_target;
out.P_36_r_th_target = out.phi_th_r * out.k_35_th_target;
out.P_36_r_target = (out.P_36_r_eth_target + out.P_36_r_th_target); % Capture x-section; reflects target composition

if in.formation_age <= 0
    % This is supposed to be the case where the formation age equals
    % the exposure age.
    % Here set nucleogenic Cl-36 to -1.
    out.N_36_r_target = -1;
elseif in.formation_age < Inf
    % Case a finite rock formation age was passed
    % Calculate finite production
    % Note formation age is in Ma
    out.N_36_r_target = out.P_36_r_target.*(1-exp(-consts.l36.*in.formation_age.*1e6))./consts.l36;
elseif isinf(in.formation_age)
    % Case infinite age, steady state form
    out.N_36_r_target = out.P_36_r_target./consts.l36; 
else
    % Not sure what ends up here
    error('generate_Cl36_parameters.m: formation age ??');
end
    
    
   
% 7. Negative muon capture yields for Ca and K
% Alfimov Eqn (10)

% Note: this part needs to not be spoofed if minimal target chemistry, e.g. only
% Ca/K/Fe/Ti, is supplied. 
 
if (sum(out.W_k_target) - (sum(out.W_k_target([15 16 17 20])))) < 0.02
    % Test is that K(15)+Ca(16)+Ti(17)+Fe(20) make up more than 98% of total
    out.flags{end+1} = 'Target appears to have only K/Ca/Ti/Fe concentrations. Negative muon capture yields likely incorrect.<br>';
    % If this is the case, we have to assume the rest of the rock is made
    % up of something that has an atomic mass and a W-value. For atomic
    % mass, use "standard rock" which has mass 22. For the W-value,
    % observe that W is more or less correlated with A and use W = 1.
    % Note that it is a terrible idea to have both W_k in this code (mass
    % fraction) and consts.W_k (negative muon yields) have the same
    % variable name. 
    missing_W_k = 1 - sum(out.W_k_target([15 16 17 20]));
    missing_N_k = N_a .* missing_W_k ./ 22; % standard rock
    out.fC_Ca_target = consts.W_k(16).*out.N_k_target(16)./(sum(out.N_k_target([15 16 17 20]).*consts.W_k([15 16 17 20])) + 1.*missing_N_k);
    out.fC_K_target = consts.W_k(15).*out.N_k_target(15)./(sum(out.N_k_target([15 16 17 20]).*consts.W_k([15 16 17 20])) + 1.*missing_N_k);   
else
    % Case appears to be complete chemistry, use usual formula
    out.fC_Ca_target = consts.W_k(16).*out.N_k_target(16)./(sum(out.N_k_target.*consts.W_k));
    out.fC_K_target = consts.W_k(15).*out.N_k_target(15)./(sum(out.N_k_target.*consts.W_k));
end



% Also total Cl-36 yield from negative muon capture on Ca, K
% This is Y^{36}_{Ca}, etc. in Alfimov but k_neg in P_mu_total.m and
% related code

% 0.969 is isotopic abundance of Ca-40 and 0.9326 is same of K-39
out.Y_Ca_target = out.fC_Ca_target.*0.969.*consts.f_surv_k(16).*consts.f_star_k(16);
out.Y_K_target = out.fC_K_target.*0.9326.*consts.f_surv_k(15).*consts.f_star_k(15);

% 8. Neutron yield from negative muon capture

out.Y_n_mu = sum(consts.y_n_k.*consts.f_surv_k.*(out.N_k_rock.*consts.W_k./sum(out.N_k_rock.*consts.W_k)));


		
