ImpactChron

AsteroidHistory(typeseed<:Number; nnodes, Δt, tmax, downscale_factor)

struct containing Arrays that record the evolution of a (bombarded) asteroid.

Constructor function takes a parameter from the proposals to seed type and requires the number of radial nodes nnodes (::Int), the timestep used Δt, the full time duration tmax, and the downscale_factor (::Int).

Fields in AsteroidHistory:

ImpactSite(heat<:ImpactSiteShape, C<:Number)

struct describing the scale and shape of the simulated volume of impact heating.

CONSTRUCTOR FUNCTION

ImpactSite(shape, impactor_diameter)

Providing an impactor_diameter (<:Number) calculates impact parameters based on approximate heat distribution modeled in Davison+ 2012 (GCA, http://doi.org/10.1016/j.gca.2012.08.001).

ImpactSite(shape; r, C)

If values of r and C are provided, prepares an ImpactSite that extends to the center of an asteroid with radius r (km) and has a surface diameter C times the asteroid circumference (C ∈ [0,1], C = 0.01 by default). If no r is provided this seeds an ImpactSite with zeroed ImpactSiteShape and a C value.

ImpactSiteShape

Supertype of Cone, Parabola, and Hemisphere. Each contains a field z (depth) and r (radius), both in meters.

see also: ImpactSite

Nrm(μ::Float64,σ::Float64)

Immutable struct to describe normally distributed data, reported as mean (μ) and 1σ (σ)

PetroTypes(temps, samples)

struct containing fields of type3–type6, reflecting petrologic types. Each field contains a ImpactChron.TempProp with subfields T (maximum temperature in K) and p (proportion among the chondrite record). Note that ps sum to unity.

Only build PetroTypes with its constructor function

Constructor takes a NamedTuple containing maximum temperatures as fields T3-T6 and a Vector{String} containing the petrologic types corresponding to chondrites used as a prior.

e.g. PetroTypes( (T3 = 873, T4 = 973, T5 = 1073, T6 = 1223), ["4", "6", "3", "5,6", "im"])

If no argument given – PetroTypes() – returns zeroed type_ fields and weight=false, which short-circuits weighting by petrologic type.

PriorDistribution

Supertype of Nrm, lNrm, Unf

Unf(a::Float64,b::Float64)

Immutable struct to describe uniformly distributed data, reported as minimum (a) and maximum (b).

lNrm(μ::Float64,σ::Float64)

Immutable struct to describe lognormally distributed data, reported as natural-log-space mean (μ) and 1σ (σ)

ImpactChron.agerecal(x,sig;monitor_age,n, KAr=false)

Recalculate the age and uncertainty of an Ar-Ar age x±sig (1σ) Ma with the decay constants of Steiger & Jäger 1977 (EPSL, http://doi.org/10.1016/0012-821X(77)90060-7). Quantifies uncertainty by resampling n times (default: 10⁶). Optionally accepts a monitor_age. If not given, resamples from a range of absolute monitor ages ∈ [0,3] Ga.

If KAr=true, this recalculates the sample's K-Ar age with the new decay constants (no monitor correction necessary).

asteroid_agedist!(a::AsteroidHistory, p::NamedTuple, petrotypes::PetroTypes, impactsite::ImpactSite; 
    nᵣ,Tmax,Tmin, melt_reject)

Calculates and updates an asteroid thermochronologic history, contained in a, from parameters defined in p for an impact-heating morphology of impactsite. Weights the abundance of ages in each radial shell by petrotypes.

If any petrologic types are missing (see ImpactChron.weight_petro_types!) or the radial fraction of the body exceeding Tmax ≥ melt_reject (0.1 by default), the age distribution is zeroed, which ll_dist_params rejects.

In short, this is a convenient wrapper for a series of functions used in the weighted thermochronologic model.

see also: planetesimal_cooling_timestep!, planetesimal_cooling_dates,ImpactChron.weight_petro_types!,impact_reset_array!, ImpactChron.downscale!

csv2dict(filename::String;symbol::Bool=true)

Read a .csv file into a Dict. Assumes columns reflect discrete entries, with the first row specifying that entry's key. For NaN-buffered columns, the NaNs are removed, returning a single-element entry.

A true value for symbol will convert data saved as arrays of Strings into arrays of Symbols.

see also: nt2csv, csv2nt

csv2nt(filename::String;symbol=true)

Read a .csv file into a NamedTuple. Assumes columns reflect discrete entries, with the first row specifying that element's key. For NaN-buffered columns, the NaNs are removed, returning a single-element entry.

A true value for symbol will convert data saved as arrays of Strings into arrays of Symbols.

see also: nt2csv, csv2dict

data2csv(filename::String,data)

Save data to a .csv file. data may be one of the following: Dict, NamedTuple, or an Array, Accepts single-element or Vector fields. For single-element entries, the rest of the array/table is filled with NaNs.

see also: nt2csv, dict2csv

dict2csv(filename::String,D::Dict)

Save a Dict to a .csv. Accounts for fields of any type, keys must be Symbols.

For single-element entries, the rest of the array/table is filled with NaNs.

see also: nt2csv, data2csv

ImpactChron.downscale!(B, A)

Downscales elements of Array A into smaller Array B by summing. Scales the downscaled values in B by the ratio of length(A)÷length(B) to preserve any normalizations.

Requires that length(A) % length(B) = 0.

ImpactChron.draw(x)

Make a pseudorandom draw from x, which may be a Tuple or any subtype of PriorDistribution. If x is a Number, it simply returns x.

Used in support of lognormMC.

See also: Nrm, lNrm, Unf, lognormMC

ImpactChron.histogramify!(dist::AbstractVector, domain::AbstractRange, x::AbstractVector, y::AbstractVector)

In-place histogramify that overwites a pre-allocated vector dist. See histogramify for details.

ImpactChron.histogramify(domain::AbstractRange,x::AbstractVector,y::AbstractVector)

Constructs histogram over bins defined by domain from model outputs in x with corresponding abundances in y. Does not require a constant step in x, so this calculates histograms of outputs from planetesimal_cooling_dates.

histogramify normalizes its output, such that for output dist ∑ dist * Δd = 1, so long as all x ∈ domain. If any x ∉ domain, ∑ dist * Δd = 1- (∑yₒᵤₜ / ∑yₐₗₗ ) where the xₒᵤₜ corresponding to yₒᵤₜ are ∉ domain.

Returns only the histogram masses, centers of bins must be calculated separately with rangemidpoints(domain)

impact_reset_array!(
    tₓr::AbstractArray, solartime::AbstractArray, 
    tcoolₒ::AbstractArray, Vfrxn::AbstractArray,                             
    impacts::AbstractArray, p::NamedTuple, c::NamedTuple;
    nᵣ::Integer, Δt::Number )

Simulates an impact history from -χ parameters in p (see below), and resets primary cooling dates (indices of dates in solartime in tcoolₒ) and fractional volumes (Vfraxn) based on impact/crater properties described in c.

Relative abundances of depth-age pairs are tracked in array tₓr (time x radial depth), with dimensions (length(solartime),nᵣ), where nᵣ describes the number of radial nodes, as in planetesimal_cooling_dates.

impacts and tcoolₒ are pre-allocated vectors that respectively record the number of impacts at each time step and the index of the primary cooling date in solartime.

Impact flux follows an exponential decay described by parameters in p:

ll_dist(x, dist, mu, sigma)

Calculate loglikelihood that observations in Vectors mu and sigma are drawn from modeled distribution described by Vectors x and dist, where x contains the bincenters of a normalized histogram dist and mu and sigma respectively contain the mean and 1σ of the observations.

ll_dist_params(a, p, plims, mu,sig)

Calculate the combined log-likelihood that

1. Observations with mean mu and corresponding 1σ uncertainties sig (both ::Array) are drawn from the downscaled age distribution and timesteps contained in a (::AsteroidHistory)

2. The proposal parameters p used to calculate this history are drawn from the prior distributions plims (both ::NamedTuple)

If the age distribution is zero, quickly returns a log-liklihood of -Inf.

ll_param(x::Number,D<:PriorDistribution)

Calculate the log-likelihood that x is drawn from a distribution D, which may be a Normal (Nrm), Lognormal (lNrm), or Uniform (Unf) distribution. For D::Unf, ll_param returns 0 and bounds checks are performed by ImpactChron.prior_bounds within ImpactChron.strict_priors.

ll_params(p::NamedTuple,d::NamedTuple)

Calculate log-likelihoods for the proposals in p with corresponding distributions in d

Does not evaluate impact (χ) parameters, since these were only modeled on Unf distributions.

loadArAr(maxage)

Loads and returns all Ar-Ar ages in the ImpactChron database older than maxage (in Ma) as a tuple of ages, ages_σ, types, groups, names.

see also: [`loadKAr`](@ref), [`loadKArAr`](@ref)

loadECs(maxage)

Loads and returns all Ar-Ar ages for enstatite chondrites (ECs) in the ImpactChron database older than maxage (in Ma) as a tuple of ages, ages_σ, types, groups, names.

see also: [`loadArAr`](@ref), [`loadTrie2003`](@ref)

loadKAr(maxage)

Loads and returns all K-Ar ages in the ImpactChron database older than maxage (in Ma) as a tuple of ages, ages_σ, types, groups, names.

see also: [`loadArAr`](@ref), [`loadKArAr`](@ref)

loadKArAr(maxage)

Loads and returns all K-Ar and Ar-Ar ages in the ImpactChron database older than maxage (in Ma) as a tuple of ages, ages_σ, types, groups, names.

see also: [`loadArAr`](@ref), [`loadKAr`](@ref)

loadTrie2003(maxage)

Loads and returns Ar-Ar ages for the H chondrites chosen to reflect an unperturbed onion shell by Trieloff+ (2003, doi:10.1038/nature01499) in the ImpactChron database older than maxage (in Ma) as a tuple of ages, ages_σ, types, groups, names.

see also: [`loadArAr`](@ref), [`loadTrie2003`](@ref)

ImpactChron.lognorm(x)

Calculate the lognormal distribution of x from a collection of Numbers (accepts NamedTuples). Returns a lNrm type.

ImpactChron.lognormMC(x ; n)

Calculate the lognormal distribution of a collection x by resampling the entire collection n times (default=10⁶).

x may contain data in the form of Tuples, PriorDistribution subtypes, or Numbers.

Returns a lNrm type.

See also: ImpactChron.lognorm, ImpactChron.draw

Just in case, the function has a (slow) safety net to prevent it from trying to calculate the log of any negative resamples. (This has never happened for the data I used)

ImpactChron.mcmean(x, xsig, n=10_000_000, fullpost=false)

Calculate mean of age(s) x with 1σ uncertainties xsig by Monte Carlo method. A 10⁷ (default) resample routine returns consistent results at the 10⁻⁴ level.

Returns a tuple of (mean,1σ) by default. If fullpost=true, returns posterior samples rather than summary statistics.

perturb(p::NamedTuple,k::Symbol,n::Number)

Return a NamedTuple identical to p, with one field (key k) changed to the value of n. Note that == identity is preserved only if the order of fields in p is: tss,rAlo,R,ta,cAl,Tm,Tc,ρ,Cp,k,tχα,τχα,Fχα,tχβ,τχβ,Fχβ,tχγ,τχγ,Fχγ

planetesimal_cooling_dates!(ages, Vfrxn, peakT, p::NamedTuple; nᵣ, Δt, tmax, Tmax, Tmin)

In-place planetesimal_cooling_dates that updates Arrays ages, Vfrxn, and peakT.

see also: planetesimal_cooling_dates

planetesimal_cooling_dates(p::NamedTuple; nᵣ, Δt, tmax, Tmax, Tmin)

Returns an NTuple containing (in this order) thermochronologic cooling dates (My after CAIs) and corresponding volume fractions, radial depths (km from center), and peak temperatures (K) of nᵣ evenly spaced nodes in a spherical body.

Physical and environmental parameters are described in p. Alternatively, these parameters may be individually listed in lieu of p. These parameters are outlined in the table below.

Note: several of these parameters need to be entered as the natural logarithm of the value for easy compatibility with the inversion function.

Δt gives the timestep (in My), tmax describes the duration of the model (My after CAIs = Myₛₛ), and Tmax and Tmin define the maximum and minimum temperatures (K) allowed for chondritic material in the body. Default values are only given for tmax (2000 Myₛₛ), Tmax (1500 K), and Tmin (0 K).

| Parameter                 | log?  | `NmTpl`| `func`  |
| :------------------------ | :--:  | :----: | :-----: |
| solar system age (Ma)     | no    | `tss`  | `tₛₛ`    |
| initial ²⁶Al/²⁷Al         | no    | `rAlo` | `rAlo`   |
| closure temperature (K)   | yes   | `Tc`   | `Tc`     |
| body radius (m)           | yes   | `R`    | `R`      |
| accretion date (Myₛₛ)     | yes   | `ta`   | `tₐ`     |
| disk temperature (K)      | yes   | `Tm`   | `To`     |
| [Al] (g/g)                | yes   | `cAl`  | `Al_conc`|
| density (kg/m³)           | yes   | `ρ`    | `ρ`      |
| thermal diffusivity       | yes   | `k`    | `K`      |
| specific heat capacity    | yes   | `Cp`   | `Cₚ`     |
| ------------------------- | ----  | ------ | -------- |

see also: planetesimal_cooling_dates!

ImpactChron.planetesimal_cooling_timestep!(solartime::AbstractRange, time_i::Vector, Vfrxn::Vector, peakT::Vector, p::NamedTuple; nᵣ, Tmax, Tmin)

Returns (overwrites) thermochronologic cooling dates in time_i as indices of solartime, along with corresponding volumetric fractions (Vfrxn) and peak temperatures in K (peakT) for nᵣ nodes in a body with model parameters given in p. Tmax and Tmin respectively describe the maximum and minimum peak temperatures allowed. Failing to exceed Tmin returns the date of accretion, and exceeding Tmax sets the volumetric fraction to zero (achondritic).

see also: planetesimal_cooling_dates, planetesimal_cooling_dates!

ImpactChron.planetesimal_temperature(time::AbstractArray, radii::AbstractArray; To, Ao, λ, K, κ)

Calculates the evolution of temperature over time steps at a range of radii defined for a conductively cooling sphere with thermal conductivity K and thermal diffusivity κ, given ambient temperature To, initial heat production Ao, and heat-production decay constant λ.

last(radii) defines the radius of the sphere.

Adapted from: Carlslaw & Jäger (1959, ISBN-13: 978-0198533689) and Hevey & Sanders (2006)

ImpactChron.prior_bounds(x, p<:PriorDistribution)

Evaluates whether a proposal x falls within the permissible bounds of its prior p <: PriorDistribution. Always returns true if p is of type Nrm or lNrm, tests if x ∈ (p.a,p.b) for p::Unf.

ImpactChron.radius_at_depth(rᵢ, R, x<:ImpactSiteShape)

Calculates the radius of the circle traced by a x-shaped region at a distance of rᵢ from the center of a body of radius R. Note that for x::Hemisphere, only its x.r is used.

see also: ImpactSiteShape

rangebinbounds(x)

Calculate a LinRange of the linear bounds for each "midpoint" step in x (<:AbstractRange).

rangemidpoints(x)

Calculate a LinRange of the midpoints of each step in x(<:AbstractRange).

serial2dict(file::String, vars::Tuple; n, ll=true, accept=true, perturbation=false)

Load a serialized archive file from thermochron_metropolis into a Dict for post-run analysis, only incorporating the variables in vars and covering n steps (all steps by default). Provide a Bool to incorporate log-likelihoods ll, acceptances accept, or the perturbed variables perturbation at each step.

ImpactChron.strict_priors(p, k, p_prior<:PriorDistribution)

Evaluates strict priors related to Uniform distributions and other rules for paramter proposal p (::NamedTuple) and perturbed variable k (::Symbol), where p_prior is a PriorDistribution of p[k].

Returns true if all priors are satisfied. Returns false if any priors fail.

Currently includes:

1. Ensure bounds of uniform priors are not exceeded. See ImpactChron.prior_bounds.

2. Ensure bombardment events α, β, γ begin in sequential order.

3. The ℯ-folding time of bombardment α must be longer than those of β and γ.

function thermochron_metropolis(p, pσ, pvars, mu, sigma, impactsite; nsteps, plims, downscale, petrotypes, kwargs...)

Runs a Markov chain Monte Carlo (MCMC) routine that explores the parameter space of the variables in pvars (as type Symbol, see PARAMETERS table below), constrained by observed thermochronologic ages given in mu with corresponding 1σ uncertainties in sigma. Requires an initial proposal p (::NamedTuple) for parameter values and a Gaussian jump size pσ (::NamedTuple), each containing all model parameters listed in the table below. Note that many parameters are log-normally distributed and require a natural-log-space initial guess (see table). impactsite (::ImpactSite) describes the shape of asteroidal volume reheated per impact. See ImpactSite for details.

Four kwargs are particularly important:

1. nsteps ::Int is the number of post-burn-in Markov chain steps.

2. plims ::NamedTuple – The model relies on several priors –- all non-bombardment variables (χ ∉ parameter name) are constrained by prior distributions.

Field names are as in p and pσ and values are subtypes of PriorDistribution.

1. downscale ::Integer bins the modeled age distribution to smooth it. The distribution is "downscaled" by the factor downcale. It is set to 1 by default (off). I recommend a value of 10. See ImpactChron.downscale! for details.

2. petrotypes ::PetroTypes contains the weights and upperbound temperatures of different petrologic types. This weighting is turned off by default. See PetroTypes for details.

PARAMETERS:

All the other kwargs...

ImpactChron.timemanagement(Δt, tmax, downscale_factor::Int)

Calculate (and adjust if necessary) the model timescales for a given tmax (My after CAIs), Δt timestep (in My), and a downscale_factor. May overwrite tmax to ensure downscaling works properly.

Returns a Tuple containing the timescale and the downscaled timescale.

trim_ages(maxage, data)

Helper function to remove all dates younger than maxage (in Ma) in K/Ar-Ar dataloaded internally from ImpactChron. Used under the hood in [loadArAr](@ref), [loadKAr](@ref), [loadKArAr](@ref), [loadECs](@ref), [loadTrie2003`](@ref)

ImpactChron.weight_petro_types!(v, T, petrotypes::PetroTypes)

Reweight volumetric fractions relative to the abundance of each petrologic type in the meteorite record. Takes Arrays of volumetric fraction (v) and peak temperature (T) as output by planetesimal_cooling_dates or contained within an AsteroidHistory.

Requires all petrologic types to occupy at least one radial node, otherwise returns zero in all v.

see also: PetroTypes