ET_objective.m

miss = ET_objective(E,cs,target,dFlag)

This is the objective function used by get_al_be_erosion to solve for the erosion rate.

The argument E is the erosion rate.

The argument target is the measured number of atoms to which the predicted value will be compared (atoms $\cdot$ g$^{-1}$). Thus, the value of E that yields miss = 0 solves the equation.

The argument cs is a structure with the following fields:

cs.tv	Time vector (yr)	1xn double array
cs.P_sp_t	Surface production rate due to spallation at times corresponding to those in cs.tv (atoms $\cdot$ g $^{-1}\cdot$ yr$^{-1}$)	1xn double array
cs.z_mu	Vector of depths. This is generally [0 logspace(0,5.3,100)] plus half the sample thickness. That is, we linearize $P(z)$ within the sample to get the thickness-averaged value of $P_{\mu}$ that is in cs.P_mu_z. (g $\cdot$ cm$^{-2}$)	1xm double array
cs.P_mu_z	Production rates due to muons at depths in cs.z_mu. (atoms $\cdot$ g $^{-1}\cdot$ yr$^{-1}$)	1xm double array
cs.l	Decay constant for nuclide of interest (yr$^{-1}$)	double
cs.tsf	Thickness scaling factor (nondimensional)	double
cs.L	Effective attenuation length for spallogenic production (g $\cdot$ cm$^{-2}$)	double

The argument dflag is a string variable telling the function what to return. If dflag = `no,' the output is just the objective function value, that is, the number of atoms difference between the measured value (target) and that predicted by the supplied erosion rate. If dflag = `yes,' the output is a structure containing diagnostic information, as follows:

miss.ver	Version number of this function	string
miss.Nmu	Nuclide concentration in sample attributable to production by muons at the given erosion rate (atoms $\cdot$ g$^{-1}$)	double
miss.Nsp	Nuclide concentration in sample attributable to production by spallation at the given erosion rate (atoms $\cdot$ g$^{-1}$)	double
miss.expected_Nsp, miss.PP, miss.A, miss.tsf	Assorted diagnostics. Some are optional - generally depending on whether time-dependent or steady production is being used.	various

This function calculates the following:

$$S_{thick}\int^{\infty}_{0} P_{sp,0}(t)e^{\left( -\epsilon t / \La... ...+ \int^{\infty}_{0} P_{\mu}(\epsilon t + \delta z /2) e^{-\lambda t} dt - N_{m}$$ (37)

where $z$ is depth (g $\cdot$ cm$^{-2}$), $\delta z$ is the sample thickness (g $\cdot$ cm$^{-2}$), $\epsilon$ is the erosion rate (the input argument E, here in g $\cdot$ cm $^{-2} \cdot$ yr$^{-1}$), $N_{m}$ is the measured nuclide concentration (the input argument target), $P_{sp}(t)$ is the surface production rate due to spallation at time $t$, $\lambda$ is the decay constant of the nuclide of interest, $\Lambda_{sp}$ is the effective attenuation length for production by spallation, $S_{thick}$ is the thickness scaling factor (dimensionless) and $P_{\mu}(z)$ is the thickness-averaged production rate due to muons as a function of depth.

This equation is basically the predicted nuclide concentration in the sample at the given erosion rate, less the measured nuclide concentration. Thus, the correct erosion rate is the value of E for which this function returns zero.

The first term of this equation, that is, the predicted nuclide concentration in the sample attributable to production by spallation, is integrated semi-numerically. We divide it into timesteps to account for surface production rate changes as a function of time, use trapezoidal integration for $P_{sp}(t)$ (that is, use the average production rate during each time step in that time step) and then integrate each time step analytically. So the nuclide concentration produced in a time step that extends from $T$ to $T+s$ is given by:

$$N_{T,T+s} = \frac{P_{sp}(T) + P_{sp}(T+s)}{2} \int^{T+s}_{T} e^{-\left( \lambda + \frac{\epsilon}{\Lambda} \right) t} dt$$ (38)

$$= -\frac{{P_{sp}(T) + P_{sp}(T+s)}}{2 \left( \lambda + \frac{\epsil... ... T+s \right)}- e^{-\left( \lambda + \frac{\epsilon}{\Lambda} \right) T} \right]$$ (39)

The final time step goes to infinity.

The second term, that is, the predicted nuclide concentration in the sample attributable to production by muons, must be calculated numerically. We do this by transforming the vector of depths at which we have calculated production due to muons to a vector of times (that is, dividing by the erosion rate), then using trapezoidal integration. The upper limit of integration $t_{max}$ is $(2 \times 10^{5}) / \epsilon$. The accuracy of this integration can be evaluated by similar arguments as described in the documentation for get_al_be_age.m above in section 1.3. The main difference here is that the effective attenuation length for production by muons increases with depth, so the step size can increase with time while still maintaining the same accuracy. Basically, the accuracy is set by the choice of the function argument cs.z_mu. We choose a logarithmic spacing that yields accuracy near or better than $10^{-3}$ for typical erosion rates.