Make annealing models kwdef structs, update ketcham 1999 coeffs

OpenThermochronology · Dec 10, 2024 · b402166 · b402166
1 parent 543f6c3
commit b402166
Show file tree

Hide file tree

Showing 4 changed files with 53 additions and 54 deletions.
diff --git a/src/Thermochron.jl b/src/Thermochron.jl
@@ -26,6 +26,7 @@ module Thermochron
 
     include("types.jl")
     export ZRDAAM, RDAAM, FanningCurvilinear, SimplifiedCurvilinear     # Damage and annealing model types
+    export FCKetcham1999, FCKetcham2007                                 # Aliases
     export Constraint, Unconformity, Boundary, DetailInterval           # Types used as inputs to MCMC functions
 
     include("utilities.jl")

diff --git a/src/fissiontrack.jl b/src/fissiontrack.jl
@@ -3,12 +3,6 @@
 l₀ = 16.38 # [um] initial track length
 σl₀ = 0.09 # [um] initial track length uncertainty
 
-# "Simultaneous fit" Fanning Curvilinear models 
-const FCKetcham1999 = FanningCurvilinear(-19.84402202, 0.3895104539, -51.25312954, -7.642358713, -0.12327, -11.988, 16.38, 0.09)
-export FCKetcham1999
-
-const FCKetcham2007 = SimplifiedCurvilinear(0.39528, 0.01073, -65.12969, -7.91715, 0.04672, 16.38, 0.09)
-export FCKetcham2007
 
 function equivalenttime(t::Number, T::Number, Teq::Number, fc::Union{SimplifiedCurvilinear,FanningCurvilinear})
     exp(fc.C2 + (log(t*SEC_MYR)-fc.C2)*(log(1/(Teq+273.15))-fc.C3)/(log(1/(T+273.15))-fc.C3))/SEC_MYR

diff --git a/src/types.jl b/src/types.jl
@@ -7,26 +7,30 @@ Base.length(x::AnnealingModel) = 1
 Base.iterate(x::AnnealingModel) = (x, nothing)
 Base.iterate(x::AnnealingModel, state) = nothing
 
-struct FanningCurvilinear{T<:AbstractFloat} <: AnnealingModel{T} 
-    C0::T
-    C1::T
-    C2::T
-    C3::T
-    alpha::T
-    beta::T
-    l0::T
-    l0_sigma::T
+# Fanning Curvilinear models (Ketcham 1999)
+Base. @kwdef struct FanningCurvilinear{T<:AbstractFloat} <: AnnealingModel{T} 
+    C0::T = -19.844     # "Simultaneous fit" from Ketcham et al. 1999
+    C1::T = 0.38951     # "Simultaneous fit" from Ketcham et al. 1999
+    C2::T = -51.253     # "Simultaneous fit" from Ketcham et al. 1999
+    C3::T = -7.6423     # "Simultaneous fit" from Ketcham et al. 1999
+    alpha::T = -0.12327 # "Simultaneous fit" from Ketcham et al. 1999
+    beta::T = -11.988   # "Simultaneous fit" from Ketcham et al. 1999
+    l0::T = 16.38       # Initial track length
+    l0_sigma::T = 0.09  # Initial track length unertainty
 end
-
-struct SimplifiedCurvilinear{T<:AbstractFloat} <: AnnealingModel{T} 
-    C0::T
-    C1::T
-    C2::T
-    C3::T
-    alpha::T
-    l0::T
-    l0_sigma::T
+const FCKetcham1999 = FanningCurvilinear()
+
+# Simplified Fanning Curvilinear models (Ketcham 2007)
+Base.@kwdef struct SimplifiedCurvilinear{T<:AbstractFloat} <: AnnealingModel{T} 
+    C0::T = 0.39528     # "Simultaneous fit" from Ketcham et al. 2007
+    C1::T = 0.01073     # "Simultaneous fit" from Ketcham et al. 2007
+    C2::T = -65.12969   # "Simultaneous fit" from Ketcham et al. 2007
+    C3::T = -7.91715    # "Simultaneous fit" from Ketcham et al. 2007
+    alpha::T = 0.04672  # "Simultaneous fit" from Ketcham et al. 2007
+    l0::T = 16.38       # Initial track length
+    l0_sigma::T = 0.09  # Initial track length unertainty
 end
+const FCKetcham2007 = SimplifiedCurvilinear()
 
 ## --- Define DiffusivityModel types
 

diff --git a/test/testfissiontrack.jl b/test/testfissiontrack.jl
@@ -7,20 +7,20 @@
     am = FCKetcham1999
     @test am isa Thermochron.FanningCurvilinear{Float64}
 
-    @test Thermochron.reltracklength(1, 0, am) ≈ 0.9697793598061583
-    @test Thermochron.reltracklength(1, 10, am) ≈ 0.964698241959342
-    @test Thermochron.reltracklength(1, 100, am) ≈ 0.8757409060414316
-    @test Thermochron.reltracklength(1, 200, am) ≈ 0.6865992067636948
+    @test Thermochron.reltracklength(1, 0, am) ≈ 0.969772147771495
+    @test Thermochron.reltracklength(1, 10, am) ≈ 0.964689584499856
+    @test Thermochron.reltracklength(1, 100, am) ≈ 0.8757101462242018
+    @test Thermochron.reltracklength(1, 200, am) ≈ 0.6865397486321438
     @test Thermochron.reltracklength(1, 500, am) ≈ 0.0
 
     @test Thermochron.equivalenttime.(1:10, 100, 100, am) ≈ 1:10
-    @test Thermochron.equivalenttime.(1:10, 50, 100, am) ≈ [0.0017399827659966166, 0.003298712056416963, 0.004795641117812894, 0.0062537982811099796, 0.007683787046370734, 0.009091721758817036, 0.010481577045220772, 0.01185614029716061, 0.013217472138400533, 0.01456715633288802]
-    @test Thermochron.equivalenttime.(1:10, 150, 100, am) ≈ [660.2918270322567, 1394.7742722247667, 2160.135692291735, 2946.265864298046, 3748.209693768477, 4562.984978421849, 5388.588033345723, 6223.575180578797, 7066.854409034921, 7917.569525688446]
+    @test Thermochron.equivalenttime.(1:10, 50, 100, am) ≈ [0.0017396519912747071, 0.003298079407422488, 0.004794716652673411, 0.006252588352279524, 0.007682296289640212, 0.009089953816004663, 0.010479534909963182, 0.01185382650735338, 0.013214888894288724, 0.014564305574711722]
+    @test Thermochron.equivalenttime.(1:10, 150, 100, am) ≈  [660.4429511024003, 1395.0963084362782, 2160.6369846963353, 2946.952051746238, 3749.0850835490833, 4564.053074060813, 5389.851797105461, 6225.0371840101425, 7068.5169276849765, 7919.434602704827]
 
-    @test Thermochron.reltrackdensity(1, 0, am) ≈ 0.9516469756898532
-    @test Thermochron.reltrackdensity(1, 10, am) ≈ 0.9435171871349471
-    @test Thermochron.reltrackdensity(1, 100, am) ≈ 0.8011854496662908
-    @test Thermochron.reltrackdensity(1, 200, am) ≈ 0.3212180867210108
+    @test Thermochron.reltrackdensity(1, 0, am) ≈ 0.9516354364343921
+    @test Thermochron.reltrackdensity(1, 10, am) ≈ 0.9435033351997696 
+    @test Thermochron.reltrackdensity(1, 100, am) ≈ 0.801136233958723 
+    @test Thermochron.reltrackdensity(1, 200, am) ≈ 0.3210110092650815
     @test Thermochron.reltrackdensity(1, 500, am) ≈ 0.0
 
 ## --- Ketcham et al. 2007
@@ -64,33 +64,33 @@
     display(apatite)
 
     # Isothermal residence
-    @test modelage(apatite, fill(0, 100), FCKetcham1999) ≈ 89.48198019056863
+    @test modelage(apatite, fill(0, 100), FCKetcham1999) ≈ 89.47899236366592
     @test modelage(apatite, fill(0, 100), FCKetcham2007) ≈ 91.18888272552469
 
-    @test modelage(apatite, fill(50, 100), FCKetcham1999) ≈ 71.77101682177967
+    @test modelage(apatite, fill(50, 100), FCKetcham1999) ≈ 71.75968375010271
     @test modelage(apatite, fill(50, 100), FCKetcham2007) ≈ 74.20874464530966
 
-    @test modelage(apatite, fill(75, 100), FCKetcham1999) ≈ 22.14984699770526
+    @test modelage(apatite, fill(75, 100), FCKetcham1999) ≈ 22.088815691960797
     @test modelage(apatite, fill(75, 100), FCKetcham2007) ≈ 21.6047849747377
 
-    @test modelage(apatite, fill(100, 100), FCKetcham1999) ≈ 0.4704161148027164
+    @test modelage(apatite, fill(100, 100), FCKetcham1999) ≈ 0.4682983384208703 
     @test modelage(apatite, fill(100, 100), FCKetcham2007) ≈ 0.42338708872671615
 
     # Linear cooling
-    @test modelage(apatite, reverse(1:100), FCKetcham1999) ≈ 66.77980999111053
+    @test modelage(apatite, reverse(1:100), FCKetcham1999) ≈ 66.7612470244616
     @test modelage(apatite, reverse(1:100), FCKetcham2007) ≈ 68.4228509734102
 
     apatite = ApatiteFT(agesteps=reverse(cntr(0:200)), F=1.75, Cl=0.01, OH=0.24)
     @test apatite isa ApatiteFT{Float64}
     @test apatite.rmr0 ≈ 0.8573573076438294
 
-    @test modelage(apatite, reverse(1:200), FCKetcham1999) ≈ 65.22193612822629
+    @test modelage(apatite, reverse(1:200), FCKetcham1999) ≈ 65.20383965575716
     @test modelage(apatite, reverse(1:200), FCKetcham2007) ≈ 67.08221305145062
 
-    @test modelage(apatite, reverse(1:200)./2, FCKetcham1999) ≈ 125.34528518852127
+    @test modelage(apatite, reverse(1:200)./2, FCKetcham1999) ≈ 125.30907058717082
     @test modelage(apatite, reverse(1:200)./2, FCKetcham2007) ≈ 128.59190431998582
 
-    @test modelage(apatite, reverse(1:200).*2, FCKetcham1999) ≈ 33.39910745952449
+    @test modelage(apatite, reverse(1:200).*2, FCKetcham1999) ≈ 33.38974963378803
     @test modelage(apatite, reverse(1:200).*2, FCKetcham2007) ≈ 34.52279548202092
 
     apatite = ApatiteFT(age=25, age_sigma=3, agesteps=reverse(cntr(0:28)), dpar=2.16)
@@ -104,36 +104,36 @@
     display(track)
 
     l, σ = modellength(track, fill(75, 20), FCKetcham1999) .* 16.38
-    @test l ≈ 12.434055894969841
-    @test σ ≈  0.5731491659110551 
-    @test round.(track.r, sigdigits=4) ≈ [0.7103, 0.7135, 0.7169, 0.7203, 0.7239, 0.7276, 0.7315, 0.7355, 0.7398, 0.7442, 0.749, 0.754, 0.7595, 0.7654, 0.7719, 0.7793, 0.7878, 0.798, 0.8112, 0.8311]
+    @test l ≈ 12.432312726056672
+    @test σ ≈  0.5736749779212151
+    @test round.(track.r, sigdigits=4) ≈ [0.7101, 0.7134, 0.7167, 0.7202, 0.7238, 0.7275, 0.7314, 0.7354, 0.7396, 0.7441, 0.7488, 0.7539, 0.7593, 0.7653, 0.7718, 0.7792, 0.7877, 0.7979, 0.8111, 0.8311]
     l, σ = modellength(track, fill(75, 20), FCKetcham2007) .* 16.38
     @test l ≈ 12.610086362166317
     @test σ ≈  0.6062929557510826
     @test round.(track.r, sigdigits=4) ≈ [0.7176, 0.7212, 0.725, 0.7289, 0.7329, 0.737, 0.7413, 0.7457, 0.7503, 0.7552, 0.7603, 0.7657, 0.7716, 0.7779, 0.7848, 0.7925, 0.8014, 0.8119, 0.8254, 0.8454]
 
     l, σ = modellength(track, fill(100, 20), FCKetcham1999) .* 16.38
-    @test l ≈ 10.875830406686777
-    @test σ ≈  0.6360569439228593
-    @test round.(track.r, sigdigits=4) ≈ [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.4559, 0.5506, 0.6187, 0.6878]
+    @test l ≈ 10.874500968379554 
+    @test σ ≈  0.6344191502094181
+    @test round.(track.r, sigdigits=4) ≈ [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.4544, 0.55, 0.6184, 0.6876]
     l, σ = modellength(track, fill(100, 20), FCKetcham2007) .* 16.38
     @test l ≈ 10.982218818618923
     @test σ ≈  0.608187681868191
     @test round.(track.r, sigdigits=4) ≈ [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.4888, 0.6017, 0.6895]
 
     l, σ = modellength(track, reverse(1:20).*5, FCKetcham1999) .* 16.38
-    @test l ≈ 14.353530832531991
-    @test σ ≈  1.1780179817828795
-    @test round.(track.r, sigdigits=4) ≈ [0.6251, 0.6876, 0.7304, 0.7632, 0.7897, 0.8118, 0.8307, 0.8472, 0.8616, 0.8744, 0.8858, 0.896, 0.9052, 0.9135, 0.9211, 0.928, 0.9345, 0.9406, 0.9467, 0.9536]
+    @test l ≈ 14.353253239597725
+    @test σ ≈  1.1781557181942377
+    @test round.(track.r, sigdigits=4) ≈ [0.6248, 0.6874, 0.7303, 0.7631, 0.7896, 0.8117, 0.8307, 0.8471, 0.8615, 0.8743, 0.8857, 0.896, 0.9052, 0.9135, 0.9211, 0.928, 0.9344, 0.9405, 0.9466, 0.9536]
     l, σ = modellength(track, reverse(1:20).*5, FCKetcham2007) .* 16.38
     @test l ≈ 14.548247656072547
     @test σ ≈  1.1726032551206245
     @test round.(track.r, sigdigits=4) ≈ [0.6095, 0.6893, 0.7388, 0.7749, 0.8031, 0.8261, 0.8454, 0.8618, 0.8761, 0.8885, 0.8994, 0.9091, 0.9178, 0.9256, 0.9326, 0.9389, 0.9448, 0.9503, 0.9557, 0.9619]
 
     l, σ =  modellength(track, reverse(1:20).*10, FCKetcham1999) .* 16.38 
-    @test l ≈ 14.346205934457554
-    @test σ ≈  1.1850553007509892
-    @test round.(track.r, sigdigits=4) ≈[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.6595, 0.7476, 0.8012, 0.8392, 0.8682, 0.8911, 0.9095, 0.9247, 0.9375, 0.9494]
+    @test l ≈ 14.345912487894068 
+    @test σ ≈  1.1852462926118148
+    @test round.(track.r, sigdigits=4) ≈ [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.6592, 0.7475, 0.8011, 0.8392, 0.8682, 0.891, 0.9095, 0.9247, 0.9375, 0.9494]
     l, σ = modellength(track, reverse(1:20).*10, FCKetcham2007) .* 16.38
     @test l ≈ 14.529573753525176
     @test σ ≈  1.190954235529948