function out = BCO_objective_steadyerosion(X,d,p,dFlag)

% This is the objective function for fitting exercise 2, which is assuming
% steady erosion for 14.5 Ma. 
%
% out = BCO_objective_steadyerosion(X,d,p,dFlag)
% 
% where X has the following params:
%   X(1) erosion rate (m/Myr)
%   X(2) Ne closure age (Myr)
%   X(3) Lambda_a (g/cm2) albedo e-folding length
%   X(4) f_a (nondimensional) albedo loss percent
%   X(5) fstar (nondimensonal * 1e-3) total negative muon capture likelihood
%   X(6) sigma0 (millibarns) fast muon cross-section
%
% d is data structure with core data
% p is additional struct with things like surface production rate, etc. 
% dflag is 0 for objective fctn; 1 for diagnostic struct
%
% Greg Balco - Berkeley Geochronology Center - 2018
% Not licensed for reuse or distribution

% Get in correct dimensions

ee = X(1).*p.rho.*1e-4; % to g/cm2/yr
tc = X(2).*1e6; % to Ma
La = X(3); % No conversion
fa = X(4).*0.01; % convert percent to fractional
fstar = X(5).*1e-3; % to correct fractional
sigma0 = X(6).*1e-27; % to cm2

% Nucleogenic

% Alpha production
d.a232perg = 6.*d.atsThperg.*(exp(p.l232.*tc) - 1);
d.a235perg = 7.*d.ats235perg.*(exp(p.l235.*tc) - 1);
d.a238perg = 8.*d.ats238perg.*(exp(p.l238.*tc) - 1);

% Ne-21 production
r.N21nuc = d.Ft.*( d.a232perg.*6.08e-8 + d.a235perg.*5.62e-8 + d.a238perg.*4.04e-8 );
r.delN21nuc = r.N21nuc.*p.dnfactor;

% Spallogenic

r.N21sp = exp(-d.z./p.Lsp).*(p.P0.*p.Lsp./ee).*(1-exp(-p.texp.*(ee./p.Lsp))) - exp(-d.z./La).*(fa.*p.P0.*La./ee).*(1-exp(-p.texp.*(ee./La)));

% Muons

mc.Natoms = p.Natoms;
mc.sigma0 = sigma0;
mc.k_neg = fstar.*p.fC.*p.fD;


for a = 1:length(d.N21m);
	r.N21fast(a) = integral(@(t) sigma0.*interp1(log(p.predz),p.mstub_fast,log(d.z(a)+ee.*t)),0,p.texp,'reltol',1e-3,'abstol',1e-1);
    r.N21neg(a) = integral(@(t) fstar.*interp1(log(p.predz),p.mstub_neg,log(d.z(a)+ee.*t)),0,p.texp,'reltol',1e-3,'abstol',1e-1);
end;


% Total

r.N21p = r.N21nuc + r.N21sp + r.N21fast' + r.N21neg';
r.delN21p = r.delN21nuc;
r.miss = ((r.N21p - d.N21m)./sqrt((d.delN21m.*p.dmfactor).^2 + r.delN21p.^2));
r.missquared = r.miss.^2;
r.x2 = sum(r.missquared);

if dFlag == 0;
    out = r.x2;
elseif dFlag == 1;
    out = r;
end;

