feat: add support and testing for ForwardDiff

Project.toml

@@ -34,7 +34,10 @@ julia = "1.6"
 
 [extras]
 DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
+DifferentiationInterface = "a0c0ee7d-e4b9-4e03-894e-1c5f64a51d63"
+FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41"
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test", "DelimitedFiles"]
+test = ["Test", "DelimitedFiles", "DifferentiationInterface", "FiniteDiff", "ForwardDiff"]

src/jb2008/jb2008.jl

@@ -117,20 +117,23 @@ function jb2008(
 end
 
 function jb2008(
-    jd::Number,
-    ϕ_gd::Number,
-    λ::Number,
-    h::Number,
-    F10::Number,
-    F10ₐ::Number,
-    S10::Number,
-    S10ₐ::Number,
-    M10::Number,
-    M10ₐ::Number,
-    Y10::Number,
-    Y10ₐ::Number,
-    DstΔTc::Number
-)
+    jd::JT,
+    ϕ_gd::PT,
+    λ::LT,
+    h::HT,
+    F10::FT,
+    F10ₐ::FT2,
+    S10::ST,
+    S10ₐ::ST2,
+    M10::MT,
+    M10ₐ::MT2,
+    Y10::YT,
+    Y10ₐ::YT2,
+    DstΔTc::DT
+) where {JT<:Number, PT<:Number, LT<:Number, HT<:Number, FT<:Number, FT2<:Number, ST<:Number, ST2<:Number, MT<:Number, MT2<:Number, YT<:Number, YT2<:Number, DT<:Number}
+
+    RT = promote_type(JT, PT, LT, HT, FT, FT2, ST, ST2, MT, MT2, YT, YT2, DT)
+
     ########################################################################################
     #                                      Constants                                       #
     ########################################################################################
@@ -502,7 +505,7 @@ function jb2008(
     ρ   = sum_mn / A
 
     # Create and return the output structure.
-    return JB2008Output{Float64}(
+    return JB2008Output{RT}(
         ρ,
         Tz,
         T_exo,

src/jr1971/jr1971.jl

@@ -95,14 +95,17 @@ function jr1971(
 end
 
 function jr1971(
-    jd::Number,
-    ϕ_gd::Number,
-    λ::Number,
-    h::Number,
-    F10::Number,
-    F10ₐ::Number,
-    Kp::Number
-)
+    jd::JT,
+    ϕ_gd::PT,
+    λ::LT,
+    h::HT,
+    F10::FT,
+    F10ₐ::FT2,
+    Kp::KT
+)   where {JT<:Number, PT<:Number, LT<:Number, HT<:Number, FT<:Number, FT2<:Number, KT<:Number}
+
+    RT = promote_type(JT, PT, LT, HT, FT, FT2, KT)
+
     # == Constants =========================================================================
 
     Rstar = _JR1971_CONSTANTS.Rstar
@@ -414,7 +417,7 @@ function jr1971(
             ρ = ρN₂ + ρO₂ + ρO + ρAr + ρHe
 
             # Convert to SI and return.
-            return JR1971Output(
+            return JR1971Output{RT}(
                 1000ρ,
                 Tz,
                 T∞,
@@ -495,7 +498,7 @@ function jr1971(
         ρ = ρN₂ + ρO₂ + ρO + ρAr + ρHe + ρH
 
         # Convert to SI and return.
-        return JR1971Output(
+        return JR1971Output{RT}(
             1000ρ,
             Tz,
             T∞,

src/nrlmsise00/auxiliary.jl

@@ -16,10 +16,13 @@ Compute the chemistry / dissociation correction for MSIS models.
 - `h₁::Number`: Transition scale length.
 - `zh::Number`: Altitude of `1/2 r`.
 """
-function _ccor(h::T, r::T, h₁::T, zh::T) where T<:Number
+function _ccor(h::HT, r::T, h₁::HT2, zh::ZT) where {HT<:Number, T<:Number, HT2<:Number, ZT<:Number}
+
+    RT = promote_type(HT, T, HT2, ZT)
+
     e = (h - zh) / h₁
 
-    (e > +70) && return exp(T(0))
+    (e > +70) && return exp(RT(0))
     (e < -70) && return exp(r)
 
     return exp(r / (1 + exp(e)))
@@ -38,7 +41,7 @@ Compute the O and O₂ chemistry / dissociation correction for MSIS models.
 - `zh::Number`: Altitude of `1/2 r`.
 - `h₂::Number`: Transition scale length 2.
 """
-function _ccor2(h::T, r::T, h₁::T, zh::T, h₂::T) where T<:Number
+function _ccor2(h::HT, r::T, h₁::HT2, zh::ZT, h₂::HT3) where {HT<:Number, T<:Number, HT2<:Number, ZT<:Number, HT3<:Number}
     e1 = (h - zh) / h₁
     e2 = (h - zh) / h₂
 
@@ -61,12 +64,15 @@ Compute the turbopause correction for MSIS models, returning the combined densit
 - `xmm::T`: Full mixed molecular weight.
 - `xm::T`: Species molecular weight.
 """
-function _dnet(dd::T, dm::T, zhm::T, xmm::T, xm::T) where T<:Number
+function _dnet(dd::DT, dm::DT2, zhm::ZT, xmm::XT, xm::XT2) where {DT<:Number, DT2<:Number, ZT<:Number, XT<:Number, XT2<:Number}
+
+    RT = promote_type(DT, DT2, ZT, XT, XT2)
+
     a  = zhm / (xmm - xm)
 
     if !((dm > 0) && (dd > 0))
         if (dd == 0) && (dm == 0)
-            dd = T(1)
+            dd = RT(1)
         end
 
         (dm == 0) && return dd
@@ -112,12 +118,15 @@ Compute the scale height.
 - `g_lat::T`: Reference gravity at desired latitude [cm / s²].
 - `r_lat::T`: Reference radius at desired latitude [km].
 """
-function _scale_height(h::T, xm::T, temp::T, g_lat::T, r_lat::T) where T<:Number
+function _scale_height(h::HT, xm::XT, temp::TT, g_lat::GT, r_lat::RLT) where {HT<:Number, XT<:Number, TT<:Number, GT<:Number, RLT<:Number}
+
+    RT = promote_type(HT, XT, TT, GT, RLT)
+
     # Compute the gravity at the selected altitude.
     gh = g_lat / (1 + h / r_lat)^2
 
     # Compute the scale height
-    sh = T(_RGAS) * temp / (gh * xm)
+    sh = RT(_RGAS) * temp / (gh * xm)
 
     return sh
 end

src/nrlmsise00/math.jl

@@ -173,10 +173,10 @@ function _spline(
 end
 
 """
-    _ζ(r_lat::T, zz::T, zl::T) where T<:Number -> float(T)
+    _ζ(r_lat::Number, zz::Number, zl::Number) -> Number
 
 Compute the zeta function.
 """
-function _ζ(r_lat::T, zz::T, zl::T) where T<:Number
+function _ζ(r_lat::Number, zz::Number, zl::Number)
     return (zz - zl) * (r_lat + zl) / (r_lat + zz)
 end