From a846d8c90871c88e32b07a71b1516cbcfa1ccd24 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jorge=20A=2E=20P=C3=A9rez-Hern=C3=A1ndez?= Date: Thu, 25 Jan 2024 23:06:30 +0100 Subject: [PATCH 01/10] WIP: try fix type parameters in TaylorInterpolant --- src/interpolation.jl | 144 +++++++++++++++++++++---------------------- 1 file changed, 72 insertions(+), 72 deletions(-) diff --git a/src/interpolation.jl b/src/interpolation.jl index 01c74d2..886a5d7 100644 --- a/src/interpolation.jl +++ b/src/interpolation.jl @@ -2,62 +2,62 @@ @doc raw""" TaylorInterpolant{T, U, N} - -Collection of Taylor polynomials that interpolate a dependent variable as a function of -an independent variable. For example, the ``x``-axis position of the Earth as a function of + +Collection of Taylor polynomials that interpolate a dependent variable as a function of +an independent variable. For example, the ``x``-axis position of the Earth as a function of time ``x(t)``; or a lunar core Euler angle as a function of time ``\theta_c(t)``. # Fields - `t0::T`: Start time. -- `t::AbstractVector{T}`: Vector of time instances when the timespan of the ``i``-th element of `x` ends and the ``(i+1)``-th element of `x` starts being valid. +- `t::AbstractVector{T}`: Vector of time instances when the timespan of the ``i``-th element of `x` ends and the ``(i+1)``-th element of `x` starts being valid. - `x::AbstractArray{Taylor1{U},N}`: Array of Taylor polynomials that interpolate the dependent variable as a function of the independent variable. """ -@auto_hash_equals struct TaylorInterpolant{T, U, N} - t0::T - t::AbstractVector{T} - x::AbstractArray{Taylor1{U}, N} +@auto_hash_equals struct TaylorInterpolant{T0<:Number, T<:AbstractVector{T0}, U<:Number, N, X<:AbstractArray{Taylor1{U}, N}} + t0::T0 + t::T + x::X # Inner constructor - function TaylorInterpolant{T, U, N}(t0::T, t::AbstractVector{T}, x::AbstractArray{Taylor1{U}, N}) where {T <: Real, U <: Number, N} + function TaylorInterpolant{T0, T, U, N, X}(t0::T0, t::T, x::X) where {T0<:Number, T<:AbstractVector{T0}, U<:Number, N, X<:AbstractArray{Taylor1{U}, N}} @assert size(x)[1] == length(t)-1 @assert issorted(t) || issorted(t, rev = true) - return new{T, U, N}(t0, t, x) + return new{T0, T, U, N, X}(t0, t, x) end end # Outer constructors -function TaylorInterpolant(t0::T, t::AbstractVector{T}, x::AbstractArray{Taylor1{U}, N}) where {T <: Real, U <: Number, N} - return TaylorInterpolant{T, U, N}(t0, t, x) +function TaylorInterpolant(t0::T0, t::T, x::X) where {T0<:Number, T<:AbstractVector{T0}, U<:Number, N, X<:AbstractArray{Taylor1{U}, N}} + return TaylorInterpolant{T0, T, U, N, X}(t0, t, x) end function TaylorInterpolant(t0::T, t::SubArray{T, 1}, x::SubArray{Taylor1{U}, N}) where {T <: Real, U <: Number, N} - return TaylorInterpolant{T, U, N}(t0, t.parent[t.indices...], x.parent[x.indices...]) -end + return TaylorInterpolant{T0, T, U, N, X}(t0, t.parent[t.indices...], x.parent[x.indices...]) +end -# Custom print -function show(io::IO, interp::TaylorInterpolant{T, U, N}) where {T, U, N} +# Custom print +function show(io::IO, interp::TaylorInterpolant{T0, T, U, N, X}) where {T0, T, U, N, X} t_range = minmax(interp.t0 + interp.t[1], interp.t0 + interp.t[end]) S = eltype(interp.x) if isone(N) print(io, "t: ", t_range, ", x: 1 ", S, " variable") - else + else L = size(interp.x, 2) print(io, "t: ", t_range, ", x: ", L, " ", S, " variables") end -end +end @doc raw""" convert(::Type{T}, interp::TaylorInterpolant) where {T <: Real} -Convert `inter.t0`, `inter.t` and coefficients of `interp.x` to type `T`. +Convert `inter.t0`, `inter.t` and coefficients of `interp.x` to type `T`. """ function convert(::Type{T}, interp::TaylorInterpolant) where {T <: Real} return TaylorInterpolant( - T(interp.t0), - T.(interp.t), + T(interp.t0), + T.(interp.t), map( x -> Taylor1( T.(x.coeffs) ), interp.x) ) -end +end @doc raw""" getinterpindex(tinterp::TaylorInterpolant{T,U,N}, t::V) where {T<:Real, U<:Number, V<:Number, N} @@ -65,7 +65,7 @@ end Return the index of `tinterp.t` corresponding to `t` and the time elapsed from `tinterp.t0` to `t`. """ -function getinterpindex(tinterp::TaylorInterpolant{T, U, N}, t::V) where {T<:Real, U<:Number, V<:Number, N} +function getinterpindex(tinterp::TaylorInterpolant{T0, T, U, N, X}, t::V) where {T0<:Number, T<:AbstractVector{T0}, U<:Number, N, X<:AbstractArray{Taylor1{U}, N}, V<:Union{T0, Taylor1{U}}} t00 = constant_term(constant_term(t)) # Current time tmin, tmax = minmax(tinterp.t[end], tinterp.t[1]) # Min and max time in tinterp Δt = t-tinterp.t0 # Time since start of tinterp @@ -82,18 +82,18 @@ function getinterpindex(tinterp::TaylorInterpolant{T, U, N}, t::V) where {T<:Rea end # Return index and elapsed time - return ind, Δt + return ind::Int, Δt::V end @doc raw""" numberofbodies(interp::TaylorInterpolant{T, U, 2}) where {T, U} -Return the number of bodies saved in a `TaylorInterpolant` produced by `PlanetaryEphemeris`. +Return the number of bodies saved in a `TaylorInterpolant` produced by `PlanetaryEphemeris`. """ numberofbodies(L::Int) = (L - 13) ÷ 6 numberofbodies(v::Vector{T}) where {T} = numberofbodies(length(v)) numberofbodies(m::Matrix{T}) where {T} = numberofbodies(size(m, 2)) -numberofbodies(interp::TaylorInterpolant{T, U, 2}) where {T, U} = numberofbodies(size(interp.x, 2)) +numberofbodies(interp::TaylorInterpolant{T0, T, U, 2, X}) where {T0, T, U, X} = numberofbodies(size(interp.x, 2)) # Function-like (callability) methods @@ -107,49 +107,49 @@ Evaluate `tinterp.x` at time `t`. See also [`getinterpindex`](@ref). """ -function (tinterp::TaylorInterpolant{T,U,1})(t::V) where {T<:Real, U<:Number, V<:Number} +function (tinterp::TaylorInterpolant{T0, T, U, 1, X})(t::V) where {T0, T, U, X, V<:Union{T0, Taylor1{U}}} # Get index of tinterp.x that interpolates at time t - ind, Δt = getinterpindex(tinterp, t) + ind::Int, Δt::V = getinterpindex(tinterp, t) # Time since the start of the ind-th timespan - δt = Δt-tinterp.t[ind] + δt::V = Δt-(tinterp.t[ind]) # Evaluate tinterp.x[ind] at δt - return tinterp.x[ind](δt) + return (tinterp.x[ind]::Taylor1{U})(δt) end -function (tinterp::TaylorInterpolant{T,U,2})(t::V) where {T<:Real, U<:Number, V<:Number} +function (tinterp::TaylorInterpolant{T0, T, U, 2, X})(t::V) where {T0, T, U, X, V<:Union{T0, Taylor1{U}}} # Get index of tinterp.x that interpolates at time t - ind, Δt = getinterpindex(tinterp, t) + ind::Int, Δt::V = getinterpindex(tinterp, t) # Time since the start of the ind-th timespan - δt = Δt-tinterp.t[ind] + δt::V = Δt-tinterp.t[ind] # Evaluate tinterp.x[ind] at δt return tinterp.x[ind,:](δt) end -function (tinterp::TaylorInterpolant{T,U,2})(target::Int, observer::Int, t::V) where {T<:Real, U<:Number, V<:Number} +function (tinterp::TaylorInterpolant{T0,T,U,2,X})(target::Int, observer::Int, t::T0) where {T0, T, U, X} # Number of bodies in tinterp N = numberofbodies(tinterp) # Ephemeris at time t eph_t = tinterp(t) - # Relative state vector + # Relative state vector if observer == 0 return eph_t[nbodyind(N, target)] else return eph_t[nbodyind(N, target)] - eph_t[nbodyind(N, observer)] - end + end end -(tinterp::TaylorInterpolant{T,U,2})(target::Int, t::V) where {T<:Real, U<:Number, V<:Number} = tinterp(target, 0, t) +(tinterp::TaylorInterpolant{T0,T,U,2,X})(target::Int, t::V) where {T0, T, U, X, V<:Union{T0, Taylor1{U}}} = tinterp(target, 0, t) @doc raw""" reverse(tinterp::TaylorInterpolant{T,U,N}) where {T<:Real, U<:Number, N} -Return a `TaylorInterpolant` object with the same information as `tinterp` but +Return a `TaylorInterpolant` object with the same information as `tinterp` but the independent variable reversed. See also [`TaylorInterpolant`](@ref). """ function reverse(tinterp::TaylorInterpolant{T,U,N}) where {T<:Real, U<:Number, N} - # tinterp end time is the new start time + # tinterp end time is the new start time tinterp_rev_t0 = tinterp.t[end] # reverse independent variable vector tinterp.t tinterp_rev_t = tinterp.t[end:-1:1] .- tinterp_rev_t0 @@ -164,26 +164,26 @@ end ttmtdb::Bool = false) where {T <: Real, U <: Number} selecteph(eph::TaylorInterpolant{T, U, 2}, i::Int) where {T <: Real, U <: Number} -Return a `TaylorInterpolant` with only the ephemeris of the bodies with indices `bodyind/i`. The keyword arguments allow to -include lunar euler angles and/or TT-TDB. +Return a `TaylorInterpolant` with only the ephemeris of the bodies with indices `bodyind/i`. The keyword arguments allow to +include lunar euler angles and/or TT-TDB. """ -function selecteph(eph::TaylorInterpolant{T, U, 2}, bodyind::AbstractVector{Int}; euler::Bool = false, - ttmtdb::Bool = false) where {T <: Real, U <: Number} +function selecteph(eph::TaylorInterpolant{T0, T, U, 2, X}, bodyind::AbstractVector{Int}; euler::Bool = false, + ttmtdb::Bool = false) where {T0, T, U, X} N = numberofbodies(eph) idxs = nbodyind(N, bodyind) - if euler + if euler idxs = union(idxs, 6N+1:6N+12) - end + end if ttmtdb idxs = union(idxs, 6N+13) - end + end x = eph.x[:, idxs] - return TaylorInterpolant{T, U, 2}(eph.t0, eph.t, x) -end + return TaylorInterpolant{T0, T, U, 2, X}(eph.t0, eph.t, x) +end -selecteph(eph::TaylorInterpolant{T, U, 2}, i::Int) where {T <: Real, U <: Number} = selecteph(eph, i:i) +selecteph(eph::TaylorInterpolant{T0, T, U, 2, X}, i::Int) where {T0, T, U, X} = selecteph(eph, i:i) -function join(bwd::TaylorInterpolant{T, U, 2}, fwd::TaylorInterpolant{T, U, 2}) where {T, U} +function join(bwd::TaylorInterpolant{T0, T, U, 2, X}, fwd::TaylorInterpolant{T0, T, U, 2, X}) where {T0, T, U, X} @assert bwd.t0 == fwd.t0 "Initial time must be the same for both TaylorInterpolant" order_bwd = get_order(bwd.x[1, 1]) order_fwd = get_order(fwd.x[1, 1]) @@ -195,7 +195,7 @@ function join(bwd::TaylorInterpolant{T, U, 2}, fwd::TaylorInterpolant{T, U, 2}) x = vcat(reverse(bwd.x, dims = 1), fwd.x) return TaylorInterpolant(t0, t, x) -end +end @doc raw""" kmsec2auday(pv) @@ -219,19 +219,19 @@ function auday2kmsec(pv) return pv end -# Custom serialization +# Custom serialization @doc raw""" PlanetaryEphemerisSerialization{T} -Custom serialization struct to save a `TaylorInterpolant{T, T, 2}` to a `.jld2` file. +Custom serialization struct to save a `TaylorInterpolant{T, T, 2}` to a `.jld2` file. # Fields - `order::Int`: order of Taylor polynomials. -- `dims::Tuple{Int, Int}`: matrix dimensions. +- `dims::Tuple{Int, Int}`: matrix dimensions. - `t0::T`: initial time. -- `t::Vector{T}`: vector of times. -- `x::Vector{T}`: vector of coefficients. +- `t::Vector{T}`: vector of times. +- `x::Vector{T}`: vector of coefficients. """ struct PlanetaryEphemerisSerialization{T} order::Int @@ -239,13 +239,13 @@ struct PlanetaryEphemerisSerialization{T} t0::T t::Vector{T} x::Vector{T} -end +end # Tell JLD2 to save TaylorInterpolant{T, T, 2} as PlanetaryEphemerisSerialization{T} -writeas(::Type{TaylorInterpolant{T, T, 2}}) where {T} = PlanetaryEphemerisSerialization{T} +writeas(::Type{TaylorInterpolant{T0, T, U, 2, X}}) where {T0, T, U, X} = PlanetaryEphemerisSerialization{T} # Convert method to write .jld2 files -function convert(::Type{PlanetaryEphemerisSerialization{T}}, eph::TaylorInterpolant{T, T, 2}) where {T} +function convert(::Type{PlanetaryEphemerisSerialization{T}}, eph::TaylorInterpolant{T0, T, U, 2, X}) where {T0, T, U, X} # Taylor polynomials order order = eph.x[1, 1].order # Number of coefficients in each polynomial @@ -254,18 +254,18 @@ function convert(::Type{PlanetaryEphemerisSerialization{T}}, eph::TaylorInterpol dims = size(eph.x) # Number of elements in matrix N = dims[1] * dims[2] - # Vector of coefficients + # Vector of coefficients x = Vector{T}(undef, k * N) - # Save coefficients + # Save coefficients for i in 1:N x[(i-1)*k+1 : i*k] = eph.x[i].coeffs - end + end return PlanetaryEphemerisSerialization{T}(order, dims, eph.t0, eph.t, x) -end +end # Convert method to read .jld2 files -function convert(::Type{TaylorInterpolant{T, T, 2}}, eph::PlanetaryEphemerisSerialization{T}) where {T} +function convert(::Type{TaylorInterpolant{T0, T, U, 2, X}}, eph::PlanetaryEphemerisSerialization{T}) where {T0, T, U, X} # Taylor polynomials order order = eph.order # Number of coefficients in each polynomial @@ -273,22 +273,22 @@ function convert(::Type{TaylorInterpolant{T, T, 2}}, eph::PlanetaryEphemerisSeri # Matrix dimensions dims = eph.dims # Number of elements in matrix - N = dims[1] * dims[2] - # Matrix of Taylor polynomials + N = dims[1] * dims[2] + # Matrix of Taylor polynomials x = Matrix{Taylor1{T}}(undef, dims[1], dims[2]) - # Reconstruct Taylor polynomials + # Reconstruct Taylor polynomials for i in 1:N x[i] = Taylor1{T}(eph.x[(i-1)*k+1 : i*k], order) - end + end - return TaylorInterpolant{T, T, 2}(eph.t0, eph.t, x) -end + return TaylorInterpolant{T0, T, U, 2, X}(eph.t0, eph.t, x) +end @doc raw""" Taylor1Serialization{T} -Custom serialization struct used in previous versions (<= 0.4) of `PlanetaryEphemeris`. Currently, it is only used to -read old `.jld2` files. +Custom serialization struct used in previous versions (<= 0.4) of `PlanetaryEphemeris`. Currently, it is only used to +read old `.jld2` files. """ struct Taylor1Serialization{T} x::Vector{T} From 65f431d1855284feaa32a699b61b0dc01b230455 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jorge=20A=2E=20P=C3=A9rez-Hern=C3=A1ndez?= Date: Thu, 25 Jan 2024 23:08:22 +0100 Subject: [PATCH 02/10] Update CI --- .github/workflows/CI.yml | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 014f9fa..b59ae5f 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -1,8 +1,9 @@ name: CI on: push: - branches: - - main + paths-ignore: + - 'LICENSE.md' + - 'README.md' tags: ['*'] pull_request: concurrency: From 8ede385ee9b711e16fb88d37b8eaaa0759acfe9b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jorge=20A=2E=20P=C3=A9rez-Hern=C3=A1ndez?= Date: Sat, 27 Jan 2024 20:20:43 +0100 Subject: [PATCH 03/10] Fix typing issues in TaylorInterpolant and update tests --- src/PlanetaryEphemeris.jl | 6 +- src/{ => dynamics}/dynamical_model.jl | 0 src/{ => dynamics}/jetcoeffs.jl | 0 src/dynamics/trivial.jl | 6 + src/interpolation.jl | 66 ++++---- src/propagation.jl | 6 +- test/propagation.jl | 224 +++++++++++++------------- 7 files changed, 160 insertions(+), 148 deletions(-) rename src/{ => dynamics}/dynamical_model.jl (100%) rename src/{ => dynamics}/jetcoeffs.jl (100%) create mode 100644 src/dynamics/trivial.jl diff --git a/src/PlanetaryEphemeris.jl b/src/PlanetaryEphemeris.jl index 5beae4b..b2700d2 100644 --- a/src/PlanetaryEphemeris.jl +++ b/src/PlanetaryEphemeris.jl @@ -26,12 +26,14 @@ import JLD2: writeas include("constants.jl") include("jpl-de-430-431-earth-orientation-model.jl") include("initial_conditions.jl") -include("dynamical_model.jl") -include("jetcoeffs.jl") include("interpolation.jl") include("propagation.jl") include("osculating.jl") include("barycenter.jl") #include("precompile.jl") +include("dynamics/trivial.jl") +include("dynamics/dynamical_model.jl") +include("dynamics/jetcoeffs.jl") + end # module diff --git a/src/dynamical_model.jl b/src/dynamics/dynamical_model.jl similarity index 100% rename from src/dynamical_model.jl rename to src/dynamics/dynamical_model.jl diff --git a/src/jetcoeffs.jl b/src/dynamics/jetcoeffs.jl similarity index 100% rename from src/jetcoeffs.jl rename to src/dynamics/jetcoeffs.jl diff --git a/src/dynamics/trivial.jl b/src/dynamics/trivial.jl new file mode 100644 index 0000000..80ae0b9 --- /dev/null +++ b/src/dynamics/trivial.jl @@ -0,0 +1,6 @@ +@taylorize function trivialdynamics!(dq, q, params, t) + Threads.@threads for j in eachindex(dq) + dq[j] = zero(q[j]) + end + nothing +end diff --git a/src/interpolation.jl b/src/interpolation.jl index 886a5d7..399d2bc 100644 --- a/src/interpolation.jl +++ b/src/interpolation.jl @@ -13,29 +13,29 @@ time ``x(t)``; or a lunar core Euler angle as a function of time ``\theta_c(t)`` - `t::AbstractVector{T}`: Vector of time instances when the timespan of the ``i``-th element of `x` ends and the ``(i+1)``-th element of `x` starts being valid. - `x::AbstractArray{Taylor1{U},N}`: Array of Taylor polynomials that interpolate the dependent variable as a function of the independent variable. """ -@auto_hash_equals struct TaylorInterpolant{T0<:Number, T<:AbstractVector{T0}, U<:Number, N, X<:AbstractArray{Taylor1{U}, N}} - t0::T0 - t::T +@auto_hash_equals struct TaylorInterpolant{T<:Number, U<:Number, N, VT<:AbstractVector{T}, X<:AbstractArray{Taylor1{U}, N}} + t0::T + t::VT x::X # Inner constructor - function TaylorInterpolant{T0, T, U, N, X}(t0::T0, t::T, x::X) where {T0<:Number, T<:AbstractVector{T0}, U<:Number, N, X<:AbstractArray{Taylor1{U}, N}} + function TaylorInterpolant{T, U, N, VT, X}(t0::T, t::VT, x::X) where {T<:Number, U<:Number, N, VT<:AbstractVector{T}, X<:AbstractArray{Taylor1{U}, N}} @assert size(x)[1] == length(t)-1 @assert issorted(t) || issorted(t, rev = true) - return new{T0, T, U, N, X}(t0, t, x) + return new{T, U, N, VT, X}(t0, t, x) end end # Outer constructors -function TaylorInterpolant(t0::T0, t::T, x::X) where {T0<:Number, T<:AbstractVector{T0}, U<:Number, N, X<:AbstractArray{Taylor1{U}, N}} - return TaylorInterpolant{T0, T, U, N, X}(t0, t, x) +function TaylorInterpolant(t0::T, t::VT, x::X) where {T<:Number, U<:Number, N, VT<:AbstractVector{T}, X<:AbstractArray{Taylor1{U}, N}} + return TaylorInterpolant{T, U, N, VT, X}(t0, t, x) end -function TaylorInterpolant(t0::T, t::SubArray{T, 1}, x::SubArray{Taylor1{U}, N}) where {T <: Real, U <: Number, N} - return TaylorInterpolant{T0, T, U, N, X}(t0, t.parent[t.indices...], x.parent[x.indices...]) +function TaylorInterpolant(t0::T, t::SubArray{T, 1}, x::SubArray{Taylor1{U}, N}) where {T<:Number, U<:Number, N} + return TaylorInterpolant{T, U, N, Vector{T}, Array{Taylor1{U}, N}}(t0, t.parent[t.indices...], x.parent[x.indices...]) end # Custom print -function show(io::IO, interp::TaylorInterpolant{T0, T, U, N, X}) where {T0, T, U, N, X} +function show(io::IO, interp::T) where {U, V, N, T<:TaylorInterpolant{U,V,N}} t_range = minmax(interp.t0 + interp.t[1], interp.t0 + interp.t[end]) S = eltype(interp.x) if isone(N) @@ -65,7 +65,7 @@ end Return the index of `tinterp.t` corresponding to `t` and the time elapsed from `tinterp.t0` to `t`. """ -function getinterpindex(tinterp::TaylorInterpolant{T0, T, U, N, X}, t::V) where {T0<:Number, T<:AbstractVector{T0}, U<:Number, N, X<:AbstractArray{Taylor1{U}, N}, V<:Union{T0, Taylor1{U}}} +function getinterpindex(tinterp::TaylorInterpolant{T,U,N}, t::TT) where {T,U,N,TT<:Union{T, Taylor1{U}}} t00 = constant_term(constant_term(t)) # Current time tmin, tmax = minmax(tinterp.t[end], tinterp.t[1]) # Min and max time in tinterp Δt = t-tinterp.t0 # Time since start of tinterp @@ -82,7 +82,7 @@ function getinterpindex(tinterp::TaylorInterpolant{T0, T, U, N, X}, t::V) where end # Return index and elapsed time - return ind::Int, Δt::V + return ind::Int, Δt::TT end @doc raw""" @@ -93,7 +93,7 @@ Return the number of bodies saved in a `TaylorInterpolant` produced by `Planetar numberofbodies(L::Int) = (L - 13) ÷ 6 numberofbodies(v::Vector{T}) where {T} = numberofbodies(length(v)) numberofbodies(m::Matrix{T}) where {T} = numberofbodies(size(m, 2)) -numberofbodies(interp::TaylorInterpolant{T0, T, U, 2, X}) where {T0, T, U, X} = numberofbodies(size(interp.x, 2)) +numberofbodies(interp::TaylorInterpolant) = numberofbodies(size(interp.x, 2)) # Function-like (callability) methods @@ -107,25 +107,25 @@ Evaluate `tinterp.x` at time `t`. See also [`getinterpindex`](@ref). """ -function (tinterp::TaylorInterpolant{T0, T, U, 1, X})(t::V) where {T0, T, U, X, V<:Union{T0, Taylor1{U}}} +function (tinterp::TaylorInterpolant{T, U, 1})(t::TT) where {T, U, TT<:Union{T, Taylor1{U}}} # Get index of tinterp.x that interpolates at time t - ind::Int, Δt::V = getinterpindex(tinterp, t) + ind::Int, Δt::TT = getinterpindex(tinterp, t) # Time since the start of the ind-th timespan - δt::V = Δt-(tinterp.t[ind]) + δt::TT = Δt-(tinterp.t[ind]) # Evaluate tinterp.x[ind] at δt return (tinterp.x[ind]::Taylor1{U})(δt) end -function (tinterp::TaylorInterpolant{T0, T, U, 2, X})(t::V) where {T0, T, U, X, V<:Union{T0, Taylor1{U}}} +function (tinterp::TaylorInterpolant{T, U, 2})(t::TT) where {T, U, TT<:Union{T, Taylor1{U}}} # Get index of tinterp.x that interpolates at time t - ind::Int, Δt::V = getinterpindex(tinterp, t) + ind::Int, Δt::TT = getinterpindex(tinterp, t) # Time since the start of the ind-th timespan - δt::V = Δt-tinterp.t[ind] + δt::TT = Δt-tinterp.t[ind] # Evaluate tinterp.x[ind] at δt return tinterp.x[ind,:](δt) end -function (tinterp::TaylorInterpolant{T0,T,U,2,X})(target::Int, observer::Int, t::T0) where {T0, T, U, X} +function (tinterp::TaylorInterpolant{T, U, 2})(target::Int, observer::Int, t::T) where {T, U} # Number of bodies in tinterp N = numberofbodies(tinterp) # Ephemeris at time t @@ -138,7 +138,7 @@ function (tinterp::TaylorInterpolant{T0,T,U,2,X})(target::Int, observer::Int, t: end end -(tinterp::TaylorInterpolant{T0,T,U,2,X})(target::Int, t::V) where {T0, T, U, X, V<:Union{T0, Taylor1{U}}} = tinterp(target, 0, t) +(tinterp::TaylorInterpolant{T,U,2})(target::Int, t::TT) where {T, U, TT<:Union{T, Taylor1{U}}} = tinterp(target, 0, t) @doc raw""" reverse(tinterp::TaylorInterpolant{T,U,N}) where {T<:Real, U<:Number, N} @@ -167,8 +167,8 @@ end Return a `TaylorInterpolant` with only the ephemeris of the bodies with indices `bodyind/i`. The keyword arguments allow to include lunar euler angles and/or TT-TDB. """ -function selecteph(eph::TaylorInterpolant{T0, T, U, 2, X}, bodyind::AbstractVector{Int}; euler::Bool = false, - ttmtdb::Bool = false) where {T0, T, U, X} +function selecteph(eph::TaylorInterpolant{T, U, 2}, bodyind::AbstractVector{Int}; euler::Bool = false, + ttmtdb::Bool = false) where {T, U} N = numberofbodies(eph) idxs = nbodyind(N, bodyind) if euler @@ -177,21 +177,21 @@ function selecteph(eph::TaylorInterpolant{T0, T, U, 2, X}, bodyind::AbstractVect if ttmtdb idxs = union(idxs, 6N+13) end - x = eph.x[:, idxs] - return TaylorInterpolant{T0, T, U, 2, X}(eph.t0, eph.t, x) + x = view(eph.x, :, idxs) + return TaylorInterpolant(eph.t0, eph.t, x) end -selecteph(eph::TaylorInterpolant{T0, T, U, 2, X}, i::Int) where {T0, T, U, X} = selecteph(eph, i:i) +selecteph(eph::TaylorInterpolant, i::Int) = selecteph(eph, i:i) -function join(bwd::TaylorInterpolant{T0, T, U, 2, X}, fwd::TaylorInterpolant{T0, T, U, 2, X}) where {T0, T, U, X} +function join(bwd::TaylorInterpolant{T, U, 2}, fwd::TaylorInterpolant{T, U, 2}) where {T, U} @assert bwd.t0 == fwd.t0 "Initial time must be the same for both TaylorInterpolant" order_bwd = get_order(bwd.x[1, 1]) order_fwd = get_order(fwd.x[1, 1]) @assert order_bwd == order_fwd "Expansion order must be the same for both TaylorInterpolant" t0 = bwd.t0 + bwd.t[end] - t_1 = abs.(bwd.t) - t = vcat(t_1, t_1[end] .+ fwd.t[2:end]) + t1 = abs.(bwd.t) + t = vcat(t1, t1[end] .+ fwd.t[2:end]) x = vcat(reverse(bwd.x, dims = 1), fwd.x) return TaylorInterpolant(t0, t, x) @@ -242,10 +242,10 @@ struct PlanetaryEphemerisSerialization{T} end # Tell JLD2 to save TaylorInterpolant{T, T, 2} as PlanetaryEphemerisSerialization{T} -writeas(::Type{TaylorInterpolant{T0, T, U, 2, X}}) where {T0, T, U, X} = PlanetaryEphemerisSerialization{T} +writeas(::Type{TaylorInterpolant{T, T, 2, Vector{T}, Matrix{Taylor1{T}}}}) where {T<:Real} = PlanetaryEphemerisSerialization{T} # Convert method to write .jld2 files -function convert(::Type{PlanetaryEphemerisSerialization{T}}, eph::TaylorInterpolant{T0, T, U, 2, X}) where {T0, T, U, X} +function convert(::Type{PlanetaryEphemerisSerialization{T}}, eph::TaylorInterpolant{T, T, 2, Vector{T}, Matrix{Taylor1{T}}}) where {T <: Real} # Taylor polynomials order order = eph.x[1, 1].order # Number of coefficients in each polynomial @@ -265,7 +265,7 @@ function convert(::Type{PlanetaryEphemerisSerialization{T}}, eph::TaylorInterpol end # Convert method to read .jld2 files -function convert(::Type{TaylorInterpolant{T0, T, U, 2, X}}, eph::PlanetaryEphemerisSerialization{T}) where {T0, T, U, X} +function convert(::Type{TaylorInterpolant{T, T, 2, Vector{T}, Array{Taylor1{T},2}}}, eph::PlanetaryEphemerisSerialization{T}) where {T<:Real} # Taylor polynomials order order = eph.order # Number of coefficients in each polynomial @@ -281,7 +281,7 @@ function convert(::Type{TaylorInterpolant{T0, T, U, 2, X}}, eph::PlanetaryEpheme x[i] = Taylor1{T}(eph.x[(i-1)*k+1 : i*k], order) end - return TaylorInterpolant{T0, T, U, 2, X}(eph.t0, eph.t, x) + return TaylorInterpolant{T, T, 2, Vector{T}, Matrix{Taylor1{T}}}(eph.t0, eph.t, x) end @doc raw""" diff --git a/src/propagation.jl b/src/propagation.jl index db134c8..52557ac 100644 --- a/src/propagation.jl +++ b/src/propagation.jl @@ -193,13 +193,9 @@ for V_dense in (:(Val{true}), :(Val{false})) maxsteps, parse_eqs) if $V_dense == Val{true} - - return TaylorInterpolant{T, T, 2}(jd0 - J2000, sol[1], sol[3]) - + return TaylorInterpolant{T, T, 2, SubArray{T,1}, SubArray{Taylor1{T},2}}(jd0 - J2000, sol[1], sol[3]) else - return sol - end end diff --git a/test/propagation.jl b/test/propagation.jl index 43ca1a6..aca9986 100644 --- a/test/propagation.jl +++ b/test/propagation.jl @@ -8,50 +8,13 @@ using Test using PlanetaryEphemeris: initialcond, ssic_1969, ssic_2000, astic_1969, astic_2000 -@testset "Initial conditions" begin - - # Number of asteroids - local N_ast = 343 - # Total number of bodies (solar system + asteroids) - local N = 11 + N_ast - # Check initial conditions global variables - @test size(ssic_1969) == (11, 6) - @test size(ssic_2000) == (11, 6) - - @test isa(ssic_1969, Matrix{Float64}) - @test isa(ssic_2000, Matrix{Float64}) - - @test size(astic_1969) == (N_ast, 6) - @test size(astic_2000) == (N_ast, 6) - - @test isa(astic_1969, Matrix{Float64}) - @test isa(astic_2000, Matrix{Float64}) - - # Number of bodies - local N = 11 + 343 - # Check output of initialcond - for jd0 in datetime2julian.( ( DateTime(1969,6,28,0,0,0), DateTime(2000,1,1,12) ) ) - - q0_64 = initialcond(N, jd0) - q0_128 = initialcond(N, convert(Float128, jd0) ) - - @test isa(q0_64, Vector{Float64}) - @test isa(q0_128, Vector{Float128}) - @test length(q0_64) == length(q0_128) == 6*N+13 - @test q0_64 == q0_128 - - end - -end - using Downloads using SPICE: furnsh, spkgeo using PlanetaryEphemeris: order, abstol using LinearAlgebra: norm -@testset "Propagation" begin - @show Threads.nthreads() +@testset "Propagation" begin # Number of asteroids local N_ast = 343 @@ -68,76 +31,121 @@ using LinearAlgebra: norm # Dynamical function local dynamics = DE430! - # Float64 - - # Test integration - sol64 = propagate(100, jd0, nyears, dense; dynamics, order, abstol) - # Save solution - filename = selecteph2jld2(sol64, bodyind, nyears) - # Recovered solution - recovered_sol64 = JLD2.load(filename, "ss16ast_eph") - - @test selecteph(sol64, bodyind, euler = true, ttmtdb = true) == recovered_sol64 - - LSK = "https://naif.jpl.nasa.gov/pub/naif/generic_kernels/lsk/naif0012.tls" - TTmTDBK = "https://ssd.jpl.nasa.gov/ftp/eph/planets/bsp/TTmTDB.de430.19feb2015.bsp" - SPK = "https://ssd.jpl.nasa.gov/ftp/eph/planets/bsp/de430_1850-2150.bsp" - - # Download kernels - Downloads.download(LSK, "naif0012.tls") - Downloads.download(SPK, "de430_1850-2150.bsp") - Downloads.download(TTmTDBK, "TTmTDB.de430.19feb2015.bsp") - - # Load kernels - furnsh("naif0012.tls", "de430_1850-2150.bsp", "TTmTDB.de430.19feb2015.bsp") - - ttmtdb_pe = TaylorInterpolant(sol64.t0, sol64.t, sol64.x[:, 6N+13]) # TT-TDB - posvel_pe_su = selecteph(sol64,su) # Sun - posvel_pe_ea = selecteph(sol64,ea) # Earth - posvel_pe_mo = selecteph(sol64,mo) # Moon - posvel_pe_ma = selecteph(sol64,6) # Mars - posvel_pe_ju = selecteph(sol64,7) # Jupiter - - ttmtdb_jpl(et) = spkgeo(1000000001, et, "J2000", 1000000000)[1][1] # TT-TDB - posvel_jpl_su(et) = kmsec2auday(spkgeo(10, et, "J2000", 0)[1]) # Sun - posvel_jpl_ea(et) = kmsec2auday(spkgeo(399, et, "J2000", 0)[1]) # Earth - posvel_jpl_mo(et) = kmsec2auday(spkgeo(301, et, "J2000", 0)[1]) # Moon - posvel_jpl_ma(et) = kmsec2auday(spkgeo(4, et, "J2000", 0)[1]) # Mars - posvel_jpl_ju(et) = kmsec2auday(spkgeo(5, et, "J2000", 0)[1]) # Jupiter - - tv = range(sol64.t0, sol64.t[end], 10) - for t in tv - et = t * daysec - @show t, et - @show abs(ttmtdb_jpl(et) - ttmtdb_pe(t)) - @show norm(posvel_jpl_su(et) - posvel_pe_su(t), Inf) - @show norm(posvel_jpl_ea(et) - posvel_pe_ea(t), Inf) - @show norm(posvel_jpl_mo(et) - posvel_pe_mo(t), Inf) - @show norm(posvel_jpl_ma(et) - posvel_pe_ma(t), Inf) - @show norm(posvel_jpl_ju(et) - posvel_pe_ju(t), Inf) - - @test abs(ttmtdb_jpl(et) - ttmtdb_pe(t)) < 1e-12 - @test norm(posvel_jpl_su(et) - posvel_pe_su(t), Inf) < 1e-14 - @test norm(posvel_jpl_ea(et) - posvel_pe_ea(t), Inf) < 1e-11 - @test norm(posvel_jpl_mo(et) - posvel_pe_mo(t), Inf) < 1e-11 - @test norm(posvel_jpl_ma(et) - posvel_pe_ma(t), Inf) < 1e-12 - @test norm(posvel_jpl_ju(et) - posvel_pe_ju(t), Inf) < 1e-13 + @testset "Initial conditions" begin + # Check initial conditions global variables + @test size(ssic_1969) == (11, 6) + @test size(ssic_2000) == (11, 6) + + @test isa(ssic_1969, Matrix{Float64}) + @test isa(ssic_2000, Matrix{Float64}) + + @test size(astic_1969) == (N_ast, 6) + @test size(astic_2000) == (N_ast, 6) + + @test isa(astic_1969, Matrix{Float64}) + @test isa(astic_2000, Matrix{Float64}) + + # Check output of initialcond + for jd0 in datetime2julian.( ( DateTime(1969,6,28,0,0,0), DateTime(2000,1,1,12) ) ) + + q0_64 = initialcond(N, jd0) + q0_128 = initialcond(N, convert(Float128, jd0) ) + + @test isa(q0_64, Vector{Float64}) + @test isa(q0_128, Vector{Float128}) + @test length(q0_64) == length(q0_128) == 6*N+13 + @test q0_64 == q0_128 + + end + end + + @testset "TaylorInterpolant" begin + # Test integration + sol = propagate(5, jd0, nyears, dense; dynamics=PlanetaryEphemeris.trivialdynamics!, order, abstol) + @test sol isa TaylorInterpolant{Float64,Float64,2} + q0 = initialcond(N, jd0) + @test sol(sol.t0) == q0 + @test sol.t0 == 0.0 + @test length(sol.t) == size(sol.x, 1) + 1 + @test length(q0) == size(sol.x, 2) + end + + @testset "Propagation: DE430 dynamical model" begin + + @show Threads.nthreads() + + # Float64 + + # Test integration + sol64 = propagate(100, jd0, nyears, dense; dynamics, order, abstol) + # Save solution + filename = selecteph2jld2(sol64, bodyind, nyears) + # Recovered solution + recovered_sol64 = JLD2.load(filename, "ss16ast_eph") + + @test selecteph(sol64, bodyind, euler = true, ttmtdb = true) == recovered_sol64 + + LSK = "https://naif.jpl.nasa.gov/pub/naif/generic_kernels/lsk/naif0012.tls" + TTmTDBK = "https://ssd.jpl.nasa.gov/ftp/eph/planets/bsp/TTmTDB.de430.19feb2015.bsp" + SPK = "https://ssd.jpl.nasa.gov/ftp/eph/planets/bsp/de430_1850-2150.bsp" + + # Download kernels + Downloads.download(LSK, "naif0012.tls") + Downloads.download(SPK, "de430_1850-2150.bsp") + Downloads.download(TTmTDBK, "TTmTDB.de430.19feb2015.bsp") + + # Load kernels + furnsh("naif0012.tls", "de430_1850-2150.bsp", "TTmTDB.de430.19feb2015.bsp") + + ttmtdb_pe = TaylorInterpolant(sol64.t0, sol64.t, sol64.x[:, 6N+13]) # TT-TDB + posvel_pe_su = selecteph(sol64,su) # Sun + posvel_pe_ea = selecteph(sol64,ea) # Earth + posvel_pe_mo = selecteph(sol64,mo) # Moon + posvel_pe_ma = selecteph(sol64,6) # Mars + posvel_pe_ju = selecteph(sol64,7) # Jupiter + + ttmtdb_jpl(et) = spkgeo(1000000001, et, "J2000", 1000000000)[1][1] # TT-TDB + posvel_jpl_su(et) = kmsec2auday(spkgeo(10, et, "J2000", 0)[1]) # Sun + posvel_jpl_ea(et) = kmsec2auday(spkgeo(399, et, "J2000", 0)[1]) # Earth + posvel_jpl_mo(et) = kmsec2auday(spkgeo(301, et, "J2000", 0)[1]) # Moon + posvel_jpl_ma(et) = kmsec2auday(spkgeo(4, et, "J2000", 0)[1]) # Mars + posvel_jpl_ju(et) = kmsec2auday(spkgeo(5, et, "J2000", 0)[1]) # Jupiter + + tv = range(sol64.t0, sol64.t[end], 10) + for t in tv + et = t * daysec + @show t, et + @show abs(ttmtdb_jpl(et) - ttmtdb_pe(t)) + @show norm(posvel_jpl_su(et) - posvel_pe_su(t), Inf) + @show norm(posvel_jpl_ea(et) - posvel_pe_ea(t), Inf) + @show norm(posvel_jpl_mo(et) - posvel_pe_mo(t), Inf) + @show norm(posvel_jpl_ma(et) - posvel_pe_ma(t), Inf) + @show norm(posvel_jpl_ju(et) - posvel_pe_ju(t), Inf) + + @test abs(ttmtdb_jpl(et) - ttmtdb_pe(t)) < 1e-12 + @test norm(posvel_jpl_su(et) - posvel_pe_su(t), Inf) < 1e-14 + @test norm(posvel_jpl_ea(et) - posvel_pe_ea(t), Inf) < 1e-11 + @test norm(posvel_jpl_mo(et) - posvel_pe_mo(t), Inf) < 1e-11 + @test norm(posvel_jpl_ma(et) - posvel_pe_ma(t), Inf) < 1e-12 + @test norm(posvel_jpl_ju(et) - posvel_pe_ju(t), Inf) < 1e-13 + end + + # Remove files + rm.((filename, "naif0012.tls", "de430_1850-2150.bsp", "TTmTDB.de430.19feb2015.bsp")) + + # Float 128 + #= + # Test integration + sol128 = propagate(1, Float128(jd0), nyears, dense; dynamics = dynamics, order = order, abstol = abstol) + # Save solution + filename = selecteph2jld2(sol128, bodyind, nyears) + # Recovered solution + recovered_sol128 = JLD2.load(filename, "ss16ast_eph") + # Remove file + rm(filename) + + @test selecteph(sol128, bodyind, euler = true, ttmtdb = true) == recovered_sol128 + =# end - # Remove files - rm.((filename, "naif0012.tls", "de430_1850-2150.bsp", "TTmTDB.de430.19feb2015.bsp")) - - # Float 128 - #= - # Test integration - sol128 = propagate(1, Float128(jd0), nyears, dense; dynamics = dynamics, order = order, abstol = abstol) - # Save solution - filename = selecteph2jld2(sol128, bodyind, nyears) - # Recovered solution - recovered_sol128 = JLD2.load(filename, "ss16ast_eph") - # Remove file - rm(filename) - - @test selecteph(sol128, bodyind, euler = true, ttmtdb = true) == recovered_sol128 - =# end From 687d220333c57de99e64f3eb4b10e7b75efb21ea Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jorge=20A=2E=20P=C3=A9rez-Hern=C3=A1ndez?= Date: Sat, 27 Jan 2024 22:21:22 +0100 Subject: [PATCH 04/10] Add TaylorInterpolant function-like behavior for Taylor1{TaylorN{Float64}} --- src/dynamics/trivial.jl | 2 + src/initial_conditions.jl | 2 + src/interpolation.jl | 47 +++-- src/jpl-de-430-431-earth-orientation-model.jl | 90 +++++----- src/osculating.jl | 72 ++++---- src/propagation.jl | 2 + test/Project.toml | 2 + test/propagation.jl | 161 +++++++++--------- 8 files changed, 203 insertions(+), 175 deletions(-) diff --git a/src/dynamics/trivial.jl b/src/dynamics/trivial.jl index 80ae0b9..4a62abb 100644 --- a/src/dynamics/trivial.jl +++ b/src/dynamics/trivial.jl @@ -1,3 +1,5 @@ +# This file is part of the PlanetaryEphemeris.jl package; MIT licensed + @taylorize function trivialdynamics!(dq, q, params, t) Threads.@threads for j in eachindex(dq) dq[j] = zero(q[j]) diff --git a/src/initial_conditions.jl b/src/initial_conditions.jl index e555d2a..4bffa6f 100644 --- a/src/initial_conditions.jl +++ b/src/initial_conditions.jl @@ -1,3 +1,5 @@ +# This file is part of the PlanetaryEphemeris.jl package; MIT licensed + # Planets (+ Sun & Moon) initial conditions file const ssic_1969_fname = joinpath( src_path, "ss11ic_1969Jun28.txt" ) const ssic_1969 = readdlm( ssic_1969_fname ) diff --git a/src/interpolation.jl b/src/interpolation.jl index 399d2bc..7e3dc6e 100644 --- a/src/interpolation.jl +++ b/src/interpolation.jl @@ -65,7 +65,7 @@ end Return the index of `tinterp.t` corresponding to `t` and the time elapsed from `tinterp.t0` to `t`. """ -function getinterpindex(tinterp::TaylorInterpolant{T,U,N}, t::TT) where {T,U,N,TT<:Union{T, Taylor1{U}}} +function getinterpindex(tinterp::TaylorInterpolant{T,U,N}, t::TT) where {T,U,N,TT<:Union{T, Taylor1{U}, Taylor1{TaylorN{U}}}} t00 = constant_term(constant_term(t)) # Current time tmin, tmax = minmax(tinterp.t[end], tinterp.t[1]) # Min and max time in tinterp Δt = t-tinterp.t0 # Time since start of tinterp @@ -81,8 +81,11 @@ function getinterpindex(tinterp::TaylorInterpolant{T,U,N}, t::TT) where {T,U,N,T ind = searchsortedlast(tinterp.t, Δt00, rev=true) end - # Return index and elapsed time - return ind::Int, Δt::TT + # Time since the start of the ind-th timestep + δt = Δt-tinterp.t[ind] + + # Return index and elapsed time since i-th timestep + return ind::Int, δt::TT end @doc raw""" @@ -98,31 +101,39 @@ numberofbodies(interp::TaylorInterpolant) = numberofbodies(size(interp.x, 2)) # Function-like (callability) methods @doc raw""" - (tinterp::TaylorInterpolant{T,U,1})(t::V) where {T<:Real, U<:Number, V<:Number} - (tinterp::TaylorInterpolant{T,U,2})(t::V) where {T<:Real, U<:Number, V<:Number} - (tinterp::TaylorInterpolant{T,U,2})(target::Int, t::V) where {T<:Real, U<:Number, V<:Number} - (tinterp::TaylorInterpolant{T,U,2})(target::Int, observer::Int, t::V) where {T<:Real, U<:Number, V<:Number} + (tinterp::TaylorInterpolant{T, U, 1})(t::T) where {T, U} + (tinterp::TaylorInterpolant{T, U, 1})(t::TT) where {T, U, TT<:Union{Taylor1{U},Taylor1{TaylorN{U}}}} + (tinterp::TaylorInterpolant{T, U, 2})(t::T) where {T, U} + (tinterp::TaylorInterpolant{T, U, 2})(t::TT) where {T, U, TT<:Union{Taylor1{U},Taylor1{TaylorN{U}}}} Evaluate `tinterp.x` at time `t`. See also [`getinterpindex`](@ref). """ -function (tinterp::TaylorInterpolant{T, U, 1})(t::TT) where {T, U, TT<:Union{T, Taylor1{U}}} +function (tinterp::TaylorInterpolant{T, U, 1})(t::T) where {T, U} + # Get index of tinterp.x that interpolates at time t + ind::Int, δt::T = getinterpindex(tinterp, t) + # Evaluate tinterp.x[ind] at δt + return (tinterp.x[ind])(δt)::U +end +function (tinterp::TaylorInterpolant{T, U, 1})(t::TT) where {T, U, TT<:Union{Taylor1{U},Taylor1{TaylorN{U}}}} # Get index of tinterp.x that interpolates at time t - ind::Int, Δt::TT = getinterpindex(tinterp, t) - # Time since the start of the ind-th timespan - δt::TT = Δt-(tinterp.t[ind]) + ind::Int, δt::TT = getinterpindex(tinterp, t) # Evaluate tinterp.x[ind] at δt - return (tinterp.x[ind]::Taylor1{U})(δt) + return (tinterp.x[ind])(δt)::Taylor1{TT} end -function (tinterp::TaylorInterpolant{T, U, 2})(t::TT) where {T, U, TT<:Union{T, Taylor1{U}}} +function (tinterp::TaylorInterpolant{T, U, 2})(t::T) where {T, U} + # Get index of tinterp.x that interpolates at time t + ind::Int, δt::T = getinterpindex(tinterp, t) + # Evaluate tinterp.x[ind] at δt + return view(tinterp.x, ind, :)(δt)::Vector{U} +end +function (tinterp::TaylorInterpolant{T, U, 2})(t::TT) where {T, U, TT<:Union{Taylor1{U},Taylor1{TaylorN{U}}}} # Get index of tinterp.x that interpolates at time t - ind::Int, Δt::TT = getinterpindex(tinterp, t) - # Time since the start of the ind-th timespan - δt::TT = Δt-tinterp.t[ind] + ind::Int, δt::TT = getinterpindex(tinterp, t) # Evaluate tinterp.x[ind] at δt - return tinterp.x[ind,:](δt) + return view(tinterp.x, ind, :)(δt)::Vector{TT} end function (tinterp::TaylorInterpolant{T, U, 2})(target::Int, observer::Int, t::T) where {T, U} @@ -138,7 +149,7 @@ function (tinterp::TaylorInterpolant{T, U, 2})(target::Int, observer::Int, t::T) end end -(tinterp::TaylorInterpolant{T,U,2})(target::Int, t::TT) where {T, U, TT<:Union{T, Taylor1{U}}} = tinterp(target, 0, t) +(tinterp::TaylorInterpolant{T,U,2})(target::Int, t::TT) where {T, U, TT<:Union{T, Taylor1{U}, Taylor1{TaylorN{U}}}} = tinterp(target, 0, t) @doc raw""" reverse(tinterp::TaylorInterpolant{T,U,N}) where {T<:Real, U<:Number, N} diff --git a/src/jpl-de-430-431-earth-orientation-model.jl b/src/jpl-de-430-431-earth-orientation-model.jl index fd2db31..855fcf6 100644 --- a/src/jpl-de-430-431-earth-orientation-model.jl +++ b/src/jpl-de-430-431-earth-orientation-model.jl @@ -1,3 +1,5 @@ +# This file is part of the PlanetaryEphemeris.jl package; MIT licensed + # JPL-DE-430/431 model for the orientation of the Earth # From JPL DE430 documentation (Folkner et al., 2014): @@ -10,8 +12,8 @@ @doc raw""" Ω(t) -Return the longitude (in radians) of the mean ascending node of the lunar orbit on the -ecliptic, measured from the mean equinox of date +Return the longitude (in radians) of the mean ascending node of the lunar orbit on the +ecliptic, measured from the mean equinox of date ```math \Omega(t) = 125^\circ 02' 40''.280 - \left(1934^\circ 8' 10''.539\right) T + 7''.455 T^2 + 0''.008 T^3, ``` @@ -24,12 +26,12 @@ See equation (5-64) in page 5-27 of https://doi.org/10.1002/0471728470. @doc raw""" Delta_psi(t) -Return the nutation in longitude (in radians) +Return the nutation in longitude (in radians) ```math \Delta\psi(t) = -17''.1996 \sin\Omega, ``` where ``t = 36,525 T`` is the TDB time in Julian days from J2000.0 and ``\Omega`` is the -longitude of the mean ascending node of the lunar orbit. +longitude of the mean ascending node of the lunar orbit. See equation (5-190) in page 5-72 of https://doi.org/10.1002/0471728470. @@ -40,7 +42,7 @@ Delta_psi(t) = deg2rad( (-17.1996/3600)*sin(Ω(t)) ) @doc raw""" Delta_epsilon(t) -Return the nutation in obliquity (in radians) +Return the nutation in obliquity (in radians) ```math \Delta\epsilon(t) = 9''.2025 \cos\Omega, ``` @@ -78,7 +80,7 @@ end @doc raw""" ϵ̄(t) -Return the mean obliquity (in radians) +Return the mean obliquity (in radians) ```math \bar{\epsilon}(t) = 84,381''.448 - 46''.8150 T - 0''.00059 T^2 + 0''.001813 T^3, ``` @@ -90,14 +92,14 @@ See equation (5-153) in page 5-61 of https://doi.org/10.1002/0471728470. @doc raw""" pole_frame(t) - -Return the pole unit vector in the inertial frame ``\mathbf{p}_\mathrm{E}``, computed by + +Return the pole unit vector in the inertial frame ``\mathbf{p}_\mathrm{E}``, computed by precessing the pole of date with an estimated linear correction ```math \mathbf{p}_\mathrm{E} = R_z(\zeta_A)R_y(-\theta_A)R_z(z_A)R_x(-\phi_x)R_y(-\phi_y)\mathbf{p}_\mathrm{d}, ``` -where ``t`` is the TDB time in Julian days from J2000.0, ``\mathbf{p}_d`` is the true pole -of date unit vector, ``(\zeta_A, \theta_A, z_A)`` are the equatorial precession angles and +where ``t`` is the TDB time in Julian days from J2000.0, ``\mathbf{p}_d`` is the true pole +of date unit vector, ``(\zeta_A, \theta_A, z_A)`` are the equatorial precession angles and ``(\phi_x, \phi_y)`` are the linear corrections. See equation (25) in page 11 of https://ui.adsabs.harvard.edu/abs/2014IPNPR.196C...1F%2F/abstract. @@ -123,8 +125,8 @@ Return the X-axis linear correction to precession (in radians) ```math \phi_x = \phi_{x0} + 100T\times \frac{d\phi_x}{dt}, ``` -where ``t = 36,525 T`` is the TDB time in Julian days from J2000.0, ``\phi_{x0}`` is the -X–axis rotation at J2000.0 (in arcseconds) and ``\frac{d\phi_x}{dt}`` is the negative +where ``t = 36,525 T`` is the TDB time in Julian days from J2000.0, ``\phi_{x0}`` is the +X–axis rotation at J2000.0 (in arcseconds) and ``\frac{d\phi_x}{dt}`` is the negative obliquity rate correction (in arcseconds per year). See equation (25) in page 11 and Table 10 in page 50 of https://ui.adsabs.harvard.edu/abs/2014IPNPR.196C...1F%2F/abstract. @@ -143,8 +145,8 @@ Return the Y-axis linear correction to precession (in radians) ```math \phi_y = \phi_{y0} + 100T\times \frac{d\phi_y}{dt}, ``` -where ``t = 36,525 T`` is the TDB time in Julian days from J2000.0, ``\phi_{y0}`` is the -Y–axis rotation at J2000.0 (in arcseconds) and ``\frac{d\phi_y}{dt}`` is precession rate +where ``t = 36,525 T`` is the TDB time in Julian days from J2000.0, ``\phi_{y0}`` is the +Y–axis rotation at J2000.0 (in arcseconds) and ``\frac{d\phi_y}{dt}`` is precession rate correction times sine of obliquity (in arcseconds per year). See equation (25) in page 11 and Table 10 in page 50 of https://ui.adsabs.harvard.edu/abs/2014IPNPR.196C...1F%2F/abstract. @@ -158,7 +160,7 @@ const Dt_phi_y = -1.2118591216559240E-03 # Precession rate correction times si @doc raw""" Zeta(t) - + Return the ``\zeta_A`` equatorial precession angle (in radians) ```math \zeta_A(t) = 2306''.2181 T + 0''.30188 T^2 + 0''.017998 T^3, @@ -179,7 +181,7 @@ end @doc raw""" Theta(t) - + Return the ``\theta_A`` equatorial precession angle (in radians) ```math \theta_A(t) = 2004''.3109 T - 0''.42665 T^2 - 0''.041833 T^3, @@ -200,7 +202,7 @@ end @doc raw""" zeta(t) - + Return the ``z_A`` equatorial precession angle (in radians) ```math z_A(t) = 2306''.2181 T + 1''.09468 T^2 + 0''.018203 T^3, @@ -226,16 +228,16 @@ end Return the rotation matrix around the x-axis ```math -R_x(\alpha) = +R_x(\alpha) = \left[ \begin{array}{ccc} 1 & 0 & 0 \\ 0 & \cos\alpha & \sin\alpha \\ - 0 & -\sin\alpha & \cos\alpha \\ + 0 & -\sin\alpha & \cos\alpha \\ \end{array} \right], ``` -where ``\alpha`` is an angle in radians. +where ``\alpha`` is an angle in radians. See equation (11) in page 9 of https://ui.adsabs.harvard.edu/abs/2014IPNPR.196C...1F%2F/abstract. @@ -244,12 +246,12 @@ See also [`Ry`](@ref) and [`Rz`](@ref). function Rx(alpha::T) where {T<:Number} # Allocate memmory res = Array{T}(undef, 3, 3) - # Local variables + # Local variables one_alpha = one(alpha) zero_alpha = zero(alpha) cos_alpha = cos(alpha) sin_alpha = sin(alpha) - # Matrix elements + # Matrix elements res[1, 1] = one_alpha res[2, 1] = zero_alpha res[3, 1] = zero_alpha @@ -267,22 +269,22 @@ end Return the rotation matrix around the y-axis ```math -R_y(\alpha) = +R_y(\alpha) = \left[ \begin{array}{ccc} \cos\alpha & 0 & -\sin\alpha \\ 0 & 1 & 0 \\ - \sin\alpha & 0 & \cos\alpha \\ + \sin\alpha & 0 & \cos\alpha \\ \end{array} \right], ``` -where ``\alpha`` is an angle in radians. +where ``\alpha`` is an angle in radians. See equation (12) in page 9 of https://ui.adsabs.harvard.edu/abs/2014IPNPR.196C...1F%2F/abstract. -``\textbf{Note:}`` Lieske et al. (1977, 1979) introduces the convention that this matrix -has opposite signs outside the sine terms (with respect to `Rx`, `Rz`). See equation (4) -in page 3 of https://ui.adsabs.harvard.edu/abs/1977A%26A....58....1L/abstract and equation +``\textbf{Note:}`` Lieske et al. (1977, 1979) introduces the convention that this matrix +has opposite signs outside the sine terms (with respect to `Rx`, `Rz`). See equation (4) +in page 3 of https://ui.adsabs.harvard.edu/abs/1977A%26A....58....1L/abstract and equation (3) in page 283 of https://ui.adsabs.harvard.edu/abs/1979A%26A....73..282L/abstract. See also [`Rx`](@ref) and [`Rz`](@ref). @@ -313,16 +315,16 @@ end Return the rotation matrix around the z-axis ```math -R_z(\alpha) = +R_z(\alpha) = \left[ \begin{array}{ccc} \cos\alpha & \sin\alpha & 0 \\ -\sin\alpha & \cos\alpha & 0 \\ - 0 & 0 & 1 \\ + 0 & 0 & 1 \\ \end{array} \right], ``` -where ``\alpha`` is an angle in radians. +where ``\alpha`` is an angle in radians. See equation (13) in page 9 of https://ui.adsabs.harvard.edu/abs/2014IPNPR.196C...1F%2F/abstract. @@ -408,7 +410,7 @@ end @doc raw""" pole_rotation(α::T, δ::T, W::T=zero(α)) where {T <: Number} -Return the rotation matrix from the inertial frame to the frame with pole at right +Return the rotation matrix from the inertial frame to the frame with pole at right ascension ``\alpha``, declination ``\delta`` and prime meridian at ``W`` ```math A = R_z(W + \Delta W)R_x\left(\frac{\pi}{2} - \delta - \Delta\delta\right)R_z\left(\alpha + \Delta\alpha + \frac{\pi}{2}\right). @@ -416,7 +418,7 @@ A = R_z(W + \Delta W)R_x\left(\frac{\pi}{2} - \delta - \Delta\delta\right)R_z\le See equation (6-3) in page 6-5 of https://doi.org/10.1002/0471728470. -``\textbf{Note:}`` Rotation matrices ``(R_x, R_y, R_z)`` convention is the same as +``\textbf{Note:}`` Rotation matrices ``(R_x, R_y, R_z)`` convention is the same as equations (11), (12) and (13) in page 9 of https://ui.adsabs.harvard.edu/abs/2014IPNPR.196C...1F%2F/abstract. See also [`Rx`](@ref) and [`Rz`](@ref). @@ -446,9 +448,9 @@ Return the IAU 1980 nutation matrix (Explanatory Supplement to the Astronomical ```math N = R_x(-\bar{\epsilon}(t) - \Delta\epsilon(t))R_z(-\Delta\psi(t))R_x(\bar{\epsilon}(t)), ``` -where ``t`` is the TDB time in Julian days from J2000.0, ``\bar{\epsilon}`` is the mean -obliquity, ``\Delta \psi`` is the nutation in longitude and ``\Delta\epsilon`` is the -nutation in obliquity. All angles are in radians. +where ``t`` is the TDB time in Julian days from J2000.0, ``\bar{\epsilon}`` is the mean +obliquity, ``\Delta \psi`` is the nutation in longitude and ``\Delta\epsilon`` is the +nutation in obliquity. All angles are in radians. See equation (5-152) in page 5-60 of https://doi.org/10.1002/0471728470. @@ -486,7 +488,7 @@ end @doc raw""" c2t_jpl_de430(t) -Return the matrix for celestial-to-terrestrial coordinate transformation, according to +Return the matrix for celestial-to-terrestrial coordinate transformation, according to JPL DE 430/431 Earth orientation model ```math \texttt{c2t_jpl_de430}(t) = A\times C\times N, @@ -512,8 +514,8 @@ end @doc raw""" moon_omega(ϕ::Taylor1, θ::Taylor1, ψ::Taylor1) - -Return the Moon's angular velocity, computed by differentiating the Euler angles + +Return the Moon's angular velocity, computed by differentiating the Euler angles ``(\phi, \theta, \psi)`` ```math \begin{align*} @@ -547,7 +549,7 @@ Return the first term of the time-dependent part of lunar total moment of intert \end{array} \right], ``` -where ``\mathbf{r} = (x, y, z)`` is the position of the Moon relative to Earth referred to +where ``\mathbf{r} = (x, y, z)`` is the position of the Moon relative to Earth referred to the mantle frame, ``k_{2,M}`` is the lunar potential Love number, ``m_E`` is the mass of the Earth and ``R_M`` is the equatorial radius of the Moon. @@ -588,8 +590,8 @@ Return the second term of the time-dependent part of lunar total moment of inter \end{array} \right], ``` -where ``\mathbf{\omega}_m = (\omega_{m,x}, \omega_{m,y}, \omega_{m,z})`` is the angular -velocity of the mantle in the mantle frame, ``k_{2,M}`` is the lunar potential Love number, +where ``\mathbf{\omega}_m = (\omega_{m,x}, \omega_{m,y}, \omega_{m,z})`` is the angular +velocity of the mantle in the mantle frame, ``k_{2,M}`` is the lunar potential Love number, ``R_M`` is the equatorial radius of the Moon and ``n`` is the lunar mean motion. ``\textbf{Note:}`` Euler angles must be evaluated at time ``t-τ_M``. @@ -601,7 +603,7 @@ See also [`ITM1`](@ref) and [`ITM`](@ref). function ITM2(ωx::T, ωy::T, ωz::T) where {T <: Number} # Compute corresponding matrix elements m = Matrix{T}(undef, 3, 3) - ωx2 = ωx*ωx + ωx2 = ωx*ωx ωy2 = ωy*ωy ωz2 = ωz*ωz # Angular speed @@ -626,7 +628,7 @@ end Return lunar mantle inertia tensor ```math -\mathbf{I}_m(t) = \tilde{\mathbf{I}}_m +\mathbf{I}_m(t) = \tilde{\mathbf{I}}_m - \frac{k_{2,M} m_E R_M^5}{r^5} \left[ diff --git a/src/osculating.jl b/src/osculating.jl index f538d55..00fdca6 100644 --- a/src/osculating.jl +++ b/src/osculating.jl @@ -1,41 +1,43 @@ +# This file is part of the PlanetaryEphemeris.jl package; MIT licensed + @doc raw""" semimajoraxis(x, y, z, u, v, w, m1, m2) -Calculate semimajor axis for the two body problem defined by the relative position +Calculate semimajor axis for the two body problem defined by the relative position `(x,y,z)` and velocity `(u,v,w)` vectors between two bodies with masses `m1` and `m2`. """ function semimajoraxis(x, y, z, u, v, w, m1, m2) r = sqrt(x^2+y^2+z^2) # Position magnitude vsq = u^2+v^2+w^2 # Velocity magnitude squared M = m1+m2 # Total mass - return 1/( (2/r)-(vsq/M) ) # Semimajor axis + return 1/( (2/r)-(vsq/M) ) # Semimajor axis end @doc raw""" eccentricity(x, y, z, u, v, w, m1, m2) -Calculate eccentricity for the two body problem defined by the relative position +Calculate eccentricity for the two body problem defined by the relative position `(x,y,z)` and velocity `(u,v,w)` vectors between two bodies with masses `m1` and `m2`. """ function eccentricity(x, y, z, u, v, w, m1, m2) # h: Angular momentum per unit mass hsq = (y*w-z*v)^2 + (z*u-x*w)^2 + (x*v-y*u)^2 # h magnitude squared a = semimajoraxis(x, y, z, u, v, w, m1, m2) # Semimajor axis - M = m1+m2 # Total mass - quotient = hsq/( M*a ) - return sqrt(1 - quotient) # Eccentricity + M = m1+m2 # Total mass + quotient = hsq/( M*a ) + return sqrt(1 - quotient) # Eccentricity end @doc raw""" inclination(x, y, z, u, v, w) -Calculate inclination for the two body problem defined by the relative position +Calculate inclination for the two body problem defined by the relative position `(x,y,z)` and velocity `(u,v,w)` vectors between two bodies.""" function inclination(x, y, z, u, v, w) # h: Angular momentum per unit mass hz = x*v-y*u # z-coordinate of h h = sqrt( (y*w-z*v)^2 + (z*u-x*w)^2 + (x*v-y*u)^2 ) # h magnitude squared - return acos(hz/h) # Inclination + return acos(hz/h) # Inclination end @doc raw""" @@ -50,15 +52,15 @@ See also [`semimajoraxis`](@ref), [`eccentricity`](@ref) and [`inclination`](@re function aei(x, y, z, u, v, w, m1, m2) r = sqrt(x^2+y^2+z^2) # Position magnitude vsq = u^2+v^2+w^2 # Velocity magnitude squared - M = m1+m2 # Total mass + M = m1+m2 # Total mass a = 1/( (2/r)-(vsq/M) ) # Semimajor axis # h: Angular momentum per unit mass hsq = (y*w-z*v)^2 + (z*u-x*w)^2 + (x*v-y*u)^2 # h magnitude squared h = sqrt(hsq) # h magnitude - quotient = hsq/( M*a ) + quotient = hsq/( M*a ) e = sqrt(1 - quotient) # Eccentricity hz = x*v-y*u # z-coordinate of h - inc = acos(hz/h) # Inclination + inc = acos(hz/h) # Inclination return a, e, inc end @@ -76,10 +78,10 @@ function ae(x, y, z, u, v, w, m1, m2) r = sqrt(x^2+y^2+z^2) # Position magnitude vsq = u^2+v^2+w^2 # Velocity magnitude squared M = m1+m2 # Total mass - a = 1/( (2/r)-(vsq/M) ) # Semimajor axis + a = 1/( (2/r)-(vsq/M) ) # Semimajor axis # h: Angular momentum per unit mass hsq = (y*w-z*v)^2 + (z*u-x*w)^2 + (x*v-y*u)^2 # h magnitude squared - quotient = hsq/( M*a ) + quotient = hsq/( M*a ) e = sqrt(1 - quotient) # Eccentricity return a, e end @@ -87,12 +89,12 @@ end @doc raw""" longascnode(x, y, z, u, v, w) -Calculate longitude of ascending node for the two body problem defined by the relative +Calculate longitude of ascending node for the two body problem defined by the relative position `(x,y,z)` and velocity `(u,v,w)` vectors between two bodies. """ function longascnode(x, y, z, u, v, w) - # The longitude of ascending node is computed as the angle between x-axis and the + # The longitude of ascending node is computed as the angle between x-axis and the # vector n = (-hy, hx, 0) where hx, hy, are resp., the x and y comps. of angular # momentum vector per unit mass, h @@ -110,7 +112,7 @@ end argperi(x, y, z, u, v, w, m1, m2) Calculate the instantaneous (osculating) argument of pericentre for the two -body problem defined by the relative position `(x,y,z)` and velocity `(u,v,w)` vectors +body problem defined by the relative position `(x,y,z)` and velocity `(u,v,w)` vectors between two bodies with masses `m1` and `m2`. """ function argperi(x, y, z, u, v, w, m1, m2) @@ -134,7 +136,7 @@ end longperi(x, y, z, u, v, w, m1, m2) Calculate the instantaneous (osculating) longitude of pericentre for the two -body problem defined by the relative position `(x,y,z)` and velocity `(u,v,w)` vectors +body problem defined by the relative position `(x,y,z)` and velocity `(u,v,w)` vectors between two bodies with masses `m1` and `m2`. """ function longperi(x, y, z, u, v, w, m1, m2) @@ -146,20 +148,20 @@ end trueanomaly(x, y, z, u, v, w, m1, m2) Calculate the instantaneous (osculating) true anomaly for the two -body problem defined by the relative position `(x,y,z)` and velocity `(u,v,w)` vectors +body problem defined by the relative position `(x,y,z)` and velocity `(u,v,w)` vectors between two bodies with masses `m1` and `m2`. """ function trueanomaly(x, y, z, u, v, w, m1, m2) R = sqrt(x^2+y^2+z^2) # Position magnitude V2 = u^2+v^2+w^2 # Velocity mangitude squared - V = sqrt(V2) # Velocity magnitude + V = sqrt(V2) # Velocity magnitude a = 1.0/( (2.0/R)-(V2/(m1+m2)) ) # Semimajor axis # h: Angular momentum per unit mass - h2 = (y*w-z*v)^2 + (z*u-x*w)^2 + (x*v-y*u)^2 # h magnitude squared - h = sqrt(h2) # h magnitude + h2 = (y*w-z*v)^2 + (z*u-x*w)^2 + (x*v-y*u)^2 # h magnitude squared + h = sqrt(h2) # h magnitude quotient=h2/( (m1+m2)*a ) @@ -183,7 +185,7 @@ end ecanomaly(x, y, z, u, v, w, m1, m2) Calculate the instantaneous (osculating) eccentric anomaly for the two -body problem defined by the relative position `(x,y,z)` and velocity `(u,v,w)` vectors +body problem defined by the relative position `(x,y,z)` and velocity `(u,v,w)` vectors between two bodies with masses `m1` and `m2`. See also [`semimajoraxis`](@ref) and [`eccentricity`](@ref). @@ -199,7 +201,7 @@ function ecanomaly(x, y, z, u, v, w, m1, m2) cosE = (a-R)/(a*e) # sin( eccentric anomaly ) sinE = Rdot*R/( e*sqrt(mu*a) ) # cos( eccentric anomaly ) - res = atan(sinE, cosE) # Eccentric anomaly + res = atan(sinE, cosE) # Eccentric anomaly return res @@ -208,8 +210,8 @@ end @doc """ meananomaly(x, y, z, u, v, w, m1, m2) -Calculate the instantaneous (osculating) mean anomaly for the two -body problem defined by the relative position `(x,y,z)` and velocity `(u,v,w)` vectors +Calculate the instantaneous (osculating) mean anomaly for the two +body problem defined by the relative position `(x,y,z)` and velocity `(u,v,w)` vectors between two bodies with masses `m1` and `m2`. See also [`ecanomaly`](@ref) and [`eccentricity`](@ref). @@ -235,7 +237,7 @@ function timeperipass(t, x, y, z, u, v, w, m1, m2) me_an = meananomaly(x, y, z, u, v, w, m1, m2) # Mean anomaly a = semimajoraxis(x, y, z, u, v, w, m1, m2) # Semimajor axis - me_mo = meanmotion(m1+m2, a) # Mean motion + me_mo = meanmotion(m1+m2, a) # Mean motion # Time of pericentre passage return t-me_an/me_mo end @@ -244,7 +246,7 @@ end rungelenzx(x, y, z, u, v, w, m1, m2) Calculate the instantaneous (osculating) x-component of the Laplace-Runge-Lenz vector for the -two body problem defined by the relative position `(x,y,z)` and velocity `(u,v,w)` vectors +two body problem defined by the relative position `(x,y,z)` and velocity `(u,v,w)` vectors between two bodies with masses `m1` and `m2`. See also [`rungelenzy`](@ref), [`rungelenzz`](@ref), [`rungelenzmag`](@ref) and [`lrlvec`](@ref). @@ -254,7 +256,7 @@ function rungelenzx(x, y, z, u, v, w, m1, m2) # Velocity v = [u, v, w] # Angular momentum per unit mass h = [y.*w-z.*v, z.*u-x.*w, x.*v-y.*u] # Total mass - mu = m1+m2 + mu = m1+m2 # Position r = [x, y, z] rmag = sqrt(x^2+y^2+z^2) # Position magnitude # x-component of Laplace-Runge-Lenz vector @@ -267,7 +269,7 @@ end rungelenzy(x, y, z, u, v, w, m1, m2) Calculate the instantaneous (osculating) y-component of the Laplace-Runge-Lenz vector for the -two body problem defined by the relative position `(x,y,z)` and velocity `(u,v,w)` vectors +two body problem defined by the relative position `(x,y,z)` and velocity `(u,v,w)` vectors between two bodies with masses `m1` and `m2`. See also [`rungelenzx`](@ref), [`rungelenzz`](@ref), [`rungelenzmag`](@ref) and [`lrlvec`](@ref). @@ -290,7 +292,7 @@ end rungelenzz(x, y, z, u, v, w, m1, m2) Calculate the instantaneous (osculating) z-component of the Laplace-Runge-Lenz vector for the -two body problem defined by the relative position `(x,y,z)` and velocity `(u,v,w)` vectors +two body problem defined by the relative position `(x,y,z)` and velocity `(u,v,w)` vectors between two bodies with masses `m1` and `m2`. See also [`rungelenzx`](@ref), [`rungelenzy`](@ref), [`rungelenzmag`](@ref) and [`lrlvec`](@ref). @@ -313,7 +315,7 @@ end rungelenzmag(x, y, z, u, v, w, m1, m2) Calculate the instantaneous (osculating) magnitude of the Laplace-Runge-Lenz vector for the -two body problem defined by the relative position `(x,y,z)` and velocity `(u,v,w)` vectors +two body problem defined by the relative position `(x,y,z)` and velocity `(u,v,w)` vectors between two bodies with masses `m1` and `m2`. See also [`rungelenzx`](@ref), [`rungelenzy`](@ref), [`rungelenzz`](@ref) and [`lrlvec`](@ref). @@ -326,7 +328,7 @@ function rungelenzmag(x, y, z, u, v, w, m1, m2) mu = m1+m2 # Position r = [x, y, z] rmag = sqrt(x^2+y^2+z^2) # Position magnitude - + # Components of Laplace-Runge-Lenz vector lrl_x = ( -(z*u-x*w)*w+(x*v-y*u)*v )/(mu)-x/rmag # x-coordinate lrl_y = ( -(x*v-y*u)*u+(y*w-z*v)*w )/(mu)-y/rmag # y-coordinate @@ -339,7 +341,7 @@ end lrlvec(x, y, z, u, v, w, m1, m2) Calculate the instantaneous (osculating) cartesian components of the Laplace-Runge-Lenz vector -for the two body problem defined by the relative position `(x,y,z)` and velocity `(u,v,w)` +for the two body problem defined by the relative position `(x,y,z)` and velocity `(u,v,w)` vectors between two bodies with masses `m1` and `m2`. See also [`rungelenzx`](@ref), [`rungelenzy`](@ref), [`rungelenzz`](@ref) and [`rungelenzmag`](@ref). @@ -360,7 +362,7 @@ end @doc raw""" eccentricanomaly(e::T, M::T) where {T <: Number} -Compute eccentric anomaly (`E`) from eccentricity (`e`) and mean anomaly (`M`). +Compute eccentric anomaly (`E`) from eccentricity (`e`) and mean anomaly (`M`). See pages 32-37 of https://doi.org/10.1017/CBO9781139174817. """ function eccentricanomaly(e::T, M::T) where {T <: Number} @@ -417,7 +419,7 @@ end @doc raw""" time2truean(a, e, mu, t, taup) -Compute true anomaly from time (`t`), semimajor axis (`a`), eccentricity (`e`), +Compute true anomaly from time (`t`), semimajor axis (`a`), eccentricity (`e`), mass parameter (`mu`) and time of pericenter passage (`taup`). See also [`meanan2truean`](@ref), [`meananomaly`](@ref) and [`meanmotion`](@ref). diff --git a/src/propagation.jl b/src/propagation.jl index 52557ac..73ef53c 100644 --- a/src/propagation.jl +++ b/src/propagation.jl @@ -1,3 +1,5 @@ +# This file is part of the PlanetaryEphemeris.jl package; MIT licensed + @doc raw""" loadeph(ss16asteph::TaylorInterpolant, μ::Vector{<:Real}) diff --git a/test/Project.toml b/test/Project.toml index 403e0c3..1b368b0 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -5,9 +5,11 @@ JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Quadmath = "be4d8f0f-7fa4-5f49-b795-2f01399ab2dd" SPICE = "5bab7191-041a-5c2e-a744-024b9c3a5062" +TaylorSeries = "6aa5eb33-94cf-58f4-a9d0-e4b2c4fc25ea" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [compat] JLD2 = "0.4" Quadmath = "0.5" SPICE = "0.2" +TaylorSeries = "0.15" diff --git a/test/propagation.jl b/test/propagation.jl index aca9986..ae51a21 100644 --- a/test/propagation.jl +++ b/test/propagation.jl @@ -4,6 +4,7 @@ using PlanetaryEphemeris using Dates using Quadmath using JLD2 +using TaylorSeries using Test using PlanetaryEphemeris: initialcond, ssic_1969, ssic_2000, astic_1969, astic_2000 @@ -13,7 +14,6 @@ using SPICE: furnsh, spkgeo using PlanetaryEphemeris: order, abstol using LinearAlgebra: norm - @testset "Propagation" begin # Number of asteroids @@ -68,84 +68,89 @@ using LinearAlgebra: norm @test sol.t0 == 0.0 @test length(sol.t) == size(sol.x, 1) + 1 @test length(q0) == size(sol.x, 2) + dq = TaylorSeries.set_variables("dq", order=2, numvars=2) + tmid = sol.t0 + sol.t[1]/2 + @test sol(tmid) isa Vector{Float64} + @test sol(tmid + Taylor1(order)) isa Vector{Taylor1{Float64}} + @test sol(tmid + Taylor1([dq[1],dq[1]*dq[2]],order)) isa Vector{Taylor1{TaylorN{Float64}}} end - @testset "Propagation: DE430 dynamical model" begin - - @show Threads.nthreads() - - # Float64 - - # Test integration - sol64 = propagate(100, jd0, nyears, dense; dynamics, order, abstol) - # Save solution - filename = selecteph2jld2(sol64, bodyind, nyears) - # Recovered solution - recovered_sol64 = JLD2.load(filename, "ss16ast_eph") - - @test selecteph(sol64, bodyind, euler = true, ttmtdb = true) == recovered_sol64 - - LSK = "https://naif.jpl.nasa.gov/pub/naif/generic_kernels/lsk/naif0012.tls" - TTmTDBK = "https://ssd.jpl.nasa.gov/ftp/eph/planets/bsp/TTmTDB.de430.19feb2015.bsp" - SPK = "https://ssd.jpl.nasa.gov/ftp/eph/planets/bsp/de430_1850-2150.bsp" - - # Download kernels - Downloads.download(LSK, "naif0012.tls") - Downloads.download(SPK, "de430_1850-2150.bsp") - Downloads.download(TTmTDBK, "TTmTDB.de430.19feb2015.bsp") - - # Load kernels - furnsh("naif0012.tls", "de430_1850-2150.bsp", "TTmTDB.de430.19feb2015.bsp") - - ttmtdb_pe = TaylorInterpolant(sol64.t0, sol64.t, sol64.x[:, 6N+13]) # TT-TDB - posvel_pe_su = selecteph(sol64,su) # Sun - posvel_pe_ea = selecteph(sol64,ea) # Earth - posvel_pe_mo = selecteph(sol64,mo) # Moon - posvel_pe_ma = selecteph(sol64,6) # Mars - posvel_pe_ju = selecteph(sol64,7) # Jupiter - - ttmtdb_jpl(et) = spkgeo(1000000001, et, "J2000", 1000000000)[1][1] # TT-TDB - posvel_jpl_su(et) = kmsec2auday(spkgeo(10, et, "J2000", 0)[1]) # Sun - posvel_jpl_ea(et) = kmsec2auday(spkgeo(399, et, "J2000", 0)[1]) # Earth - posvel_jpl_mo(et) = kmsec2auday(spkgeo(301, et, "J2000", 0)[1]) # Moon - posvel_jpl_ma(et) = kmsec2auday(spkgeo(4, et, "J2000", 0)[1]) # Mars - posvel_jpl_ju(et) = kmsec2auday(spkgeo(5, et, "J2000", 0)[1]) # Jupiter - - tv = range(sol64.t0, sol64.t[end], 10) - for t in tv - et = t * daysec - @show t, et - @show abs(ttmtdb_jpl(et) - ttmtdb_pe(t)) - @show norm(posvel_jpl_su(et) - posvel_pe_su(t), Inf) - @show norm(posvel_jpl_ea(et) - posvel_pe_ea(t), Inf) - @show norm(posvel_jpl_mo(et) - posvel_pe_mo(t), Inf) - @show norm(posvel_jpl_ma(et) - posvel_pe_ma(t), Inf) - @show norm(posvel_jpl_ju(et) - posvel_pe_ju(t), Inf) - - @test abs(ttmtdb_jpl(et) - ttmtdb_pe(t)) < 1e-12 - @test norm(posvel_jpl_su(et) - posvel_pe_su(t), Inf) < 1e-14 - @test norm(posvel_jpl_ea(et) - posvel_pe_ea(t), Inf) < 1e-11 - @test norm(posvel_jpl_mo(et) - posvel_pe_mo(t), Inf) < 1e-11 - @test norm(posvel_jpl_ma(et) - posvel_pe_ma(t), Inf) < 1e-12 - @test norm(posvel_jpl_ju(et) - posvel_pe_ju(t), Inf) < 1e-13 - end - - # Remove files - rm.((filename, "naif0012.tls", "de430_1850-2150.bsp", "TTmTDB.de430.19feb2015.bsp")) - - # Float 128 - #= - # Test integration - sol128 = propagate(1, Float128(jd0), nyears, dense; dynamics = dynamics, order = order, abstol = abstol) - # Save solution - filename = selecteph2jld2(sol128, bodyind, nyears) - # Recovered solution - recovered_sol128 = JLD2.load(filename, "ss16ast_eph") - # Remove file - rm(filename) - - @test selecteph(sol128, bodyind, euler = true, ttmtdb = true) == recovered_sol128 - =# - end + # @testset "Propagation: DE430 dynamical model" begin + + # @show Threads.nthreads() + + # # Float64 + + # # Test integration + # sol64 = propagate(100, jd0, nyears, dense; dynamics, order, abstol) + # # Save solution + # filename = selecteph2jld2(sol64, bodyind, nyears) + # # Recovered solution + # recovered_sol64 = JLD2.load(filename, "ss16ast_eph") + + # @test selecteph(sol64, bodyind, euler = true, ttmtdb = true) == recovered_sol64 + + # LSK = "https://naif.jpl.nasa.gov/pub/naif/generic_kernels/lsk/naif0012.tls" + # TTmTDBK = "https://ssd.jpl.nasa.gov/ftp/eph/planets/bsp/TTmTDB.de430.19feb2015.bsp" + # SPK = "https://ssd.jpl.nasa.gov/ftp/eph/planets/bsp/de430_1850-2150.bsp" + + # # Download kernels + # Downloads.download(LSK, "naif0012.tls") + # Downloads.download(SPK, "de430_1850-2150.bsp") + # Downloads.download(TTmTDBK, "TTmTDB.de430.19feb2015.bsp") + + # # Load kernels + # furnsh("naif0012.tls", "de430_1850-2150.bsp", "TTmTDB.de430.19feb2015.bsp") + + # ttmtdb_pe = TaylorInterpolant(sol64.t0, sol64.t, sol64.x[:, 6N+13]) # TT-TDB + # posvel_pe_su = selecteph(sol64,su) # Sun + # posvel_pe_ea = selecteph(sol64,ea) # Earth + # posvel_pe_mo = selecteph(sol64,mo) # Moon + # posvel_pe_ma = selecteph(sol64,6) # Mars + # posvel_pe_ju = selecteph(sol64,7) # Jupiter + + # ttmtdb_jpl(et) = spkgeo(1000000001, et, "J2000", 1000000000)[1][1] # TT-TDB + # posvel_jpl_su(et) = kmsec2auday(spkgeo(10, et, "J2000", 0)[1]) # Sun + # posvel_jpl_ea(et) = kmsec2auday(spkgeo(399, et, "J2000", 0)[1]) # Earth + # posvel_jpl_mo(et) = kmsec2auday(spkgeo(301, et, "J2000", 0)[1]) # Moon + # posvel_jpl_ma(et) = kmsec2auday(spkgeo(4, et, "J2000", 0)[1]) # Mars + # posvel_jpl_ju(et) = kmsec2auday(spkgeo(5, et, "J2000", 0)[1]) # Jupiter + + # tv = range(sol64.t0, sol64.t[end], 10) + # for t in tv + # et = t * daysec + # @show t, et + # @show abs(ttmtdb_jpl(et) - ttmtdb_pe(t)) + # @show norm(posvel_jpl_su(et) - posvel_pe_su(t), Inf) + # @show norm(posvel_jpl_ea(et) - posvel_pe_ea(t), Inf) + # @show norm(posvel_jpl_mo(et) - posvel_pe_mo(t), Inf) + # @show norm(posvel_jpl_ma(et) - posvel_pe_ma(t), Inf) + # @show norm(posvel_jpl_ju(et) - posvel_pe_ju(t), Inf) + + # @test abs(ttmtdb_jpl(et) - ttmtdb_pe(t)) < 1e-12 + # @test norm(posvel_jpl_su(et) - posvel_pe_su(t), Inf) < 1e-14 + # @test norm(posvel_jpl_ea(et) - posvel_pe_ea(t), Inf) < 1e-11 + # @test norm(posvel_jpl_mo(et) - posvel_pe_mo(t), Inf) < 1e-11 + # @test norm(posvel_jpl_ma(et) - posvel_pe_ma(t), Inf) < 1e-12 + # @test norm(posvel_jpl_ju(et) - posvel_pe_ju(t), Inf) < 1e-13 + # end + + # # Remove files + # rm.((filename, "naif0012.tls", "de430_1850-2150.bsp", "TTmTDB.de430.19feb2015.bsp")) + + # # Float 128 + # #= + # # Test integration + # sol128 = propagate(1, Float128(jd0), nyears, dense; dynamics = dynamics, order = order, abstol = abstol) + # # Save solution + # filename = selecteph2jld2(sol128, bodyind, nyears) + # # Recovered solution + # recovered_sol128 = JLD2.load(filename, "ss16ast_eph") + # # Remove file + # rm(filename) + + # @test selecteph(sol128, bodyind, euler = true, ttmtdb = true) == recovered_sol128 + # =# + # end end From 38d2febce50b5a207a9b4bda97cd5c2cf952a5d0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jorge=20A=2E=20P=C3=A9rez-Hern=C3=A1ndez?= Date: Sat, 27 Jan 2024 22:25:50 +0100 Subject: [PATCH 05/10] Uncomment tests --- test/propagation.jl | 154 ++++++++++++++++++++++---------------------- 1 file changed, 77 insertions(+), 77 deletions(-) diff --git a/test/propagation.jl b/test/propagation.jl index ae51a21..28160ba 100644 --- a/test/propagation.jl +++ b/test/propagation.jl @@ -75,82 +75,82 @@ using LinearAlgebra: norm @test sol(tmid + Taylor1([dq[1],dq[1]*dq[2]],order)) isa Vector{Taylor1{TaylorN{Float64}}} end - # @testset "Propagation: DE430 dynamical model" begin - - # @show Threads.nthreads() - - # # Float64 - - # # Test integration - # sol64 = propagate(100, jd0, nyears, dense; dynamics, order, abstol) - # # Save solution - # filename = selecteph2jld2(sol64, bodyind, nyears) - # # Recovered solution - # recovered_sol64 = JLD2.load(filename, "ss16ast_eph") - - # @test selecteph(sol64, bodyind, euler = true, ttmtdb = true) == recovered_sol64 - - # LSK = "https://naif.jpl.nasa.gov/pub/naif/generic_kernels/lsk/naif0012.tls" - # TTmTDBK = "https://ssd.jpl.nasa.gov/ftp/eph/planets/bsp/TTmTDB.de430.19feb2015.bsp" - # SPK = "https://ssd.jpl.nasa.gov/ftp/eph/planets/bsp/de430_1850-2150.bsp" - - # # Download kernels - # Downloads.download(LSK, "naif0012.tls") - # Downloads.download(SPK, "de430_1850-2150.bsp") - # Downloads.download(TTmTDBK, "TTmTDB.de430.19feb2015.bsp") - - # # Load kernels - # furnsh("naif0012.tls", "de430_1850-2150.bsp", "TTmTDB.de430.19feb2015.bsp") - - # ttmtdb_pe = TaylorInterpolant(sol64.t0, sol64.t, sol64.x[:, 6N+13]) # TT-TDB - # posvel_pe_su = selecteph(sol64,su) # Sun - # posvel_pe_ea = selecteph(sol64,ea) # Earth - # posvel_pe_mo = selecteph(sol64,mo) # Moon - # posvel_pe_ma = selecteph(sol64,6) # Mars - # posvel_pe_ju = selecteph(sol64,7) # Jupiter - - # ttmtdb_jpl(et) = spkgeo(1000000001, et, "J2000", 1000000000)[1][1] # TT-TDB - # posvel_jpl_su(et) = kmsec2auday(spkgeo(10, et, "J2000", 0)[1]) # Sun - # posvel_jpl_ea(et) = kmsec2auday(spkgeo(399, et, "J2000", 0)[1]) # Earth - # posvel_jpl_mo(et) = kmsec2auday(spkgeo(301, et, "J2000", 0)[1]) # Moon - # posvel_jpl_ma(et) = kmsec2auday(spkgeo(4, et, "J2000", 0)[1]) # Mars - # posvel_jpl_ju(et) = kmsec2auday(spkgeo(5, et, "J2000", 0)[1]) # Jupiter - - # tv = range(sol64.t0, sol64.t[end], 10) - # for t in tv - # et = t * daysec - # @show t, et - # @show abs(ttmtdb_jpl(et) - ttmtdb_pe(t)) - # @show norm(posvel_jpl_su(et) - posvel_pe_su(t), Inf) - # @show norm(posvel_jpl_ea(et) - posvel_pe_ea(t), Inf) - # @show norm(posvel_jpl_mo(et) - posvel_pe_mo(t), Inf) - # @show norm(posvel_jpl_ma(et) - posvel_pe_ma(t), Inf) - # @show norm(posvel_jpl_ju(et) - posvel_pe_ju(t), Inf) - - # @test abs(ttmtdb_jpl(et) - ttmtdb_pe(t)) < 1e-12 - # @test norm(posvel_jpl_su(et) - posvel_pe_su(t), Inf) < 1e-14 - # @test norm(posvel_jpl_ea(et) - posvel_pe_ea(t), Inf) < 1e-11 - # @test norm(posvel_jpl_mo(et) - posvel_pe_mo(t), Inf) < 1e-11 - # @test norm(posvel_jpl_ma(et) - posvel_pe_ma(t), Inf) < 1e-12 - # @test norm(posvel_jpl_ju(et) - posvel_pe_ju(t), Inf) < 1e-13 - # end - - # # Remove files - # rm.((filename, "naif0012.tls", "de430_1850-2150.bsp", "TTmTDB.de430.19feb2015.bsp")) - - # # Float 128 - # #= - # # Test integration - # sol128 = propagate(1, Float128(jd0), nyears, dense; dynamics = dynamics, order = order, abstol = abstol) - # # Save solution - # filename = selecteph2jld2(sol128, bodyind, nyears) - # # Recovered solution - # recovered_sol128 = JLD2.load(filename, "ss16ast_eph") - # # Remove file - # rm(filename) - - # @test selecteph(sol128, bodyind, euler = true, ttmtdb = true) == recovered_sol128 - # =# - # end + @testset "Propagation: DE430 dynamical model" begin + + @show Threads.nthreads() + + # Float64 + + # Test integration + sol64 = propagate(100, jd0, nyears, dense; dynamics, order, abstol) + # Save solution + filename = selecteph2jld2(sol64, bodyind, nyears) + # Recovered solution + recovered_sol64 = JLD2.load(filename, "ss16ast_eph") + + @test selecteph(sol64, bodyind, euler = true, ttmtdb = true) == recovered_sol64 + + LSK = "https://naif.jpl.nasa.gov/pub/naif/generic_kernels/lsk/naif0012.tls" + TTmTDBK = "https://ssd.jpl.nasa.gov/ftp/eph/planets/bsp/TTmTDB.de430.19feb2015.bsp" + SPK = "https://ssd.jpl.nasa.gov/ftp/eph/planets/bsp/de430_1850-2150.bsp" + + # Download kernels + Downloads.download(LSK, "naif0012.tls") + Downloads.download(SPK, "de430_1850-2150.bsp") + Downloads.download(TTmTDBK, "TTmTDB.de430.19feb2015.bsp") + + # Load kernels + furnsh("naif0012.tls", "de430_1850-2150.bsp", "TTmTDB.de430.19feb2015.bsp") + + ttmtdb_pe = TaylorInterpolant(sol64.t0, sol64.t, sol64.x[:, 6N+13]) # TT-TDB + posvel_pe_su = selecteph(sol64,su) # Sun + posvel_pe_ea = selecteph(sol64,ea) # Earth + posvel_pe_mo = selecteph(sol64,mo) # Moon + posvel_pe_ma = selecteph(sol64,6) # Mars + posvel_pe_ju = selecteph(sol64,7) # Jupiter + + ttmtdb_jpl(et) = spkgeo(1000000001, et, "J2000", 1000000000)[1][1] # TT-TDB + posvel_jpl_su(et) = kmsec2auday(spkgeo(10, et, "J2000", 0)[1]) # Sun + posvel_jpl_ea(et) = kmsec2auday(spkgeo(399, et, "J2000", 0)[1]) # Earth + posvel_jpl_mo(et) = kmsec2auday(spkgeo(301, et, "J2000", 0)[1]) # Moon + posvel_jpl_ma(et) = kmsec2auday(spkgeo(4, et, "J2000", 0)[1]) # Mars + posvel_jpl_ju(et) = kmsec2auday(spkgeo(5, et, "J2000", 0)[1]) # Jupiter + + tv = range(sol64.t0, sol64.t[end], 10) + for t in tv + et = t * daysec + @show t, et + @show abs(ttmtdb_jpl(et) - ttmtdb_pe(t)) + @show norm(posvel_jpl_su(et) - posvel_pe_su(t), Inf) + @show norm(posvel_jpl_ea(et) - posvel_pe_ea(t), Inf) + @show norm(posvel_jpl_mo(et) - posvel_pe_mo(t), Inf) + @show norm(posvel_jpl_ma(et) - posvel_pe_ma(t), Inf) + @show norm(posvel_jpl_ju(et) - posvel_pe_ju(t), Inf) + + @test abs(ttmtdb_jpl(et) - ttmtdb_pe(t)) < 1e-12 + @test norm(posvel_jpl_su(et) - posvel_pe_su(t), Inf) < 1e-14 + @test norm(posvel_jpl_ea(et) - posvel_pe_ea(t), Inf) < 1e-11 + @test norm(posvel_jpl_mo(et) - posvel_pe_mo(t), Inf) < 1e-11 + @test norm(posvel_jpl_ma(et) - posvel_pe_ma(t), Inf) < 1e-12 + @test norm(posvel_jpl_ju(et) - posvel_pe_ju(t), Inf) < 1e-13 + end + + # Remove files + rm.((filename, "naif0012.tls", "de430_1850-2150.bsp", "TTmTDB.de430.19feb2015.bsp")) + + # Float 128 + #= + # Test integration + sol128 = propagate(1, Float128(jd0), nyears, dense; dynamics = dynamics, order = order, abstol = abstol) + # Save solution + filename = selecteph2jld2(sol128, bodyind, nyears) + # Recovered solution + recovered_sol128 = JLD2.load(filename, "ss16ast_eph") + # Remove file + rm(filename) + + @test selecteph(sol128, bodyind, euler = true, ttmtdb = true) == recovered_sol128 + =# + end end From a747faf925e117e2d90831fbb85e4af858f06af4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jorge=20A=2E=20P=C3=A9rez-Hern=C3=A1ndez?= Date: Sun, 28 Jan 2024 22:28:56 +0100 Subject: [PATCH 06/10] Fixes in TaylorInterpolant --- src/interpolation.jl | 20 +++++++++++++------- test/propagation.jl | 1 + 2 files changed, 14 insertions(+), 7 deletions(-) diff --git a/src/interpolation.jl b/src/interpolation.jl index 7e3dc6e..b80deeb 100644 --- a/src/interpolation.jl +++ b/src/interpolation.jl @@ -25,7 +25,13 @@ time ``x(t)``; or a lunar core Euler angle as a function of time ``\theta_c(t)`` end end +const TaylorInterpCallingArgs{T,U} = Union{T, Taylor1{U}, TaylorN{U}, Taylor1{TaylorN{U}}} where {T,U} + # Outer constructors +function TaylorInterpolant{T, U, N}(t0::T, t::VT, x::X) where {T<:Number, U<:Number, N, VT<:AbstractVector{T}, X<:AbstractArray{Taylor1{U}, N}} + return TaylorInterpolant{T, U, N, VT, X}(t0, t, x) +end + function TaylorInterpolant(t0::T, t::VT, x::X) where {T<:Number, U<:Number, N, VT<:AbstractVector{T}, X<:AbstractArray{Taylor1{U}, N}} return TaylorInterpolant{T, U, N, VT, X}(t0, t, x) end @@ -65,7 +71,7 @@ end Return the index of `tinterp.t` corresponding to `t` and the time elapsed from `tinterp.t0` to `t`. """ -function getinterpindex(tinterp::TaylorInterpolant{T,U,N}, t::TT) where {T,U,N,TT<:Union{T, Taylor1{U}, Taylor1{TaylorN{U}}}} +function getinterpindex(tinterp::TaylorInterpolant{T,U,N}, t::TT) where {T,U,N,TT<:TaylorInterpCallingArgs{T,U}} t00 = constant_term(constant_term(t)) # Current time tmin, tmax = minmax(tinterp.t[end], tinterp.t[1]) # Min and max time in tinterp Δt = t-tinterp.t0 # Time since start of tinterp @@ -102,9 +108,9 @@ numberofbodies(interp::TaylorInterpolant) = numberofbodies(size(interp.x, 2)) @doc raw""" (tinterp::TaylorInterpolant{T, U, 1})(t::T) where {T, U} - (tinterp::TaylorInterpolant{T, U, 1})(t::TT) where {T, U, TT<:Union{Taylor1{U},Taylor1{TaylorN{U}}}} + (tinterp::TaylorInterpolant{T, U, 1})(t::TT) where {T, U, TT<:TaylorInterpCallingArgs{T,U}} (tinterp::TaylorInterpolant{T, U, 2})(t::T) where {T, U} - (tinterp::TaylorInterpolant{T, U, 2})(t::TT) where {T, U, TT<:Union{Taylor1{U},Taylor1{TaylorN{U}}}} + (tinterp::TaylorInterpolant{T, U, 2})(t::TT) where {T, U, TT<:TaylorInterpCallingArgs{T,U}} Evaluate `tinterp.x` at time `t`. @@ -116,11 +122,11 @@ function (tinterp::TaylorInterpolant{T, U, 1})(t::T) where {T, U} # Evaluate tinterp.x[ind] at δt return (tinterp.x[ind])(δt)::U end -function (tinterp::TaylorInterpolant{T, U, 1})(t::TT) where {T, U, TT<:Union{Taylor1{U},Taylor1{TaylorN{U}}}} +function (tinterp::TaylorInterpolant{T, U, 1})(t::TT) where {T, U, TT<:TaylorInterpCallingArgs{T,U}} # Get index of tinterp.x that interpolates at time t ind::Int, δt::TT = getinterpindex(tinterp, t) # Evaluate tinterp.x[ind] at δt - return (tinterp.x[ind])(δt)::Taylor1{TT} + return (tinterp.x[ind])(δt)::TT end function (tinterp::TaylorInterpolant{T, U, 2})(t::T) where {T, U} @@ -129,7 +135,7 @@ function (tinterp::TaylorInterpolant{T, U, 2})(t::T) where {T, U} # Evaluate tinterp.x[ind] at δt return view(tinterp.x, ind, :)(δt)::Vector{U} end -function (tinterp::TaylorInterpolant{T, U, 2})(t::TT) where {T, U, TT<:Union{Taylor1{U},Taylor1{TaylorN{U}}}} +function (tinterp::TaylorInterpolant{T, U, 2})(t::TT) where {T, U, TT<:TaylorInterpCallingArgs{T,U}} # Get index of tinterp.x that interpolates at time t ind::Int, δt::TT = getinterpindex(tinterp, t) # Evaluate tinterp.x[ind] at δt @@ -149,7 +155,7 @@ function (tinterp::TaylorInterpolant{T, U, 2})(target::Int, observer::Int, t::T) end end -(tinterp::TaylorInterpolant{T,U,2})(target::Int, t::TT) where {T, U, TT<:Union{T, Taylor1{U}, Taylor1{TaylorN{U}}}} = tinterp(target, 0, t) +(tinterp::TaylorInterpolant{T,U,2})(target::Int, t::TT) where {T, U, TT<:TaylorInterpCallingArgs{T,U}} = tinterp(target, 0, t) @doc raw""" reverse(tinterp::TaylorInterpolant{T,U,N}) where {T<:Real, U<:Number, N} diff --git a/test/propagation.jl b/test/propagation.jl index 28160ba..948c838 100644 --- a/test/propagation.jl +++ b/test/propagation.jl @@ -72,6 +72,7 @@ using LinearAlgebra: norm tmid = sol.t0 + sol.t[1]/2 @test sol(tmid) isa Vector{Float64} @test sol(tmid + Taylor1(order)) isa Vector{Taylor1{Float64}} + @test sol(tmid + dq[1] + dq[1]*dq[2]) isa Vector{TaylorN{Float64}} @test sol(tmid + Taylor1([dq[1],dq[1]*dq[2]],order)) isa Vector{Taylor1{TaylorN{Float64}}} end From 034e3239b30bac17712055049cb966a69d6aa4f5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jorge=20A=2E=20P=C3=A9rez-Hern=C3=A1ndez?= Date: Sun, 28 Jan 2024 22:59:47 +0100 Subject: [PATCH 07/10] Add U to TaylorInterpCallingArgs{T,U} --- src/interpolation.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/interpolation.jl b/src/interpolation.jl index b80deeb..44c45a1 100644 --- a/src/interpolation.jl +++ b/src/interpolation.jl @@ -25,7 +25,7 @@ time ``x(t)``; or a lunar core Euler angle as a function of time ``\theta_c(t)`` end end -const TaylorInterpCallingArgs{T,U} = Union{T, Taylor1{U}, TaylorN{U}, Taylor1{TaylorN{U}}} where {T,U} +const TaylorInterpCallingArgs{T,U} = Union{T, U, Taylor1{U}, TaylorN{U}, Taylor1{TaylorN{U}}} where {T,U} # Outer constructors function TaylorInterpolant{T, U, N}(t0::T, t::VT, x::X) where {T<:Number, U<:Number, N, VT<:AbstractVector{T}, X<:AbstractArray{Taylor1{U}, N}} From 38146a98aec12af572dd8762c2bb612fb1dc4250 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jorge=20A=2E=20P=C3=A9rez-Hern=C3=A1ndez?= Date: Mon, 29 Jan 2024 16:00:10 +0100 Subject: [PATCH 08/10] Try backwards-compatible fix for TaylorInterpolant serialization --- src/interpolation.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/interpolation.jl b/src/interpolation.jl index 44c45a1..48d8902 100644 --- a/src/interpolation.jl +++ b/src/interpolation.jl @@ -29,7 +29,7 @@ const TaylorInterpCallingArgs{T,U} = Union{T, U, Taylor1{U}, TaylorN{U}, Taylor1 # Outer constructors function TaylorInterpolant{T, U, N}(t0::T, t::VT, x::X) where {T<:Number, U<:Number, N, VT<:AbstractVector{T}, X<:AbstractArray{Taylor1{U}, N}} - return TaylorInterpolant{T, U, N, VT, X}(t0, t, x) + return TaylorInterpolant{T, U, N, Vector{T}, Array{Taylor1{U},N}}(t0, convert(Vector{T},t), convert(Array{Taylor1{U},N},x)) end function TaylorInterpolant(t0::T, t::VT, x::X) where {T<:Number, U<:Number, N, VT<:AbstractVector{T}, X<:AbstractArray{Taylor1{U}, N}} @@ -37,7 +37,7 @@ function TaylorInterpolant(t0::T, t::VT, x::X) where {T<:Number, U<:Number, N, V end function TaylorInterpolant(t0::T, t::SubArray{T, 1}, x::SubArray{Taylor1{U}, N}) where {T<:Number, U<:Number, N} - return TaylorInterpolant{T, U, N, Vector{T}, Array{Taylor1{U}, N}}(t0, t.parent[t.indices...], x.parent[x.indices...]) + return TaylorInterpolant{T, U, N}(t0, t, x) end # Custom print @@ -262,7 +262,7 @@ end writeas(::Type{TaylorInterpolant{T, T, 2, Vector{T}, Matrix{Taylor1{T}}}}) where {T<:Real} = PlanetaryEphemerisSerialization{T} # Convert method to write .jld2 files -function convert(::Type{PlanetaryEphemerisSerialization{T}}, eph::TaylorInterpolant{T, T, 2, Vector{T}, Matrix{Taylor1{T}}}) where {T <: Real} +function convert(::Type{PlanetaryEphemerisSerialization{T}}, eph::TaylorInterpolant{T, T, 2}) where {T <: Real} # Taylor polynomials order order = eph.x[1, 1].order # Number of coefficients in each polynomial @@ -282,7 +282,7 @@ function convert(::Type{PlanetaryEphemerisSerialization{T}}, eph::TaylorInterpol end # Convert method to read .jld2 files -function convert(::Type{TaylorInterpolant{T, T, 2, Vector{T}, Array{Taylor1{T},2}}}, eph::PlanetaryEphemerisSerialization{T}) where {T<:Real} +function convert(::Type{TaylorInterpolant{T, T, 2}}, eph::PlanetaryEphemerisSerialization{T}) where {T<:Real} # Taylor polynomials order order = eph.order # Number of coefficients in each polynomial @@ -298,7 +298,7 @@ function convert(::Type{TaylorInterpolant{T, T, 2, Vector{T}, Array{Taylor1{T},2 x[i] = Taylor1{T}(eph.x[(i-1)*k+1 : i*k], order) end - return TaylorInterpolant{T, T, 2, Vector{T}, Matrix{Taylor1{T}}}(eph.t0, eph.t, x) + return TaylorInterpolant{T, T, 2}(eph.t0, eph.t, x) end @doc raw""" From fe550b28ebee5ba0d1a633eeea0f2a82dd7cbd90 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jorge=20A=2E=20P=C3=A9rez-Hern=C3=A1ndez?= Date: Mon, 29 Jan 2024 16:39:48 +0100 Subject: [PATCH 09/10] Fix writeas? --- src/interpolation.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/interpolation.jl b/src/interpolation.jl index 48d8902..eca7587 100644 --- a/src/interpolation.jl +++ b/src/interpolation.jl @@ -259,7 +259,7 @@ struct PlanetaryEphemerisSerialization{T} end # Tell JLD2 to save TaylorInterpolant{T, T, 2} as PlanetaryEphemerisSerialization{T} -writeas(::Type{TaylorInterpolant{T, T, 2, Vector{T}, Matrix{Taylor1{T}}}}) where {T<:Real} = PlanetaryEphemerisSerialization{T} +writeas(::Type{TaylorInterpolant{T, T, 2}}) where {T<:Real} = PlanetaryEphemerisSerialization{T} # Convert method to write .jld2 files function convert(::Type{PlanetaryEphemerisSerialization{T}}, eph::TaylorInterpolant{T, T, 2}) where {T <: Real} From 3b09239cddca73b6dba1b1958763876ee990f57b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jorge=20A=2E=20P=C3=A9rez-Hern=C3=A1ndez?= Date: Mon, 29 Jan 2024 20:44:46 +0100 Subject: [PATCH 10/10] Bump patch version --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index ed3c42e..78e4928 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "PlanetaryEphemeris" uuid = "d83715d0-7e5f-11e9-1a59-4137b20d8363" authors = ["Jorge A. Pérez Hernández", "Luis Benet", "Luis Eduardo Ramírez Montoya"] -version = "0.7.5" +version = "0.7.6" [deps] ArgParse = "c7e460c6-2fb9-53a9-8c5b-16f535851c63"