function ages = get_Cl36_ages_v3(in,control,consts)

% This is the code that is supposed to calculate Cl-36 ages.
% in is the input data structure, which is the output of validate_..
%
% Development, not licensed for use or distribution
% Greg Balco -- Berkeley Geochronology Center -- December 2019
% Various updates and bug fixes to Sept 2021

% Punch list:
%   Must generate error flag if formation age < exposure age. 
%   NOTE: need to do something about LOI/water content..in.LOI
%   Related: renormalize compositions to 100%?
%   NOTE: nucleogenic Cl-36 needs to be aware of formation age...that needs
%       to get passed to and incorporated in generate_Cl36_params.m...is
%       this done already?
%   Must classify as Ca,K,etc. dominated for 2-nuclide plots?
%   Should accept XRF data as text input lines...this goes in validate_
%   Should do something about muon likelihoods when only target element
%       concentrations are known for target fraction. Already generates
%       error...should assume standard rock instead? 
%   Calculate production rate uncerts for other scaling besides St
%   Actually calculate ext uncert for Nnorm when saturated


if nargin < 3
    load Cl36_consts;
end

% Defaults
if nargin < 2
    % I don't think the 'control' input structure does anything at the
    % moment.
    control.cal = 0; % Not a calibration data set, calculate ages
end

% Version
in.version_ages = 'get_Cl36_ages_v3 - 0.1-dev';


% ---------------- COMPUTE SAMPLE-SPECIFIC PARAMETERS ------------------

% 1. Compute sample-specific parameters: pressure, various SFs

% Calculate pressures with minimal number of calls to ERA40atm.m.
% In this section is where we would pull an elevation history from an 
% isostatic rebound model. 

calc_ps_ERA40 = find(strcmp(in.s.aa,'std')); % Find indices to calculate
calc_ps_ant = find(strcmp(in.s.aa,'ant'));

% Main point here is to call ERA40.m only once
in.s.pressure(calc_ps_ERA40) = ERA40atm(in.s.lat(calc_ps_ERA40),in.s.long(calc_ps_ERA40),in.s.elv(calc_ps_ERA40));
in.s.pressure(calc_ps_ant) = antatm(in.s.elv(calc_ps_ant));

% Note if atmospheric pressure were made time-dependent, then one would 
% need a cell array rather than a vector, which will incur lots of down-
% stream code changes. 

% ------------ COMPUTE NUCLIDE-SPECIFIC PARAMETERS  --------------------

% Get L; this is placeholder for now but may be spatially variable, so it 
% is a vector. It belongs to n not s so that it can also vary by nuclide if
% needed, although this isn't relevant for Cl-36 code. 

in.n.Lsp = zeros(size(in.n.index)) + Lsp();

% Thickness correction; by extension this is also a property of n
in.n.sf_thick = (thickness(in.s.thick(in.n.index),in.n.Lsp,in.s.rho(in.n.index)))';    
% Transpose is because thickness.m returns wrong shape...??

% St SF. Then on the other hand this one is technically a nuclide-
% independent property of the sample, but we put it in the nuclide 
% properties for convenience. 

in.n.sf_St = (stone2000(in.s.lat(in.n.index),in.s.pressure(in.n.index),1))'; 
% Note transpose because stone2000.m returns wrong shape.

% Switch for pathological situations
in.n.ok_to_calculate = ones(size(in.n.N));

% ------------ COMPUTE PRODUCTION PROFILES -----------------------------

% Inputs are in.wpe, in.wpetarget, in.normalized_sp, in.lat, in.long,
% in.pressure, in.zbot
%


for a = 1:length(in.n.index)
    in.flags{a} = [];
    % Note this is stepping through n, so if you want a parameter from
    % structure s, you need to get s.whatever(in.n.index(a)), not
    % s.whatever(a). 
    this_input.zbot = 1000; 
    % Assign element concentrations
    this_input.wpe = in.n.wpe{a};
    if isempty(find(in.n.wpetarget{a} > 0)) 
        % Case no data for target, assume same as WR and flag
        in.flags{a} = [in.flags{a} 'Sample ' in.s.sample_name{in.n.index(a)} ': get_Cl36_ages.m assumes WR and target chemistry are the same.<br>'];
        this_input.wpetarget = in.n.wpe{a};
    else
        this_input.wpetarget = in.n.wpetarget{a};
    end
    % Assign lat, long, pressure
    this_input.lat = in.s.lat(in.n.index(a)); this_input.lon = in.s.long(in.n.index(a));
    this_input.pressure = in.s.pressure(in.n.index(a));
    % Assign target Cl concentration
    this_input.wpetarget(13) = in.n.Cl(a).*1e-6;
    % Also put that back in the data structure...have to do it at the end
    % to deal with WR=target case
    in.n.wpetarget{a}(13) = in.n.Cl(a).*1e-6;
    in.n.delwpetarget{a}(13) = in.n.delCl(a).*1e-6;
    % Assign whole rock LOI.
    % This needs to be fixed to split LOI for target/WR. Oops.
    % Also, this needs to do something. What does it do? 
    this_input.LOI = in.n.LOI(a);
    % Assign host rock formation age
    this_input.formation_age = in.s.rock_age_Ma(in.n.index(a));
    % Now calculate production profile
    in.profile{a} = calculate_production_profile_web(this_input,consts,0);
    % Initial calculation of surface production fractions using St scaling.
    % This is also used to estimate production rate uncertainty for error
    % propagation later. 
    
    % These include uncertainties in calibrated pathways as well as in
    % target element concentrations. 
    % Note: extra effort to make sure this doesn't produce NaNs or Infs if there is zero
    % of the element in question. 
    relErrConcCa = in.n.delwpetarget{a}(16)./this_input.wpetarget(16); if ~isfinite(relErrConcCa); relErrConcCa = 0; end
    relErrConcK = in.n.delwpetarget{a}(15)./this_input.wpetarget(15); if ~isfinite(relErrConcK); relErrConcK = 0; end
    relErrConcTi = in.n.delwpetarget{a}(17)./this_input.wpetarget(17); if ~isfinite(relErrConcTi); relErrConcTi = 0; end
    relErrConcFe = in.n.delwpetarget{a}(20)./this_input.wpetarget(20);  if ~isfinite(relErrConcFe); relErrConcFe = 0; end
    
    % Fixed a bug here where the second term wasn't squared 20220620. 
    in.n.approxP.Psp_Ca(a) = in.profile{a}.P36sp_Ca(1).*consts.Psp_0_Ca_St.*in.n.sf_St(a);
    in.n.approxP.delPsp_Ca(a) = sqrt( (in.profile{a}.P36sp_Ca(1).*relErrConcCa.*consts.Psp_0_Ca_St.*in.n.sf_St(a)).^2 + (in.profile{a}.P36sp_Ca(1).*consts.delPsp_0_Ca_St.*in.n.sf_St(a)).^2 );
    
    in.n.approxP.Psp_K(a) = in.profile{a}.P36sp_K(1).*consts.Psp_0_K_St.*in.n.sf_St(a);
    in.n.approxP.delPsp_K(a) = sqrt( (in.profile{a}.P36sp_K(1).*relErrConcK.*consts.Psp_0_K_St.*in.n.sf_St(a)).^2 + (in.profile{a}.P36sp_K(1).*consts.delPsp_0_K_St.*in.n.sf_St(a)).^2 );
    
    in.n.approxP.Psp_Ti(a) = in.profile{a}.P36sp_Ti(1).*consts.Psp_0_Ti_St.*in.n.sf_St(a);
    in.n.approxP.delPsp_Ti(a) = sqrt( (in.profile{a}.P36sp_Ti(1).*relErrConcTi.*consts.Psp_0_Ti_St.*in.n.sf_St(a)).^2 + (in.profile{a}.P36sp_Ti(1).*consts.delPsp_0_Ti_St.*in.n.sf_St(a)).^2 );
    
    in.n.approxP.Psp_Fe(a) = in.profile{a}.P36sp_Fe(1).*consts.Psp_0_Fe_St.*in.n.sf_St(a);
    in.n.approxP.delPsp_Fe(a) = sqrt( (in.profile{a}.P36sp_Fe(1).*relErrConcFe.*consts.Psp_0_Fe_St.*in.n.sf_St(a)).^2 + (in.profile{a}.P36sp_Fe(1).*consts.delPsp_0_Fe_St.*in.n.sf_St(a)).^2 );
    
    temp_P_n = (in.profile{a}.P36sp_eth_Cl(1) + in.profile{a}.P36sp_th_Cl(1));
    in.n.approxP.Pn(a) = temp_P_n.*consts.P_f_0_St.*in.n.sf_St(a);
    in.n.approxP.delPn(a) = sqrt( ((temp_P_n.*in.n.delwpetarget{a}(13)./in.n.wpetarget{a}(13)).*consts.P_f_0_St.*in.n.sf_St(a)).^2 + (temp_P_n.*consts.delP_f_0_St.*in.n.sf_St(a)).^2 );
        
    % Muon production assumed zero uncertainty
    in.n.approxP.Pmu(a) = in.profile{a}.P36mu(1);
    
    in.n.approxP.P(a) = in.n.approxP.Psp_Ca(a)+in.n.approxP.Psp_K(a)+in.n.approxP.Psp_Ti(a)+in.n.approxP.Psp_Fe(a)+in.n.approxP.Pn(a)+in.n.approxP.Pmu(a);
    in.n.approxP.delP(a) = sqrt(in.n.approxP.delPsp_Ca(a).^2 + in.n.approxP.delPsp_K(a).^2 + in.n.approxP.delPsp_Ti(a).^2 + in.n.approxP.delPsp_Fe(a).^2 + in.n.approxP.delPn(a).^2);
    
    in.n.approxP.fn(a) = in.n.approxP.Pn(a)./in.n.approxP.P(a);
    
    % Checking uncertainty calcs
    %disp(['Ca ' sprintf('%0.2f',in.n.approxP.Psp_Ca(a)) '+/-' sprintf('%0.2f',in.n.approxP.delPsp_Ca(a)) '; K ' sprintf('%0.2f',in.n.approxP.Psp_K(a)) '+/-' sprintf('%0.2f',in.n.approxP.delPsp_K(a)) '; n: ' sprintf('%0.2f',in.n.approxP.Pn(a)) '+/-' sprintf('%0.2f',in.n.approxP.delPn(a)) '; total ' sprintf('%0.2f',in.n.approxP.P(a)) '+/-' sprintf('%0.2f',in.n.approxP.delP(a))]);
     
    if in.n.approxP.fn > 0.15
        % Flag excessive neutron capture production
        in.flags{a} =[in.flags{a} 'Sample ' in.s.sample_name{a} ': neutron capture production is ' sprintf('%0.2f',in.n.approxP.fn(a)) ' of total production. This could cause quite a lot of inaccuracy in production rate and exposure age estimates.<br>'];
    end
      
    % Also approx calculation of age. 
    % Note: at this point apply the shielding factor, even though this is
    % approximate because it shouldn't apply to muon production. 
    in.n.approxP.A(a) = consts.l36 + (in.s.E(in.n.index(a)).*in.s.rho(in.n.index(a)))./in.n.Lsp(a);
    tempN = in.n.N(a) - in.profile{a}.params.N_36_r_target; % radiogenic correction
    if tempN > (in.n.approxP.P(a).*in.s.othercorr(in.n.index(a)))./in.n.approxP.A(a)
        % Case apparently saturated
        in.n.approxP.t(a) = 2e6; % past saturation
    else        
        in.n.approxP.t(a) = (-1/(in.n.approxP.A(a).*in.s.othercorr(in.n.index(a)))).*log(1 - tempN.*in.n.approxP.A(a)./in.n.approxP.P(a));
    end
    % Also compute depth of erosion (cm)
    in.n.approxP.totalE(a) = in.s.E(in.n.index(a)).*in.n.approxP.t(a).*in.s.rho(in.n.index(a));
    
    if in.n.approxP.totalE(a) > (this_input.zbot./1.5)
        in.flags{a} = [in.flags{a} 'Sample ' in.s.sample_name{a} ': Total erosion may exceed ' sprintf('%0.0f',1.5.*in.n.approxP.totalE(a)) ' g cm-2,  which is outside the calculation domain. Can''t calculate exposure age. Try a lower erosion rate.<br>'];   
        in.ok_to_calculate(a) = 0;
    elseif in.n.approxP.totalE(a) > 20
        % Flag excessive implied erosion not enough to halt age calculation
        in.flags{a} = [in.flags{a} 'Sample ' in.s.sample_name{a} ': Total erosion is approx. ' sprintf('%0.0f',in.n.approxP.totalE(a)) ' g cm-2, which exceeds the expected accuracy range for muon production approximations.<br>'];
    end
    
    % Pass any other flags from production profile 
    for flags = 1:length(in.profile{a}.flags)
        in.flags{a} = [in.flags{a} 'Sample ' in.s.sample_name{a} ': ' in.profile{a}.flags{flags} '<br>'];
    end
    clear this_input;
end

% -------------- DONE COMPUTING PRODUCTION PROFILES --------------------


% ------------- COMPUTE AGES --------------------------------------------

% Initialize forward calculations

rin.lat = in.s.lat(in.n.index);
rin.long = in.s.long(in.n.index);
rin.yr = in.s.yr(in.n.index);
% Add buffer
rin.t = in.n.approxP.t.*1.5; % That should be enough. 
% Get Rc vector. This builds a new structure called sfdata. 
% Currently using dipole approximation in get_DipRc.m. 

sfdata = get_DipRc(rin,consts.RCfname,consts.Sfname);

% Now the structure sfdata contains a bunch of Rc vectors, one for each
% analysis. To compute scaling factors, we'll also need to know sample
% elevations as well as which nuclide we are dealing with. 

sfdata.pressure = in.s.pressure(in.n.index);
% This also needs to know latitudes to apply St scaling
sfdata.lat = rin.lat;

% Now sfdata is a structure with:
%   cell array sfdata.t contains time vectors
%   cell array sfdata.Rc contains Rc at corresponding times
%   cell array sfdata.S contains Lifton solar parameter 
%   cell array nuclide is nuclide identifier string
%   vector pressure is atmospheric pressure. Potentially this could be
%   a cell array if there is an isostatic rebound model upstream. 

% Get SF time series from precalculated results.
% Relies on get_StSF_Cl36.m, get_LmSF_Cl36.m, get_LSDnSF_Cl36.m
% Unfortunately these have to be different from non-Cl-36 case because of
% multiple production pathways. 

sfa = {'St' 'Lm' 'LSDn'}; 

for a = 1:length(sfa)
    eval(['sfdata = get_' sfa{a} 'SF_Cl36(sfdata,consts.datadir);'])
end   

% sfdata now has:
%   cell arrays Ksp_LSDn, Casp_LSDn, Pf0sp_LSDn
%   cell arrays Ksp_Lm, etc., etc.
% etc. These have vectors of "scaling factor". 

% Add sfdata to output structure so it is available later if needed

in.f = sfdata;

for a = 1:length(sfa) 
    % Do things for all scaling methods...
    % Initialize
    eval(['in.n.t_' sfa{a} ' = zeros(size(in.n.index));']);
    eval(['in.n.delt_int_' sfa{a} ' = in.n.t_' sfa{a} ';']);
    eval(['in.n.delt_ext_' sfa{a} ' = in.n.t_' sfa{a} ';']);
    
    for b = 1:length(in.n.index)
        % Now for each measurement
        if in.n.ok_to_calculate(b) == 1
            % Case ok to calculate...
            % We want to assemble the production profile from precalculated
            % components. 
            % Obtain vector of depths from erosion rate, time vector,
            % thickness -- note avg thickness, not effective
            these_zs = (in.s.rho(in.n.index(b)).*in.s.thick(in.n.index(b))./2) + (sfdata.tmin{b}+sfdata.tmax{b}).*in.s.E(in.n.index(b))./2; % Should be g/cm2
            % 
            % Multiply Ca, Fe, K, Ti by reference production rates to obtain
            % profiles for each reaction

            eval(['this_Psp_Ca = in.profile{b}.P36sp_Ca.*consts.Psp_0_Ca_' sfa{a} ';']);
            eval(['this_Psp_K = in.profile{b}.P36sp_K.*consts.Psp_0_K_' sfa{a} ';']);
            eval(['this_Psp_Fe = in.profile{b}.P36sp_Fe.*consts.Psp_0_Fe_' sfa{a} ';']);
            eval(['this_Psp_Ti = in.profile{b}.P36sp_Ti.*consts.Psp_0_Ti_' sfa{a} ';']);

            % Interpolate and apply scaling

            eval(['this_Psp_Ca_t = sfdata.Casp_' sfa{a} '{b}.*interp1(in.profile{b}.zz,this_Psp_Ca,these_zs);']);
            eval(['this_Psp_K_t = sfdata.Ksp_' sfa{a} '{b}.*interp1(in.profile{b}.zz,this_Psp_K,these_zs);']);
            eval(['this_Psp_Fe_t = sfdata.Fesp_' sfa{a} '{b}.*interp1(in.profile{b}.zz,this_Psp_Fe,these_zs);']);
            eval(['this_Psp_Ti_t = sfdata.Tisp_' sfa{a} '{b}.*interp1(in.profile{b}.zz,this_Psp_Ti,these_zs);']);

            % Sum and apply shielding correction
            this_Psp_t = (this_Psp_Ca_t + this_Psp_K_t + this_Psp_Fe_t + this_Psp_Ti_t).*in.s.othercorr(in.n.index(b)); 

            % That should be all spallogenic. 

            % Multiply neutron capture by Pf0, interpolate, apply scaling,
            % apply shielding

            eval(['this_Pn = in.profile{b}.P36n_Cl.*consts.P_f_0_' sfa{a} ';']);
            eval(['this_Pn_t = sfdata.Pf0sp_' sfa{a} '{b}.*interp1(in.profile{b}.zz,this_Pn,these_zs);']);

            this_Pn_t = this_Pn_t.*in.s.othercorr(in.n.index(b));

            % Interpolate muons/neutrons, no scaling
            % Adding K negative and fast, Ca negative and fast, 

            this_Pmu_t = interp1(in.profile{b}.zz,in.profile{b}.P36mu,these_zs);

            %
            % Add nucleogenic, no scaling obviously, and 
            % add up to obtain total production

            this_P36_t = this_Psp_t + this_Pn_t + this_Pmu_t + in.profile{b}.P36nuc;

            % That would appear to be total production.

            % Hybrid integral in decay only, not in erosion? This needs to be
            % checked.
            int_N = -(this_P36_t./consts.l36).*(exp(-sfdata.tmax{b}.*consts.l36) - exp(-sfdata.tmin{b}.*consts.l36));
            % That gives total production in each time step, corrected
            % for decay to the present time. 
            % Cumulate
            eval(['sfdata.cum_N_' sfa{a} '{b} = cumsum(int_N);']);

            % Saturation check
            eval(['tempmax = sfdata.cum_N_' sfa{a} '{b}(end);']);
            if (in.n.N(b) - in.profile{b}.params.N_36_r_target) > tempmax
                % Case saturated, no age calculation
                eval(['in.n.t_' sfa{a} '(b) = 0;']);
                in.flags{b} = [ in.flags{b} 'Sample ' char(in.s.sample_name(in.n.index(b))) ' appears to be saturated WRT ' char(sfa{a}) ' scaling at this erosion rate.<br>'];
            else
                % Case not saturated, subtract nucleogenic Cl-36 and compute age  
                eval(['in.n.t_' sfa{a} '(b) = interp1(sfdata.cum_N_' sfa{a} '{b},sfdata.tmax{b},(in.n.N(b) - in.profile{b}.params.N_36_r_target));']);
            end
          
        else
            % Case not OK to calculate
            % Zero ages already assigned, leave as is, do nothing
        end % end if ok to calculate
    end % end loop over measurements
end % end loop over scaling methods


% Error propagation loop goes here -- remember this can be vectorized, so
% should loop only over scaling methods

for a = 1:length(sfa)
    % Uncertainty computation
    % Linearized and approximated via S.A.E.
    
    % Obtain fraction external production rate uncertainty from approximate
    % calcs above
    % Note: strictly this is applicable only to St scaling. Necessary to do
    % this more times for other scaling? 
    fdelP = (in.n.approxP.delP./in.n.approxP.P)'; % note transpose
    % extract t for current SS
    eval(['tt = in.n.t_' sfa{a} ';']); 
         
    % Flag saturated or not calculated
	ok = find(tt > 0); notok = find(tt == 0 & in.n.ok_to_calculate == 1); 
    % Note A is always > 0, i.e. no stable nuclides
    % Preload zeros
	FP = zeros(size(in.n.index));
	delFP = FP; dtdN = FP; dtdP = FP;
    
    % Note transpose to prevent this messing up.       
	FP(ok) = (in.n.N(ok).*in.n.approxP.A(ok)')./(1 - exp(-in.n.approxP.A(ok)'.*tt(ok)));
    % Also, use these to compute effective normalized nuclide concentrations for
    % time-dependent scaling schemes, for messing around with in
    % banana plots. 
    eval(['in.n.Nnorm_' sfa{a} ' = zeros(size(in.n.index));']);
	eval(['in.n.delNnorm_' sfa{a} ' = in.n.Nnorm_' sfa{1} ';']);
    
    % do easy ones that aren't saturated
	eval(['in.n.Nnorm_' sfa{a} '(ok) = in.n.N(ok)./FP(ok);']);
    eval(['in.n.delNnorm_' sfa{a} '(ok) = in.n.delN(ok)./FP(ok);']);
	eval(['in.n.delNnorm_ext_' sfa{a} '(ok) = sqrt( (in.n.delN(ok)./FP(ok)).^2 + (fdelP(ok).*FP(ok).*in.n.N(ok)./(FP(ok).^2)).^2 );']);
 
	% If saturated, the production rate is the max N achieved in the
	% time dependent calculation times lambda. 
	for c = 1:length(notok)
        eval(['temp_sat_FP = sfdata.cum_N_' sfa{a} '{notok(c)}(end).*consts.l36;']);
        eval(['in.n.Nnorm_' sfa{a} '(notok(c)) = in.n.N(notok(c))./temp_sat_FP;']);
        eval(['in.n.delNnorm_' sfa{a} '(notok(c)) = in.n.delN(notok(c))./temp_sat_FP;']);
        eval(['in.n.delNnorm_ext_' sfa{a} '(notok(c)) = in.n.delNnorm_' sfa{a} '(notok(c));']); 
        % External uncert doesn't depend on P here? Why not? That needs to
        % be fixed. 
    end
    
    % Continue with error propagation.
	% If saturated, no error propagation. 
        
	delFP(ok) = fdelP(ok) .* FP(ok);
    dtdN(ok) = 1./(FP(ok) - in.n.N(ok).*in.n.approxP.A(ok)');  
	dtdP(ok) = -in.n.N(ok)./(FP(ok).*FP(ok) - in.n.N(ok).*in.n.approxP.A(ok)'.*FP(ok));
    % Note that uncertainty in the decay consts is not considered here.
    % This is just because it is basically never relevant in exposure
	% dating. To do this, however, one could observe that delA = dell 
	% and use dtdA. 
 
	% make respective delt's
	eval(['in.n.delt_ext_' sfa{a} '(ok) = sqrt( dtdN(ok).^2 .* in.n.delN(ok).^2 + dtdP(ok).^2 .* delFP(ok).^2);']);
	eval(['in.n.delt_int_' sfa{a} '(ok) = sqrt(dtdN(ok).^2 .* in.n.delN(ok).^2);']);
    
    clear ttt fdelP ok 2 FP delFP dtdN dtdP temp_sat_FP

end

% -------------- RETURN  ------------------------------------------------

ages = in;
    
        






