✨ Forward autodiff compatibility

Project.toml

@@ -12,14 +12,18 @@ Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
 ReferenceFrameRotations = "74f56ac7-18b3-5285-802d-d4bd4f104033"
 SatelliteToolboxBase = "9e17983a-0463-41a7-9a16-1682db6d8b66"
 SatelliteToolboxLegendre = "7fa26607-a272-47b2-9aa2-a6c1d419d9d2"
+SatelliteToolboxTransformations = "6b019ec1-7a1e-4f04-96c7-a9db1ca5514d"
 Scratch = "6c6a2e73-6563-6170-7368-637461726353"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 
 [compat]
 Crayons = "4"
 Dates = "1.6"
 DelimitedFiles = "1.6"
+DifferentiationInterface = "0.5"
 Downloads = "1"
+FiniteDiff = "2.24"
+ForwardDiff = "0.10"
 LinearAlgebra = "1.6"
 ReferenceFrameRotations = "3"
 SatelliteToolboxBase = "0.2, 0.3"
@@ -32,9 +36,12 @@ 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"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 SatelliteToolboxTransformations = "6b019ec1-7a1e-4f04-96c7-a9db1ca5514d"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test", "DelimitedFiles", "LinearAlgebra", "SatelliteToolboxTransformations"]
+test = ["Test", "DelimitedFiles", "DifferentiationInterface", "FiniteDiff", "ForwardDiff", "LinearAlgebra", "SatelliteToolboxTransformations"]

src/GravityModels/API.md

@@ -6,19 +6,22 @@ overloading the functions listed here.
 ## Structure
 
 All models require a structure with supertype `AbstractGravityModel{T<:Number}`, where `T`
-is the type of the coefficients in the model.
+is the type of the coefficients in the model. For the function that passes the seconds since
+JD_J2000, the return type `RT` is the promoted type of the set {`T`, `W`}.
 
 ## API Functions
 
 ```julia
 function coefficients(model::AbstractGravityModel{T}, degree::Int, order::Int, time::DataTime) where T<:Number -> T, T
+function coefficients(model::AbstractGravityModel{T}, degree::Int, order::Int, time::W) where {T<:Number, W<:Number}> -> RT, RT
 ```
 
 This function must return the coefficients `Clm` and `Slm` of the gravity `model` for the
 specified `degree`, `order`, and `time`. Hence:
 
 ```julia
 coefficients(model, 10, 8, DateTime("2023-06-19"))
+coefficients(model, 10, 8, 7.404048e8)
 ```
 
 must return a `Tuple{T, T}` with the `Clm` and `Slm`, respectively, for the degree 10, order

src/GravityModels/accelerations.jl

@@ -11,7 +11,7 @@
 ############################################################################################
 
 """
-    gravitational_acceleration(model::AbstractGravityModel{T}, r::AbstractVector, time::DateTime = DateTime("2000-01-01"); kwargs...) where T<:Number -> NTuple{3, T}
+    gravitational_acceleration(model::AbstractGravityModel{T}, r::AbstractVector{V}, time::W = -43200; kwargs...) where {T<:Number,V<:Number,W<:Number} -> NTuple{3, RT}
 
 Compute the gravitational acceleration [m / s²] represented in ITRF using the `model` in the
 position `r` [m], also represented in ITRF, at instant `time`. If the latter argument is
@@ -46,19 +46,19 @@ omitted, the J2000.0 epoch is used.
 
 # Returns
 
-- `T`: The derivative of the gravitational field w.r.t. the radius (`∂U/∂r`).
-- `T`: The derivative of the gravitational field w.r.t. the geocentric latitude (`∂U/∂ϕ`).
-- `T`: The derivative of the gravitational field w.r.t. the longitude (`∂U/∂λ`).
+- `RT`: The derivative of the gravitational field w.r.t. the radius (`∂U/∂r`).
+- `RT`: The derivative of the gravitational field w.r.t. the geocentric latitude (`∂U/∂ϕ`).
+- `RT`: The derivative of the gravitational field w.r.t. the longitude (`∂U/∂λ`).
 """
 function gravitational_acceleration(
     model::AbstractGravityModel{T},
-    r::AbstractVector,
-    time::DateTime = DateTime("2000-01-01");
-    max_degree::Number = -1,
-    max_order::Number = -1,
+    r::AbstractVector{V},
+    time::Number = -43200;
+    max_degree::Int = -1,
+    max_order::Int = -1,
     P::Union{Nothing, AbstractMatrix} = nothing,
     dP::Union{Nothing, AbstractMatrix} = nothing
-) where T<:Number
+) where {T<:Number,V<:Number}
 
     # Compute the partial derivatives of the gravitational field w.r.t. the spherical
     # coordinates.
@@ -107,7 +107,67 @@ function gravitational_acceleration(
 end
 
 """
-    gravity_acceleration(model::AbstractGravityModel{T}, r::AbstractVector, time::DateTime = DateTime("2000-01-01"); kwargs...) where T<:Number -> NTuple{3, T}
+    gravitational_acceleration(model::AbstractGravityModel{T}, r::AbstractVector{V}, time::DateTime = DateTime("2000-01-01"); kwargs...) where {T<:Number,V<:Number} -> NTuple{3, RT}
+
+Compute the gravitational acceleration [m / s²] represented in ITRF using the `model` in the
+position `r` [m], also represented in ITRF, at instant `time`. If the latter argument is
+omitted, the J2000.0 epoch is used.
+
+!!! note
+    Gravitational acceleration is the acceleration caused by the central body mass only,
+    i.e., without considering the centrifugal potential.
+
+# Keywords
+
+- `max_degree::Int`: Maximum degree used in the spherical harmonics when computing the
+    gravitational field derivative. If it is higher than the available number of
+    coefficients in the `model`, it will be clamped. If it is lower than 0, it will be set
+    to the maximum degree available. (**Default** = -1)
+- `max_order::Int`: Maximum order used in the spherical harmonics when computing the
+    gravitational field derivative. If it is higher than `max_degree`, it will be clamped.
+    If it is lower than 0, it will be set to the same value as `max_degree`.
+    (**Default** = -1)
+- `P::Union{Nothing, AbstractMatrix}`: An optional matrix that must contain at least
+    `max_degree + 1 × max_degree + 1` real numbers that will be used to store the Legendre
+    coefficients, reducing the allocations. If it is `nothing`, the matrix will be created
+    when calling the function.
+- `dP::Union{Nothing, AbstractMatrix}`: An optional matrix that must contain at least
+    `max_degree + 1 × max_degree + 1` real numbers that will be used to store the Legendre
+    derivative coefficients, reducing the allocations. If it is `nothing`, the matrix will
+    be created when calling the function.
+
+# Returns
+
+- `T`: The derivative of the gravitational field w.r.t. the radius (`∂U/∂r`).
+- `T`: The derivative of the gravitational field w.r.t. the geocentric latitude (`∂U/∂ϕ`).
+- `T`: The derivative of the gravitational field w.r.t. the longitude (`∂U/∂λ`).
+"""
+function gravitational_acceleration(
+    model::AbstractGravityModel{T},
+    r::AbstractVector{V},
+    time::DateTime = DateTime("2000-01-01");
+    max_degree::Int = -1,
+    max_order::Int = -1,
+    P::Union{Nothing, AbstractMatrix} = nothing,
+    dP::Union{Nothing, AbstractMatrix} = nothing
+) where {T<:Number,V<:Number}
+
+    time_JD = (datetime2julian(time) - JD_J2000) * 86400
+
+    return gravitational_acceleration(
+        model,
+        r,
+        time_JD;
+        max_degree = max_degree,
+        max_order = max_order,
+        P = P,
+        dP = dP,
+    )
+
+end
+
+"""
+    gravity_acceleration(model::AbstractGravityModel{T}, r::AbstractVector{V}, time::Number = -43200; kwargs...) where {T<:Number,V<:Number,W<:Number} -> NTuple{3, RT}
 
 Compute the gravity acceleration [m / s²] represented in ITRF using the `model` in the
 position `r` [m], also represented in ITRF, at instant `time`. If the latter argument is
@@ -147,13 +207,13 @@ omitted, the J2000.0 epoch is used.
 """
 function gravity_acceleration(
     model::AbstractGravityModel{T},
-    r::AbstractVector,
-    time::DateTime = DateTime("2000-01-01");
-    max_degree::Number = -1,
-    max_order::Number = -1,
+    r::AbstractVector{V},
+    time::Number = -43200;
+    max_degree::Int = -1,
+    max_order::Int = -1,
     P::Union{Nothing, AbstractMatrix} = nothing,
     dP::Union{Nothing, AbstractMatrix} = nothing
-) where T<:Number
+) where {T<:Number,V<:Number}
 
     # == Gravitational Acceleration ========================================================
 
@@ -224,3 +284,63 @@ function gravity_acceleration(
 
     return g_itrf
 end
+
+"""
+    gravity_acceleration(model::AbstractGravityModel{T}, r::AbstractVector{V}, time::DateTime = DateTime("2000-01-01"); kwargs...) where {T<:Number,V<:Number} -> NTuple{3, RT}
+
+Compute the gravity acceleration [m / s²] represented in ITRF using the `model` in the
+position `r` [m], also represented in ITRF, at instant `time`. If the latter argument is
+omitted, the J2000.0 epoch is used.
+
+!!! note
+    Gravity acceleration is the compound acceleration caused by the central body mass and
+    the centrifugal force due to the planet's rotation.
+
+# Keywords
+
+- `max_degree::Int`: Maximum degree used in the spherical harmonics when computing the
+    gravitational field derivative. If it is higher than the available number of
+    coefficients in the `model`, it will be clamped. If it is lower than 0, it will be set
+    to the maximum degree available. (**Default** = -1)
+- `max_order::Int`: Maximum order used in the spherical harmonics when computing the
+    gravitational field derivative. If it is higher than `max_degree`, it will be clamped.
+    If it is lower than 0, it will be set to the same value as `max_degree`.
+    (**Default** = -1)
+- `P::Union{Nothing, AbstractMatrix}`: An optional matrix that must contain at least
+    `max_degree + 1 × max_degree + 1` real numbers that will be used to store the Legendre
+    coefficients, reducing the allocations. If it is `nothing`, the matrix will be created
+    when calling the function.
+- `dP::Union{Nothing, AbstractMatrix}`: An optional matrix that must contain at least
+    `max_degree + 1 × max_degree + 1` real numbers that will be used to store the Legendre
+    derivative coefficients, reducing the allocations. If it is `nothing`, the matrix will
+    be created when calling the function.
+
+# Returns
+
+- `T`: The derivative of the gravitational field w.r.t. the radius (`∂U/∂r`).
+- `T`: The derivative of the gravitational field w.r.t. the geocentric latitude (`∂U/∂ϕ`).
+- `T`: The derivative of the gravitational field w.r.t. the longitude (`∂U/∂λ`).
+"""
+function gravity_acceleration(
+    model::AbstractGravityModel{T},
+    r::AbstractVector{V},
+    time::DateTime = DateTime("2000-01-01");
+    max_degree::Int = -1,
+    max_order::Int = -1,
+    P::Union{Nothing, AbstractMatrix} = nothing,
+    dP::Union{Nothing, AbstractMatrix} = nothing
+) where {T<:Number,V<:Number,W<:Number}
+
+    time_JD = datetime2julian(time) - JD_J2000
+
+    return gravity_acceleration(
+        model,
+        r,
+        time_JD;
+        max_degree = max_degree,
+        max_order = max_order,
+        P = P,
+        dP = dP,
+    )
+
+end

src/GravityModels/api.jl

@@ -5,15 +5,25 @@
 ############################################################################################
 
 """
-    coefficients(model::AbstractGravityModel{T}, degree::Int, order::Int, time::DateTime = DateTime("2000-01-01")) where T<:Number -> T, T
+    coefficients(model::AbstractGravityModel{T}, degree::Int, order::Int, time::Number = -43200) where T<:Number -> T, T
 
 Return the `Clm` and `Slm` coefficients of the gravity `model` for the specified `degree`,
 `order`, and `time`. If the latter argument is omitted, the J2000.0 epoch is used.
 """
 function coefficients end
 
-function coefficients(model::AbstractGravityModel, degree::Int, order::Int)
-    return coefficients(model, degree, order, DateTime("2000-01-01"))
+function coefficients(model::AbstractGravityModel, degree::Int, order::Int, time::Number = -43200)
+    return coefficients(model, degree, order, time)
+end
+
+"""
+    coefficients(model::AbstractGravityModel{T}, degree::Int, order::Int, time::DateTime = DateTime("2000-01-01")) where T<:Number -> T, T
+
+Return the `Clm` and `Slm` coefficients of the gravity `model` for the specified `degree`,
+`order`, and `time`. If the latter argument is omitted, the J2000.0 epoch is used.
+"""
+function coefficients(model::AbstractGravityModel, degree::Int, order::Int, time::DateTime = DateTime("2000-01-01"))
+    return coefficients(model, degree, order, time)
 end
 
 """