diff --git a/examples/ZrnHeInversionVartCryst.jl b/examples/ZrnHeInversionVartCryst.jl
index 8405a28..eccca44 100755
--- a/examples/ZrnHeInversionVartCryst.jl
+++ b/examples/ZrnHeInversionVartCryst.jl
@@ -50,8 +50,8 @@
     dTmax = 10.0 # Maximum reheating/burial per model timestep (to prevent numerical overflow)
 
     # Other model parameters
-    CrystAgeMax_Ma = 4000.0 # Ma -- forbid anything older than this
-    TCryst = 400.0 # Temperature (in C)
+    CrystAgeMax = 4000.0 # Ma -- forbid anything older than this
+    TStart = 400.0 # Temperature (in C)
     dr = 1 # Radius step, in microns
 
     diffusionparams = (;
@@ -82,49 +82,52 @@
 
     # Populate local variables from data frame with specified options
     Halfwidth = data.Halfwidth_um
-    U_ppm = data.U238_ppm
-    Th_ppm = data.Th232_ppm
-    HeAge_Ma = data.HeAge_Ma_raw
-    HeAge_Ma_sigma = data.HeAge_Ma_sigma_raw
-    CrystAge_Ma = data.CrystAge_Ma
+    Uppm = data.U238_ppm
+    Thppm = data.Th232_ppm
+    HeAge = data.HeAge_Ma_raw
+    HeAge_sigma = data.HeAge_Ma_sigma_raw
+    CrystAge = data.CrystAge_Ma
 
-    CrystAge_Ma[CrystAge_Ma .> CrystAgeMax_Ma] .= CrystAgeMax_Ma
-    tCryst = ceil(maximum(CrystAge_Ma)/dt) * dt
+    CrystAge[CrystAge .> CrystAgeMax] .= CrystAgeMax
+    tStart = ceil(maximum(CrystAge)/dt) * dt
 
-    tSteps = Array{Float64,1}(0+dt/2 : dt : tCryst-dt/2)
+    tSteps = Array{Float64,1}(0+dt/2 : dt : tStart-dt/2)
     ageSteps = reverse(tSteps)
     ntSteps = length(tSteps) # Number of time steps
-    eU = U_ppm+.238*Th_ppm # Used only for plotting
+    eU = Uppm+.238*Thppm # Used only for plotting
 
 ## --- Test proscribed t-T paths with Neoproterozoic exhumation step
 
-    # Generate T path to test
+    # # Generate T path to test
     # Tr = 150
     # T0 = 30
-    # agePoints = Float64[tCryst, tCryst*29/30,   720, 580, 250,  0] # Age (Ma)
-    # TPoints  =  Float64[TCryst,        Tr+T0, Tr+T0,  T0,  70, 10] # Temp. (C)
+    # agePoints = Float64[tStart, tStart*29/30,   720, 580, 250,  0] # Age (Ma)
+    # TPoints  =  Float64[TStart,        Tr+T0, Tr+T0,  T0,  70, 10] # Temp. (C)
     # TSteps = linterp1s(agePoints,TPoints,ageSteps)
-
-    # Plot t-T path
-    #plot(ageSteps,TSteps,xflip=true)
-
-    # Calculate model ages
-    #CalcHeAges = Array{Float64}(undef, size(HeAge_Ma))
-    #@time pr = DamageAnnealing(dt,tSteps,TSteps)
-    #@time for i=1:length(Halfwidth)
-    #    first_index = 1 + floor(Int64,(tCryst - CrystAge_Ma[i])/dt)
-    #    CalcHeAges[i] = ZrnHeAgeSpherical(dt,ageSteps[first_index:end],TSteps[first_index:end],pr[first_index:end,first_index:end],Halfwidth[i],dr,U_ppm[i],Th_ppm[i],diffusionparams)
-    #end
-
-    # Plot Comparison of results
-    #p2 = plot(eU, CalcHeAges, seriestype=:scatter,label="Model")
-    #plot!(p2, eU, HeAge_Ma, yerror=HeAge_Ma_sigma*2, seriestype=:scatter, label="Data")
-    #xlabel!(p2,"eU"); ylabel!(p2,"Age (Ma)")
-    #display(p2)
-
-    # Check log likelihood
-    #ll = sum(-(CalcHeAges - HeAge_Ma).^2 ./ (2 .* HeAge_Ma_sigma.^2))
-    #ll = sum(-(CalcHeAges - HeAge_Ma).^2 ./ (2 .* AnnealedSigma.^2))
+    #
+    # # Plot t-T path
+    # plot(ageSteps,TSteps,xflip=true)
+    #
+    # # Calculate model ages
+    # CalcHeAges = Array{Float64}(undef, size(HeAge))
+    # @time pr = DamageAnnealing(dt,tSteps,TSteps)
+    # zircons = Array{Zircon{Float64}}(undef, length(Halfwidth))
+    # @time for i=1:length(zircons)
+    #     # Iterate through each grain, calculate the modeled age for each
+    #     first_index = 1 + floor(Int64,(tStart - CrystAge[i])/dt)
+    #     zircons[i] = Zircon(Halfwidth[i], dr, Uppm[i], Thppm[i], dt, ageSteps[first_index:end])
+    #     CalcHeAges[i] = HeAgeSpherical(zircons[i], @views(TSteps[first_index:end]), @views(pr[first_index:end,first_index:end]), diffusionparams)
+    # end
+    #
+    # # Plot Comparison of results
+    # p2 = plot(eU, CalcHeAges, seriestype=:scatter,label="Model")
+    # plot!(p2, eU, HeAge, yerror=HeAge_sigma*2, seriestype=:scatter, label="Data")
+    # xlabel!(p2,"eU"); ylabel!(p2,"Age (Ma)")
+    # display(p2)
+    #
+    # # Check log likelihood
+    # ll = sum(-(CalcHeAges - HeAge).^2 ./ (2 .* HeAge_sigma.^2))
+    # ll = sum(-(CalcHeAges - HeAge).^2 ./ (2 .* AnnealedSigma.^2))
 
 
 ## --- Invert for maximum likelihood t-T path
@@ -132,8 +135,8 @@
     # Boundary conditions (10C at present and 650 C at the time of zircon
     # formation). Optional: Apppend required unconformities at base of Cambrian,
     # or elsewhere
-    boundary_agePoints = Float64[0, tCryst] # Ma
-    boundary_TPoints = Float64[0, TCryst] # Degrees C
+    boundary_agePoints = Float64[0, tStart] # Ma
+    boundary_TPoints = Float64[0, TStart] # Degrees C
 
     # This is where the "transdimensional" part comes in
     nPoints = 0
@@ -143,7 +146,7 @@
 
     # Fill some intermediate points to give the MCMC something to work with
     Tr = 250 # Residence temperature of initial proposal (value should not matter too much)
-    for t in [tCryst/30,tCryst/4,tCryst/2,tCryst-tCryst/4,tCryst-tCryst/30]
+    for t in [tStart/30,tStart/4,tStart/2,tStart-tStart/4,tStart-tStart/30]
         global nPoints += 1
         agePoints[nPoints] = t # Ma
         TPoints[nPoints] = Tr # Degrees C
@@ -152,7 +155,7 @@
     # # (Optional) Start with something close to the expected path
     #Tr = 150
     #T0 = 30
-    #agePoints[1:4] = Float64[tCryst*29/30,   720, 510, 250]) # Age (Ma)
+    #agePoints[1:4] = Float64[tStart*29/30,   720, 510, 250]) # Age (Ma)
     #TPoints[1:4]  =  Float64[       Tr+T0, Tr+T0,  T0,  70]) # Temp. (C)
     #nPoints+=4
 
@@ -160,18 +163,20 @@
         # Calculate model ages for initial proposal
         TSteps = linterp1s([view(agePoints, 1:nPoints) ; boundary_agePoints ; unconf_agePoints],
                            [view(TPoints, 1:nPoints) ; boundary_TPoints ; unconf_TPoints], ageSteps)
-        CalcHeAges = Array{Float64}(undef, size(HeAge_Ma))
+        CalcHeAges = Array{Float64}(undef, size(HeAge))
         pr = DamageAnnealing(dt,tSteps,TSteps) # Damage annealing history
-        for i=1:length(Halfwidth)
+        zircons = Array{Zircon{Float64}}(undef, length(Halfwidth))
+        for i=1:length(zircons)
             # Iterate through each grain, calculate the modeled age for each
-            first_index = 1 + floor(Int64,(tCryst - CrystAge_Ma[i])/dt)
-            CalcHeAges[i] = ZrnHeAgeSpherical(dt,ageSteps[first_index:end],TSteps[first_index:end],pr[first_index:end,first_index:end],Halfwidth[i],dr,U_ppm[i],Th_ppm[i], diffusionparams)
+            first_index = 1 + floor(Int64,(tStart - CrystAge[i])/dt)
+            zircons[i] = Zircon(Halfwidth[i], dr, Uppm[i], Thppm[i], dt, ageSteps[first_index:end])
+            CalcHeAges[i] = HeAgeSpherical(zircons[i], @views(TSteps[first_index:end]), @views(pr[first_index:end,first_index:end]), diffusionparams)
         end
-        AnnealedSigma = simannealsigma.(1, HeAge_Ma_sigma; params=simannealparams)
-        UnAnnealedSigma = simannealsigma.(nsteps, HeAge_Ma_sigma; params=simannealparams)
+        AnnealedSigma = simannealsigma.(1, HeAge_sigma; params=simannealparams)
+        UnAnnealedSigma = simannealsigma.(nsteps, HeAge_sigma; params=simannealparams)
 
         # Log-likelihood for initial proposal
-        ll = normpdf_ll(HeAge_Ma, AnnealedSigma, CalcHeAges)
+        ll = normpdf_ll(HeAge, AnnealedSigma, CalcHeAges)
         if simplified
             ll -= log(nPoints)
         end
@@ -188,15 +193,15 @@
         CalcHeAgesₚ = similar(CalcHeAges)
 
         # Distributions to populate
-        HeAgeDist = Array{Float64}(undef, length(HeAge_Ma), nsteps)
+        HeAgeDist = Array{Float64}(undef, length(HeAge), nsteps)
         TStepDist = Array{Float64}(undef, ntSteps, nsteps)
         llDist = Array{Float64}(undef, nsteps)
         nDist = zeros(Int, nsteps)
         acceptanceDist = zeros(Bool, nsteps)
 
         # Standard deviations of Gaussian proposal distributions for temperature and time
-        t_sigma = tCryst/60
-        T_sigma = TCryst/60
+        t_sigma = tStart/60
+        T_sigma = TStart/60
 
         # Proposal probabilities (must sum to 1)
         move = 0.64
@@ -227,17 +232,17 @@
                     if agePointsₚ[k] < dt
                         # Don't let any point get too close to 0
                         agePointsₚ[k] += (dt - agePointsₚ[k])
-                    elseif agePointsₚ[k] > (tCryst - dt)
-                        # Don't let any point get too close to tCryst
-                        agePointsₚ[k] -= (agePointsₚ[k] - (tCryst - dt))
+                    elseif agePointsₚ[k] > (tStart - dt)
+                        # Don't let any point get too close to tStart
+                        agePointsₚ[k] -= (agePointsₚ[k] - (tStart - dt))
                     end
                     # Move the Temperature of one model point
                     if TPointsₚ[k] < 0
                         # Don't allow T<0
                         TPointsₚ[k] = 0
-                    elseif TPointsₚ[k] > TCryst
-                        # Don't allow T>TCryst
-                        TPointsₚ[k] = TCryst
+                    elseif TPointsₚ[k] > TStart
+                        # Don't allow T>TStart
+                        TPointsₚ[k] = TStart
                     end
 
                     # Interpolate proposed t-T path
@@ -255,8 +260,8 @@
                 # Birth: add a new model point
                 nPointsₚ += 1
                 for i=1:maxattempts # Try maxattempts times to satisfy the reheating rate limit
-                    agePointsₚ[nPointsₚ] = rand()*tCryst
-                    TPointsₚ[nPointsₚ] = rand()*TCryst
+                    agePointsₚ[nPointsₚ] = rand()*tStart
+                    TPointsₚ[nPointsₚ] = rand()*TStart
 
                     # Interpolate proposed t-T path
                     TStepsₚ = linterp1s([view(agePointsₚ, 1:nPointsₚ) ; boundary_agePoints ; unconf_agePointsₚ],
@@ -288,8 +293,8 @@
                         # Allow the present temperature to vary from 0 to 10 degrees C
                         boundary_TPointsₚ[1] = 0+rand()*10
                     else
-                        # Allow the initial temperature to vary from TCryst to TCryst-50 C
-                        boundary_TPointsₚ[2] = TCryst-rand()*50
+                        # Allow the initial temperature to vary from TStart to TStart-50 C
+                        boundary_TPointsₚ[2] = TStart-rand()*50
                     end
                     if length(unconf_agePointsₚ) > 0
                         # If there's an imposed unconformity, adjust within parameters
@@ -317,15 +322,15 @@
 
              # Calculate model ages for each grain
             DamageAnnealing!(pr, dt, tSteps, TStepsₚ)
-            for i=1:length(Halfwidth)
-                first_index = 1 + floor(Int64,(tCryst - CrystAge_Ma[i])/dt)
-                @views CalcHeAgesₚ[i] = ZrnHeAgeSpherical(dt,ageSteps[first_index:end],TStepsₚ[first_index:end],pr[first_index:end,first_index:end],Halfwidth[i],dr,U_ppm[i],Th_ppm[i],diffusionparams)
+            for i=1:length(zircons)
+                first_index = 1 + floor(Int64,(tStart - CrystAge[i])/dt)
+                CalcHeAgesₚ[i] = HeAgeSpherical(zircons[i], @views(TStepsₚ[first_index:end]), @views(pr[first_index:end,first_index:end]), diffusionparams)
             end
 
             # Calculate log likelihood of proposal
-            AnnealedSigma .= simannealsigma.(n, HeAge_Ma_sigma; params=simannealparams)
-            llₚ = normpdf_ll(HeAge_Ma, AnnealedSigma, CalcHeAgesₚ)
-            llₗ = normpdf_ll(HeAge_Ma, AnnealedSigma, CalcHeAges) # Recalulate last one too with new AnnealedSigma
+            AnnealedSigma .= simannealsigma.(n, HeAge_sigma; params=simannealparams)
+            llₚ = normpdf_ll(HeAge, AnnealedSigma, CalcHeAgesₚ)
+            llₗ = normpdf_ll(HeAge, AnnealedSigma, CalcHeAges) # Recalulate last one too with new AnnealedSigma
             if simplified # slightly penalize more complex t-T paths
                 llₚ -= log(nPointsₚ)
                 llₗ -= log(nPoints)
@@ -351,7 +356,7 @@
             end
 
             # Record results for analysis and troubleshooting
-            llDist[n] = normpdf_ll(HeAge_Ma, UnAnnealedSigma, CalcHeAges) # Recalculated to constant baseline
+            llDist[n] = normpdf_ll(HeAge, UnAnnealedSigma, CalcHeAges) # Recalculated to constant baseline
             nDist[n] = nPoints # Distribution of # of points
             HeAgeDist[:,n] = CalcHeAges # Distribution of He ages
 
@@ -366,11 +371,11 @@
     @time (TStepDist, HeAgeDist, nDist, llDist, acceptanceDist) = MCMC_vartcryst(nPoints, maxPoints, agePoints, TPoints, unconf_agePoints, unconf_TPoints, boundary_agePoints, boundary_TPoints, simannealparams, diffusionparams)
 
     # # Save results using JLD
-    @save string(name, ".jld") ageSteps tSteps TStepDist burnin nsteps TCryst tCryst
+    @save string(name, ".jld") ageSteps tSteps TStepDist burnin nsteps TStart tStart
 
     # Plot sample age-eU correlations
     h = scatter(eU,HeAgeDist[:,burnin:burnin+50], label="")
-    plot!(h, eU, HeAge_Ma, yerror=HeAge_Ma_sigma, seriestype=:scatter, label="Data")
+    plot!(h, eU, HeAge, yerror=HeAge_sigma, seriestype=:scatter, label="Data")
     xlabel!(h,"eU (ppm)"); ylabel!(h,"Age (Ma)")
     savefig(h,string(name,"Age-eU.pdf"))
 
@@ -381,14 +386,14 @@
 
     # Resize the post-burnin part of the stationary distribution
     TStepDistResized = Array{Float64}(undef, 2001, size(TStepDist,2)-burnin)
-    xq = collect(range(0,tCryst,length=2001))
+    xq = collect(range(0,tStart,length=2001))
     for i=1:size(TStepDist,2)-burnin
         TStepDistResized[:,i] = linterp1s(tSteps,TStepDist[:,i+burnin],xq)
     end
 
     # Calculate composite image
-    tTimage = zeros(ceil(Int, TCryst)*2, size(TStepDistResized,1))
-    yq = collect(0:0.5:TCryst)
+    tTimage = zeros(ceil(Int, TStart)*2, size(TStepDistResized,1))
+    yq = collect(0:0.5:TStart)
     for i=1:size(TStepDistResized,1)
         hist = fit(Histogram,TStepDistResized[i,:],yq,closed=:right)
         tTimage[:,i] = hist.weights
diff --git a/src/Thermochron.jl b/src/Thermochron.jl
index d4e1245..325ca0f 100644
--- a/src/Thermochron.jl
+++ b/src/Thermochron.jl
@@ -5,6 +5,7 @@ module Thermochron
     using LoopVectorization
     using ProgressMeter: @showprogress
 
+    include("minerals.jl")
     include("ZrnHe.jl")
     include("inversion.jl")
 
diff --git a/src/ZrnHe.jl b/src/ZrnHe.jl
index a11112a..a5434ed 100644
--- a/src/ZrnHe.jl
+++ b/src/ZrnHe.jl
@@ -1,67 +1,4 @@
 
-"""
-```julia
-intersectionfraction(r₁, r₂, d)
-```
-Calculate the fraction of the surface area of a sphere s2 with radius `r₂`
-that intersects the interior of a sphere s1 of radius `r₁` if the two are
-separated by distance `d`. Assumptions: `r₁`>0, `r₂`>0, `d`>=0
-"""
-function intersectionfraction(r₁::T, r₂::T, d::T) where T <: Number
-    dmax = r₁+r₂
-    dmin = abs(r₁-r₂)
-    if d > dmax ## If separated by more than dmax, there is no intersection
-        omega = zero(T)
-    elseif d > dmin
-        # X is the radial distance between the center of s2 and the interseciton plane
-        # See e.g., http://mathworld.wolfram.com/Sphere-SphereIntersection.html
-        # x = (d^2 - r₁^2 + r₂^2) / (2 * d)
-
-        # Let omega be is the solid angle of intersection normalized by 4pi,
-        # where the solid angle of a cone is 2pi*(1-cos(theta)) and cos(theta)
-        # is adjacent/hypotenuse = x/r₂.
-        # omega = (1 - x/r₂)/2
-
-        # Rearranged for optimization:
-        omega = one(T)/2 - (d^2 - r₁^2 + r₂^2) / (4 * d * r₂)
-
-    elseif r₁<r₂ # If r₁ is entirely within r₂
-        omega = zero(T)
-    else
-        omega = one(T) # If r₂ is entirely within r₁
-    end
-    return omega
-end
-export intersectionfraction
-
-
-"""
-```julia
-intersectiondensity(rEdges::Vector, rVolumes::Vector, ralpha, d)
-```
-Calculate the volume-nomalized fractional intersection density of an alpha
-stopping sphere of radius `ralpha` with each concentric shell (with shell edges
-`rEdges` and relative volumes `rVolumes`) of a spherical crystal where the
-two are separated by distance `d`
-"""
-function intersectiondensity(rEdges::Vector{T}, rVolumes::Vector{T}, ralpha::T, d::T) where T <: Number
-    n = length(rEdges) - 1
-    dInt = Array{Float64,1}(undef, n)
-    fIntLast = intersectionfraction(rEdges[1],ralpha,d)
-    for i=1:n
-        # Integrated intersection fraction for each concentric sphere (rEdges) of crystal
-        fInt = intersectionfraction(rEdges[i+1],ralpha,d)
-
-        # Intersection fraction for each spherical shell of the crystal (subtracting
-        # one concentric sphere from the next) normalized by shell volume
-        dInt[i] = (fInt-fIntLast) / rVolumes[i]
-        fIntLast = fInt
-    end
-    return dInt
-end
-export intersectiondensity
-
-
 """
 ```julia
 ρᵣ = DamageAnnealing(dt::Number, tSteps::Vector, TSteps::Matrix)
@@ -131,29 +68,20 @@ export DamageAnnealing!
 
 """
 ```julia
-HeAge = ZrnHeAgeSpherical(dt, ageSteps::Vector, TSteps::Vector, ρᵣ::Matrix, r, dr, Uppm, Thppm)
+HeAge = HeAgeSpherical(zircon::Zircon, TSteps::Vector, ρᵣ::Matrix, diffusionparams)
 ```
 Calculate the precdicted U-Th/He age of a zircon that has experienced a given t-T path
-(specified by `ageSteps` for time and `TSteps` for temperature, at a time resolution of `dt`)
-using a Crank-Nicholson diffusion solution for a spherical grain of radius `r` at spatial resolution `dr`.
+(specified by `zircon.ageSteps` for time and `TSteps` for temperature, at a time resolution of `zircon.dt`)
+using a Crank-Nicholson diffusion solution for a spherical grain of radius `zircon.r` at spatial resolution `zircon.dr`.
 """
-function ZrnHeAgeSpherical(dt::Number, ageSteps::AbstractVector{T}, TSteps::AbstractVector{T}, ρᵣ::AbstractMatrix{T}, r::T, dr::Number, Uppm::T, Thppm::T, diffusionparams) where T <: Number
-    # Temporal discretization
-    tSteps = reverse(ageSteps)
-    ntSteps = length(tSteps) # Number of time steps
+function HeAgeSpherical(zircon::Zircon{T}, TSteps::AbstractVector{T}, ρᵣ::AbstractMatrix{T}, diffusionparams) where T <: Number
 
     # Jaffey decay constants
     λ235U = log(2)/(7.0381*10^8)*10^6 # [1/Myr]
     λ238U = log(2)/(4.4683*10^9)*10^6 # [1/Myr]
     λ232Th = log(2)/(1.405*10^10)*10^6 # [1/Myr]
 
-    # Alpha stopping distances for each isotope in each decay chain, from
-    # Farley et al., 1996
-    alphaRadii238U = (11.78, 14.09, 13.73, 14.13, 17.32, 16.69, 28.56, 16.48,)
-    alphaRadii235U = (12.58, 15.04, 19.36, 18.06, 23.07, 26.87, 22.47,)
-    alphaRadii232Th = (10.99, 16.67, 18.16, 17.32, 23.61, 29.19,)
-
-    # Other constants
+    # Diffucsion constants
     # DzEa = 165.0 # kJ/mol
     # DzD0 = 193188.0 # cm^2/sec
     # DN17Ea = 71.0 # kJ/mol
@@ -163,7 +91,7 @@ function ZrnHeAgeSpherical(dt::Number, ageSteps::AbstractVector{T}, TSteps::Abst
     DN17Ea = diffusionparams.DN17Ea
     DN17D0 = diffusionparams.DN17D0
 
-
+    # Damage and annealing constants
     lint0 = 45920.0 # nm
     SV = 1.669 # 1/nm
     Bα = 5.48E-19 # [g/alpha] mass of amorphous material produced per alpha decay
@@ -176,77 +104,24 @@ function ZrnHeAgeSpherical(dt::Number, ageSteps::AbstractVector{T}, TSteps::Abst
     Dz = Dz*10000^2*(1E6*365.25*24*3600) # Convert to micron^2/Myr
     DN17 = DN17*10000^2*(1E6*365.25*24*3600) # Convert to micron^2/Myr
 
-    # Crystal size and spatial discretization
-    rSteps = Array{Float64}(0+dr/2 : dr: r-dr/2)
-    rEdges = Array{Float64}(0 : dr : r) # Edges of each radius element
-    nrSteps = length(rSteps)+2 # number of radial grid points
-    relVolumes = (rEdges[2:end].^3 - rEdges[1:end-1].^3)/rEdges[end]^3 # Relative volume fraction of spherical shell corresponding to each radius element
-
-
-    # Observed radial HPE profiles at present day
-    r238U = Uppm.*ones(size(rSteps)) # PPM
-    r235U = r238U/137.818.*ones(size(rSteps)) #PPM
-    r232Th = Thppm.*ones(size(rSteps)) #PPM
-
-    # Convert to atoms per gram
-    r238U *= 6.022E23 / 1E6 / 238
-    r235U *= 6.022E23 / 1E6 / 235
-    r232Th *= 6.022E23 / 1E6 / 232
-
-
-    # Calculate effective He deposition for each decay chain, corrected for alpha
-    # stopping distance
-
-    #238U
-    r238UHe = zeros(size(r238U))
-    for i=1:length(alphaRadii238U)
-        for ri = 1:length(rSteps)
-            # Effective radial alpha deposition from 238U
-            r238UHe += relVolumes[ri] * intersectiondensity(rEdges,relVolumes,alphaRadii238U[i],rSteps[ri]) * r238U[ri]
-        end
-    end
-
-    #235U
-    r235UHe = zeros(size(r235U))
-    for i=1:length(alphaRadii235U)
-        for ri = 1:length(rSteps)
-            # Effective radial alpha deposition from 235U
-            r235UHe += relVolumes[ri] * intersectiondensity(rEdges,relVolumes,alphaRadii235U[i],rSteps[ri]) * r235U[ri]
-        end
-    end
-
-    #232Th
-    r232ThHe = zeros(size(r232Th))
-    for i=1:length(alphaRadii232Th)
-        for ri = 1:length(rSteps)
-            # Effective radial alpha deposition from 232Th
-            r232ThHe += relVolumes[ri] * intersectiondensity(rEdges,relVolumes,alphaRadii232Th[i],rSteps[ri]) * r232Th[ri]
-        end
-    end
-
-    # Calculate corrected alpha deposition each time step for each radius
-    alphaDeposition = (exp.(λ238U.*(ageSteps .+ dt/2)) - exp.(λ238U.*(ageSteps .- dt/2))) * (r238UHe') +
-        (exp.(λ235U.*(ageSteps .+ dt/2)) - exp.(λ235U.*(ageSteps .- dt/2))) * (r235UHe') +
-        (exp.(λ232Th.*(ageSteps .+ dt/2)) - exp.(λ232Th.*(ageSteps .- dt/2))) * (r232ThHe')
-
-    r238Udam = 8*r238U # No smoothing of alpha damage, 8 alphas per 238U
-    r235Udam = 7*r235U # No smoothing of alpha damage, 7 alphas per 235U
-    r232Thdam = 6*r232Th # No smoothing of alpha damage, 6 alphas per 232 Th
-
-    # Calculate corrected alpha damage at each time step for each radius
-    alphaDamage = (exp.(λ238U.*(ageSteps .+ dt/2)) - exp.(λ238U.*(ageSteps .- dt/2))) * (r238Udam') +
-        (exp.(λ235U.*(ageSteps .+ dt/2)) - exp.(λ235U.*(ageSteps .- dt/2))) * (r235Udam') +
-        (exp.(λ232Th.*(ageSteps .+ dt/2)) - exp.(λ232Th.*(ageSteps .- dt/2))) * (r232Thdam')
-
+    # Get time and radius discretization
+    dr = zircon.dr
+    rSteps = zircon.rSteps
+    nrSteps = zircon.nrSteps
+    dt = zircon.dt
+    ntSteps = zircon.ntSteps
+    alphaDeposition = zircon.alphaDeposition
+    alphaDamage = zircon.alphaDamage
 
     # The annealed damage matrix is the summation of the ρᵣ for each
     # previous timestep multiplied by the the alpha dose at each
     # previous timestep; this is a linear combination, which can be
     # calculated efficiently for all radii by simple matrix multiplication.
-    annealedDamage = ρᵣ * alphaDamage
+    annealedDamage = zircon.annealedDamage
+    mul!(annealedDamage, ρᵣ, alphaDamage)
 
     # Calculate initial alpha damage
-    β = Array{Float64}(undef, nrSteps)
+    β = zircon.β
     @turbo for k = 1:(nrSteps-2)
         fₐ = 1-exp(-Bα*annealedDamage[1,k]*Phi)
         τ = (lint0/(4.2 / ((1-exp(-Bα*annealedDamage[1,k])) * SV) - 2.5))^2
@@ -256,31 +131,28 @@ function ZrnHeAgeSpherical(dt::Number, ageSteps::AbstractVector{T}, TSteps::Abst
     β[1] = β[2]
     β[end] = β[end-1]
 
-    # Allocate output matrix for all timesteps
+    # Output matrix for all timesteps
     # u = v*r is the coordinate transform (u-substitution) for the Crank-
     # Nicholson equations where v is the He profile and r is radius
-    u = Array{Float64}(undef, ntSteps, nrSteps)
+    u = zircon.u
     u[1,:] .= 0 # Initial u = v = 0 everywhere
 
-    # Declare variables for tridiagonal matrix and RHS
-    dl = ones(nrSteps-1)    # Sub-diagonal row
-    d = -2.0 .- β           # Diagonal
-    du = ones(nrSteps-1)    # Supra-diagonal row
-    du2 = ones(nrSteps-2)   # sup-sup-diagonal row for pivoting
+    # Vector for RHS of Crank-Nicholson equation with regular grid cells
+    y = zircon.y
+
+    # Tridiagonal matrix for LHS of Crank-Nicholson equation with regular grid cells
+    A = zircon.A
+    @. A.dl = 1         # Sub-diagonal row
+    @. A.d = -2 - β     # Diagonal
+    @. A.du = 1         # Supra-diagonal row
 
     # Neumann inner boundary condition (u(i,1) + u(i,2) = 0)
-    d[1] = 1
-    du[1] = 1
+    A.d[1] = 1
+    A.du[1] = 1
 
     # Dirichlet outer boundary condition (u(i,end) = u(i-1,end))
-    dl[nrSteps-1] = 0
-    d[nrSteps] = 1
-
-    # Tridiagonal matrix for LHS of Crank-Nicholson equation with regular grid cells
-    A = Tridiagonal(dl, d, du, du2)
-
-    # Vector for RHS of Crank-Nicholson equation with regular grid cells
-    y = Array{Float64}(undef, nrSteps)
+    A.dl[nrSteps-1] = 0
+    A.d[nrSteps] = 1
 
     @inbounds for i=2:ntSteps
 
@@ -324,22 +196,23 @@ function ZrnHeAgeSpherical(dt::Number, ageSteps::AbstractVector{T}, TSteps::Abst
 
     # Convert from u (coordinate-transform'd concentration) to v (real He
     # concentration)
-    vFinal = u[end,:]./[-dr/2; rSteps; r+dr/2]
-    HeAvg = mean(vFinal[2:end-1]) # Atoms/gram
+    vFinal = u[end,:]./[first(rSteps)-dr; rSteps; last(rSteps)+dr]
+    μHe = mean(vFinal[2:end-1]) # Atoms/gram
 
     # Raw Age (i.e., as measured)
-    Avg238U = mean(r238U) # Atoms/gram
-    Avg235U = mean(r235U)
-    Avg232Th = mean(r232Th)
-    He(t) = Avg238U*8*(exp(λ238U*t)-1) + Avg235U*7*(exp(λ235U*t)-1) + Avg232Th*6*(exp(λ232Th*t)-1)
+    μ238U = mean(zircon.r238U) # Atoms/gram
+    μ235U = mean(zircon.r235U)
+    μ232Th = mean(zircon.r232Th)
+
+    He(t) = μ238U*8*(exp(λ238U*t)-1) + μ235U*7*(exp(λ235U*t)-1) + μ232Th*6*(exp(λ232Th*t)-1)
 
     # Numerically solve for helium age of the grain
     HeAge = 1
     for i=1:10
         dHe = 10000*(He(HeAge+1/10000)-He(HeAge)) # Calculate derivative
-        HeAge += (HeAvg-He(HeAge))/dHe # Move towards zero (He(HeAge) == HeAvg)
+        HeAge += (μHe-He(HeAge))/dHe # Move towards zero (He(HeAge) == μHe)
     end
 
     return HeAge
 end
-export ZrnHeAgeSpherical
+export HeAgeSpherical
diff --git a/src/minerals.jl b/src/minerals.jl
new file mode 100644
index 0000000..3cf6aa7
--- /dev/null
+++ b/src/minerals.jl
@@ -0,0 +1,211 @@
+"""
+```julia
+intersectionfraction(r₁, r₂, d)
+```
+Calculate the fraction of the surface area of a sphere s2 with radius `r₂`
+that intersects the interior of a sphere s1 of radius `r₁` if the two are
+separated by distance `d`. Assumptions: `r₁`>0, `r₂`>0, `d`>=0
+"""
+function intersectionfraction(r₁::T, r₂::T, d::T) where T <: Number
+    dmax = r₁+r₂
+    dmin = abs(r₁-r₂)
+    if d > dmax ## If separated by more than dmax, there is no intersection
+        omega = zero(T)
+    elseif d > dmin
+        # X is the radial distance between the center of s2 and the interseciton plane
+        # See e.g., http://mathworld.wolfram.com/Sphere-SphereIntersection.html
+        # x = (d^2 - r₁^2 + r₂^2) / (2 * d)
+
+        # Let omega be is the solid angle of intersection normalized by 4pi,
+        # where the solid angle of a cone is 2pi*(1-cos(theta)) and cos(theta)
+        # is adjacent/hypotenuse = x/r₂.
+        # omega = (1 - x/r₂)/2
+
+        # Rearranged for optimization:
+        omega = one(T)/2 - (d^2 - r₁^2 + r₂^2) / (4 * d * r₂)
+
+    elseif r₁<r₂ # If r₁ is entirely within r₂
+        omega = zero(T)
+    else
+        omega = one(T) # If r₂ is entirely within r₁
+    end
+    return omega
+end
+export intersectionfraction
+
+
+"""
+```julia
+intersectiondensity(rEdges::Vector, rVolumes::Vector, ralpha, d)
+```
+Calculate the volume-nomalized fractional intersection density of an alpha
+stopping sphere of radius `ralpha` with each concentric shell (with shell edges
+`rEdges` and relative volumes `rVolumes`) of a spherical crystal where the
+two are separated by distance `d`
+"""
+function intersectiondensity(rEdges::Vector{T}, rVolumes::Vector{T}, ralpha::T, d::T) where T <: Number
+    n = length(rEdges) - 1
+    dInt = Array{Float64,1}(undef, n)
+    fIntLast = intersectionfraction(rEdges[1],ralpha,d)
+    for i=1:n
+        # Integrated intersection fraction for each concentric sphere (rEdges) of crystal
+        fInt = intersectionfraction(rEdges[i+1],ralpha,d)
+
+        # Intersection fraction for each spherical shell of the crystal (subtracting
+        # one concentric sphere from the next) normalized by shell volume
+        dInt[i] = (fInt-fIntLast) / rVolumes[i]
+        fIntLast = fInt
+    end
+    return dInt
+end
+export intersectiondensity
+
+abstract type Mineral end
+
+struct Zircon{T<:Number} <: Mineral
+    dt::T
+    ageSteps::Vector{T}
+    tSteps::Vector{T}
+    ntSteps::Int
+    dr::T
+    rSteps::Vector{T}
+    rEdges::Vector{T}
+    relVolumes::Vector{T}
+    nrSteps::Int
+    r238U::Vector{T}
+    r235U::Vector{T}
+    r232Th::Vector{T}
+    r238UHe::Vector{T}
+    r235UHe::Vector{T}
+    r232ThHe::Vector{T}
+    alphaDeposition::Matrix{T}
+    alphaDamage::Matrix{T}
+    annealedDamage::Matrix{T}
+    u::Matrix{T}
+    β::Vector{T}
+    A::Tridiagonal{T, Vector{T}}
+    y::Vector{T}
+end
+
+function Zircon(r::T, dr::Number, Uppm::T, Thppm::T, dt::Number, ageSteps::AbstractVector{T}) where T<:Number
+
+    # Temporal discretization
+    tSteps = reverse(ageSteps)
+    ntSteps = length(tSteps) # Number of time steps
+
+    # Crystal size and spatial discretization
+    rSteps = Array{T}(0+dr/2 : dr: r-dr/2)
+    rEdges = Array{T}(0 : dr : r) # Edges of each radius element
+    nrSteps = length(rSteps)+2 # number of radial grid points
+    relVolumes = (rEdges[2:end].^3 - rEdges[1:end-1].^3)/rEdges[end]^3 # Relative volume fraction of spherical shell corresponding to each radius element
+
+    # Alpha stopping distances for each isotope in each decay chain, from
+    # Farley et al., 1996
+    alphaRadii238U = (11.78, 14.09, 13.73, 14.13, 17.32, 16.69, 28.56, 16.48,)
+    alphaRadii235U = (12.58, 15.04, 19.36, 18.06, 23.07, 26.87, 22.47,)
+    alphaRadii232Th = (10.99, 16.67, 18.16, 17.32, 23.61, 29.19,)
+
+    # Jaffey decay constants
+    λ235U = log(2)/(7.0381*10^8)*10^6 # [1/Myr]
+    λ238U = log(2)/(4.4683*10^9)*10^6 # [1/Myr]
+    λ232Th = log(2)/(1.405*10^10)*10^6 # [1/Myr]
+
+    # Observed radial HPE profiles at present day
+    r238U = Uppm.*ones(T, size(rSteps)) # PPM
+    r235U = T.(r238U/137.818) #PPM
+    r232Th = Thppm .* ones(T, size(rSteps)) #PPM
+
+    # Convert to atoms per gram
+    r238U *= 6.022E23 / 1E6 / 238
+    r235U *= 6.022E23 / 1E6 / 235
+    r232Th *= 6.022E23 / 1E6 / 232
+
+    # Calculate effective He deposition for each decay chain, corrected for alpha
+    # stopping distance
+
+    #238U
+    r238UHe = zeros(size(r238U))
+    for i=1:length(alphaRadii238U)
+        for ri = 1:length(rSteps)
+            # Effective radial alpha deposition from 238U
+            r238UHe += relVolumes[ri] * intersectiondensity(rEdges,relVolumes,alphaRadii238U[i],rSteps[ri]) * r238U[ri]
+        end
+    end
+
+    #235U
+    r235UHe = zeros(size(r235U))
+    for i=1:length(alphaRadii235U)
+        for ri = 1:length(rSteps)
+            # Effective radial alpha deposition from 235U
+            r235UHe += relVolumes[ri] * intersectiondensity(rEdges,relVolumes,alphaRadii235U[i],rSteps[ri]) * r235U[ri]
+        end
+    end
+
+    #232Th
+    r232ThHe = zeros(size(r232Th))
+    for i=1:length(alphaRadii232Th)
+        for ri = 1:length(rSteps)
+            # Effective radial alpha deposition from 232Th
+            r232ThHe += relVolumes[ri] * intersectiondensity(rEdges,relVolumes,alphaRadii232Th[i],rSteps[ri]) * r232Th[ri]
+        end
+    end
+
+    # Calculate corrected alpha deposition each time step for each radius
+    alphaDeposition = (exp.(λ238U.*(ageSteps .+ dt/2)) - exp.(λ238U.*(ageSteps .- dt/2))) * (r238UHe') +
+        (exp.(λ235U.*(ageSteps .+ dt/2)) - exp.(λ235U.*(ageSteps .- dt/2))) * (r235UHe') +
+        (exp.(λ232Th.*(ageSteps .+ dt/2)) - exp.(λ232Th.*(ageSteps .- dt/2))) * (r232ThHe')
+
+    r238Udam = 8*r238U # No smoothing of alpha damage, 8 alphas per 238U
+    r235Udam = 7*r235U # No smoothing of alpha damage, 7 alphas per 235U
+    r232Thdam = 6*r232Th # No smoothing of alpha damage, 6 alphas per 232 Th
+
+    # Calculate corrected alpha damage at each time step for each radius
+    alphaDamage = (exp.(λ238U.*(ageSteps .+ dt/2)) - exp.(λ238U.*(ageSteps .- dt/2))) * (r238Udam') +
+        (exp.(λ235U.*(ageSteps .+ dt/2)) - exp.(λ235U.*(ageSteps .- dt/2))) * (r235Udam') +
+        (exp.(λ232Th.*(ageSteps .+ dt/2)) - exp.(λ232Th.*(ageSteps .- dt/2))) * (r232Thdam')
+
+    # Allocate additional variables that will be needed for Crank-Nicholson
+    annealedDamage = similar(alphaDamage)
+    β = Array{T}(undef, nrSteps) # First row of annealedDamage
+
+    # Allocate output matrix for all timesteps
+    u = Array{T}(undef, ntSteps, nrSteps)
+
+    # Allocate variables for tridiagonal matrix and RHS
+    dl = ones(T, nrSteps-1)    # Sub-diagonal row
+    d = ones(T, nrSteps)       # Diagonal
+    du = ones(T, nrSteps-1)    # Supra-diagonal row
+    du2 = ones(T, nrSteps-2)   # sup-sup-diagonal row for pivoting
+
+    # Tridiagonal matrix for LHS of Crank-Nicholson equation with regular grid cells
+    A = Tridiagonal(dl, d, du, du2)
+
+    # Vector for RHS of Crank-Nicholson equation with regular grid cells
+    y = Array{T}(undef, nrSteps)
+
+    return Zircon{T}(
+        dt,
+        ageSteps,
+        tSteps,
+        ntSteps,
+        dr,
+        rSteps,
+        rEdges,
+        relVolumes,
+        nrSteps,
+        r238U,
+        r235U,
+        r232Th,
+        r238UHe,
+        r235UHe,
+        r232ThHe,
+        alphaDeposition,
+        alphaDamage,
+        annealedDamage,
+        u,
+        β,
+        A,
+        y,
+    )
+end
+export Zircon
diff --git a/test/runtests.jl b/test/runtests.jl
index b5a4291..2bc3995 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -2,6 +2,6 @@ using Thermochron
 using LinearAlgebra
 using Test
 
-@testset "ZrnHe.jl" begin include("testZrnHe.jl") end
-@testset "inversion.jl" begin include("testinversion.jl") end
-@testset "Examples" begin include("testexamples.jl") end
+@testset "Zircon helium" begin include("testZrnHe.jl") end
+@testset "Inversion" begin include("testinversion.jl") end
+@testset "Integrated Examples" begin include("testexamples.jl") end
diff --git a/test/testZrnHe.jl b/test/testZrnHe.jl
index f35fd4b..50b6781 100644
--- a/test/testZrnHe.jl
+++ b/test/testZrnHe.jl
@@ -60,9 +60,11 @@
     crystalRadius = 29.26
     Uppm = 462.98
     Thppm = 177.76
-    @test round(ZrnHeAgeSpherical(dt,reverse(tSteps),TSteps,pr,crystalRadius,dr,Uppm,Thppm,diffusionparams), sigdigits=5) ≈ 520.03
+    zircon = Zircon(crystalRadius,dr,Uppm,Thppm,dt,reverse(tSteps))
+    @test round(HeAgeSpherical(zircon,TSteps,pr,diffusionparams), sigdigits=5) ≈ 520.03
 
     crystalRadius = 35.
     Uppm = 1107.
     Thppm = 351.
-    @test round(ZrnHeAgeSpherical(dt,reverse(tSteps),TSteps,pr,crystalRadius,dr,Uppm,Thppm,diffusionparams), sigdigits=5) ≈ 309.76
+    zircon = Zircon(crystalRadius,dr,Uppm,Thppm,dt,reverse(tSteps))
+    @test round(HeAgeSpherical(zircon,TSteps,pr,diffusionparams), sigdigits=5) ≈ 309.76
diff --git a/test/testexamples.jl b/test/testexamples.jl
index 1a44969..ab4313c 100644
--- a/test/testexamples.jl
+++ b/test/testexamples.jl
@@ -21,8 +21,8 @@ dt = 10 # time step size in Myr
 dTmax = 10.0 # Maximum reheating/burial per model timestep (to prevent numerical overflow)
 
 # Other model parameters
-CrystAgeMax_Ma = 4000.0 # Ma -- forbid anything older than this
-TCryst = 400.0 # Temperature (in C)
+CrystAgeMax = 4000.0 # Ma -- forbid anything older than this
+TStart = 400.0 # Temperature (in C)
 dr = 1 # Radius step, in microns
 
 diffusionparams = (;
@@ -53,19 +53,19 @@ unconf_TPoints = Float64[]
 
 # Populate local variables from data frame with specified options
 Halfwidth = data.Halfwidth_um
-U_ppm = data.U238_ppm
-Th_ppm = data.Th232_ppm
-HeAge_Ma = data.HeAge_Ma_raw
-HeAge_Ma_sigma = data.HeAge_Ma_sigma_raw
-CrystAge_Ma = data.CrystAge_Ma
+Uppm = data.U238_ppm
+Thppm = data.Th232_ppm
+HeAge = data.HeAge_Ma_raw
+HeAge_sigma = data.HeAge_Ma_sigma_raw
+CrystAge = data.CrystAge_Ma
 
-CrystAge_Ma[CrystAge_Ma .> CrystAgeMax_Ma] .= CrystAgeMax_Ma
-tCryst = ceil(maximum(CrystAge_Ma)/dt) * dt
+CrystAge[CrystAge .> CrystAgeMax] .= CrystAgeMax
+tStart = ceil(maximum(CrystAge)/dt) * dt
 
-tSteps = Array{Float64,1}(0+dt/2 : dt : tCryst-dt/2)
+tSteps = Array{Float64,1}(0+dt/2 : dt : tStart-dt/2)
 ageSteps = reverse(tSteps)
 ntSteps = length(tSteps) # Number of time steps
-eU = U_ppm+.238*Th_ppm # Used only for plotting
+eU = Uppm+.238*Thppm # Used only for plotting
 
 
 ## --- Invert for maximum likelihood t-T path
@@ -73,8 +73,8 @@ eU = U_ppm+.238*Th_ppm # Used only for plotting
 # Boundary conditions (10C at present and 650 C at the time of zircon
 # formation). Optional: Apppend required unconformities at base of Cambrian,
 # or elsewhere
-boundary_agePoints = Float64[0, tCryst] # Ma
-boundary_TPoints = Float64[0, TCryst] # Degrees C
+boundary_agePoints = Float64[0, tStart] # Ma
+boundary_TPoints = Float64[0, TStart] # Degrees C
 
 # This is where the "transdimensional" part comes in
 nPoints = 0
@@ -84,7 +84,7 @@ TPoints = Array{Float64}(undef, maxPoints+1) # Array of fixed size to hold all o
 
 # Fill some intermediate points to give the MCMC something to work with
 Tr = 250 # Residence temperature of initial proposal (value should not matter too much)
-for t in [tCryst/30,tCryst/4,tCryst/2,tCryst-tCryst/4,tCryst-tCryst/30]
+for t in [tStart/30,tStart/4,tStart/2,tStart-tStart/4,tStart-tStart/30]
     global nPoints += 1
     agePoints[nPoints] = t # Ma
     TPoints[nPoints] = Tr # Degrees C
@@ -93,7 +93,7 @@ end
 # # (Optional) Start with something close to the expected path
 #Tr = 150
 #T0 = 30
-#agePoints[1:4] = Float64[tCryst*29/30,   720, 510, 250]) # Age (Ma)
+#agePoints[1:4] = Float64[tStart*29/30,   720, 510, 250]) # Age (Ma)
 #TPoints[1:4]  =  Float64[       Tr+T0, Tr+T0,  T0,  70]) # Temp. (C)
 #nPoints+=4
 
@@ -101,18 +101,20 @@ function MCMC_vartcryst(nPoints, maxPoints, agePoints, TPoints, unconf_agePoints
     # Calculate model ages for initial proposal
     TSteps = linterp1s([view(agePoints, 1:nPoints) ; boundary_agePoints ; unconf_agePoints],
                        [view(TPoints, 1:nPoints) ; boundary_TPoints ; unconf_TPoints], ageSteps)
-    CalcHeAges = Array{Float64}(undef, size(HeAge_Ma))
+    CalcHeAges = Array{Float64}(undef, size(HeAge))
     pr = DamageAnnealing(dt,tSteps,TSteps) # Damage annealing history
-    for i=1:length(Halfwidth)
+    zircons = Array{Zircon{Float64}}(undef, length(Halfwidth))
+    for i=1:length(zircons)
         # Iterate through each grain, calculate the modeled age for each
-        first_index = 1 + floor(Int64,(tCryst - CrystAge_Ma[i])/dt)
-        CalcHeAges[i] = ZrnHeAgeSpherical(dt,ageSteps[first_index:end],TSteps[first_index:end],pr[first_index:end,first_index:end],Halfwidth[i],dr,U_ppm[i],Th_ppm[i], diffusionparams)
+        first_index = 1 + floor(Int64,(tStart - CrystAge[i])/dt)
+        zircons[i] = Zircon(Halfwidth[i], dr, Uppm[i], Thppm[i], dt, ageSteps[first_index:end])
+        CalcHeAges[i] = HeAgeSpherical(zircons[i], @views(TSteps[first_index:end]), @views(pr[first_index:end,first_index:end]), diffusionparams)
     end
-    AnnealedSigma = simannealsigma.(1, HeAge_Ma_sigma; params=simannealparams)
-    UnAnnealedSigma = simannealsigma.(nsteps, HeAge_Ma_sigma; params=simannealparams)
+    AnnealedSigma = simannealsigma.(1, HeAge_sigma; params=simannealparams)
+    UnAnnealedSigma = simannealsigma.(nsteps, HeAge_sigma; params=simannealparams)
 
     # Log-likelihood for initial proposal
-    ll = normpdf_ll(HeAge_Ma, AnnealedSigma, CalcHeAges)
+    ll = normpdf_ll(HeAge, AnnealedSigma, CalcHeAges)
     if simplified
         ll -= log(nPoints)
     end
@@ -129,15 +131,15 @@ function MCMC_vartcryst(nPoints, maxPoints, agePoints, TPoints, unconf_agePoints
     CalcHeAgesₚ = similar(CalcHeAges)
 
     # Distributions to populate
-    HeAgeDist = Array{Float64}(undef, length(HeAge_Ma), nsteps)
+    HeAgeDist = Array{Float64}(undef, length(HeAge), nsteps)
     TStepDist = Array{Float64}(undef, ntSteps, nsteps)
     llDist = Array{Float64}(undef, nsteps)
     nDist = zeros(Int, nsteps)
     acceptanceDist = zeros(Bool, nsteps)
 
     # Standard deviations of Gaussian proposal distributions for temperature and time
-    t_sigma = tCryst/60
-    T_sigma = TCryst/60
+    t_sigma = tStart/60
+    T_sigma = TStart/60
 
     # Proposal probabilities (must sum to 1)
     move = 0.64
@@ -168,17 +170,17 @@ function MCMC_vartcryst(nPoints, maxPoints, agePoints, TPoints, unconf_agePoints
                 if agePointsₚ[k] < dt
                     # Don't let any point get too close to 0
                     agePointsₚ[k] += (dt - agePointsₚ[k])
-                elseif agePointsₚ[k] > (tCryst - dt)
-                    # Don't let any point get too close to tCryst
-                    agePointsₚ[k] -= (agePointsₚ[k] - (tCryst - dt))
+                elseif agePointsₚ[k] > (tStart - dt)
+                    # Don't let any point get too close to tStart
+                    agePointsₚ[k] -= (agePointsₚ[k] - (tStart - dt))
                 end
                 # Move the Temperature of one model point
                 if TPointsₚ[k] < 0
                     # Don't allow T<0
                     TPointsₚ[k] = 0
-                elseif TPointsₚ[k] > TCryst
-                    # Don't allow T>TCryst
-                    TPointsₚ[k] = TCryst
+                elseif TPointsₚ[k] > TStart
+                    # Don't allow T>TStart
+                    TPointsₚ[k] = TStart
                 end
 
                 # Interpolate proposed t-T path
@@ -196,8 +198,8 @@ function MCMC_vartcryst(nPoints, maxPoints, agePoints, TPoints, unconf_agePoints
             # Birth: add a new model point
             nPointsₚ += 1
             for i=1:maxattempts # Try maxattempts times to satisfy the reheating rate limit
-                agePointsₚ[nPointsₚ] = rand()*tCryst
-                TPointsₚ[nPointsₚ] = rand()*TCryst
+                agePointsₚ[nPointsₚ] = rand()*tStart
+                TPointsₚ[nPointsₚ] = rand()*TStart
 
                 # Interpolate proposed t-T path
                 TStepsₚ = linterp1s([view(agePointsₚ, 1:nPointsₚ) ; boundary_agePoints ; unconf_agePointsₚ],
@@ -229,8 +231,8 @@ function MCMC_vartcryst(nPoints, maxPoints, agePoints, TPoints, unconf_agePoints
                     # Allow the present temperature to vary from 0 to 10 degrees C
                     boundary_TPointsₚ[1] = 0+rand()*10
                 else
-                    # Allow the initial temperature to vary from TCryst to TCryst-50 C
-                    boundary_TPointsₚ[2] = TCryst-rand()*50
+                    # Allow the initial temperature to vary from TStart to TStart-50 C
+                    boundary_TPointsₚ[2] = TStart-rand()*50
                 end
                 if length(unconf_agePointsₚ) > 0
                     # If there's an imposed unconformity, adjust within parameters
@@ -258,15 +260,15 @@ function MCMC_vartcryst(nPoints, maxPoints, agePoints, TPoints, unconf_agePoints
 
          # Calculate model ages for each grain
         DamageAnnealing!(pr, dt, tSteps, TStepsₚ)
-        for i=1:length(Halfwidth)
-            first_index = 1 + floor(Int64,(tCryst - CrystAge_Ma[i])/dt)
-            @views CalcHeAgesₚ[i] = ZrnHeAgeSpherical(dt,ageSteps[first_index:end],TStepsₚ[first_index:end],pr[first_index:end,first_index:end],Halfwidth[i],dr,U_ppm[i],Th_ppm[i],diffusionparams)
+        for i=1:length(zircons)
+            first_index = 1 + floor(Int64,(tStart - CrystAge[i])/dt)
+            CalcHeAgesₚ[i] = HeAgeSpherical(zircons[i], @views(TStepsₚ[first_index:end]), @views(pr[first_index:end,first_index:end]), diffusionparams)
         end
 
         # Calculate log likelihood of proposal
-        AnnealedSigma .= simannealsigma.(n, HeAge_Ma_sigma; params=simannealparams)
-        llₚ = normpdf_ll(HeAge_Ma, AnnealedSigma, CalcHeAgesₚ)
-        llₗ = normpdf_ll(HeAge_Ma, AnnealedSigma, CalcHeAges) # Recalulate last one too with new AnnealedSigma
+        AnnealedSigma .= simannealsigma.(n, HeAge_sigma; params=simannealparams)
+        llₚ = normpdf_ll(HeAge, AnnealedSigma, CalcHeAgesₚ)
+        llₗ = normpdf_ll(HeAge, AnnealedSigma, CalcHeAges) # Recalulate last one too with new AnnealedSigma
         if simplified # slightly penalize more complex t-T paths
             llₚ -= log(nPointsₚ)
             llₗ -= log(nPoints)
@@ -292,7 +294,7 @@ function MCMC_vartcryst(nPoints, maxPoints, agePoints, TPoints, unconf_agePoints
         end
 
         # Record results for analysis and troubleshooting
-        llDist[n] = normpdf_ll(HeAge_Ma, UnAnnealedSigma, CalcHeAges) # Recalculated to constant baseline
+        llDist[n] = normpdf_ll(HeAge, UnAnnealedSigma, CalcHeAges) # Recalculated to constant baseline
         nDist[n] = nPoints # Distribution of # of points
         HeAgeDist[:,n] = CalcHeAges # Distribution of He ages