From bb49d3212e1e14cdf49a6d65798743488e5eb3d7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jorge=20P=C3=A9rez?= Date: Fri, 30 Jun 2023 05:23:41 +0200 Subject: [PATCH] Update to TaylorIntegration v0.14; add `loadeph` (#23) --- Project.toml | 4 ++-- src/plephinteg.jl | 36 +++++++++++++------------------ src/propagation.jl | 52 +++++++++++++++++++++++++++++++++++++++++++++ test/Project.toml | 2 -- test/propagation.jl | 1 - 5 files changed, 69 insertions(+), 26 deletions(-) diff --git a/Project.toml b/Project.toml index 0df0f22..0e4f421 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.1" +version = "0.7.2" [deps] ArgParse = "c7e460c6-2fb9-53a9-8c5b-16f535851c63" @@ -23,6 +23,6 @@ DelimitedFiles = "1" JLD2 = "0.4" PrecompileTools = "1.1" Quadmath = "0.5" -TaylorIntegration = "0.13" +TaylorIntegration = "0.14" TaylorSeries = "0.15" julia = "1.6" diff --git a/src/plephinteg.jl b/src/plephinteg.jl index 54e702d..c828569 100644 --- a/src/plephinteg.jl +++ b/src/plephinteg.jl @@ -6,11 +6,11 @@ Threaded version of `TaylorSeries.evaluate!`. See also [`TaylorSeries.evaluate!`](@ref). """ function evaluate_threads!(x::Vector{Taylor1{T}}, δt::T, x0::Vector{T}) where { T <: Number} - + Threads.@threads for i in eachindex(x) x0[i] = evaluate( x[i], δt ) end - + nothing end @@ -46,8 +46,8 @@ end First step-size control. See section 3.2 of https://doi.org/10.1080/10586458.2005.10128904. -See also [`stepsize_threads`](@ref) and [`TaylorIntegration.stepsize`](@ref). -""" +See also [`stepsize_threads`](@ref) and [`TaylorIntegration.stepsize`](@ref). +""" function stepsize_jz05(q::AbstractArray{Taylor1{U}, N}, epsilon::T) where {T<:Real, U<:Number, N} nbodies = (length(q)-13)÷6 @@ -88,7 +88,7 @@ end # end @doc raw""" - taylorstep_threads!(f!, t::Taylor1{T}, x::Vector{Taylor1{U}}, dx::Vector{Taylor1{U}}, xaux::Vector{Taylor1{U}}, + taylorstep_threads!(f!, t::Taylor1{T}, x::Vector{Taylor1{U}}, dx::Vector{Taylor1{U}}, xaux::Vector{Taylor1{U}}, abstol::T, params, parse_eqs::Bool=true) where {T<:Real, U<:Number} taylorstep_threads!(f!, t::Taylor1{T}, x::Vector{Taylor1{U}}, dx::Vector{Taylor1{U}}, abstol::T, params, rv::TaylorIntegration.RetAlloc{Taylor1{U}}) where {T<:Real, U<:Number} @@ -97,7 +97,7 @@ Threaded version of `TaylorIntegration.taylorstep`. See also [`stepsize_threads`](@ref) and [`TaylorIntegration.taylorstep`](@ref). """ -function taylorstep_threads!(f!, t::Taylor1{T}, x::Vector{Taylor1{U}}, dx::Vector{Taylor1{U}}, xaux::Vector{Taylor1{U}}, +function taylorstep_threads!(f!, t::Taylor1{T}, x::Vector{Taylor1{U}}, dx::Vector{Taylor1{U}}, xaux::Vector{Taylor1{U}}, abstol::T, params) where {T<:Real, U<:Number} # Compute the Taylor coefficients @@ -183,10 +183,10 @@ for V in (:(Val{true}), :(Val{false})) else return _taylorinteg_threads!(f!, t, x, dx, q0, t0, tmax, abstol, $V(), params, maxsteps = maxsteps) end - + end - function _taylorinteg_threads!(f!, t::Taylor1{T}, x::Array{Taylor1{U}, 1}, dx::Array{Taylor1{U}, 1}, q0::Array{U, 1}, t0::T, + function _taylorinteg_threads!(f!, t::Taylor1{T}, x::Array{Taylor1{U}, 1}, dx::Array{Taylor1{U}, 1}, q0::Array{U, 1}, t0::T, tmax::T, abstol::T, ::$V, params; maxsteps::Int = 500) where {T <: Real, U <: Number} # Initialize the vector of Taylor1 expansions @@ -196,7 +196,7 @@ for V in (:(Val{true}), :(Val{false})) tv = Array{T}(undef, maxsteps+1) xv = Array{U}(undef, dof, maxsteps+1) if $V == Val{true} - polynV = Array{Taylor1{U}}(undef, dof, maxsteps+1) + psol = Array{Taylor1{U}}(undef, dof, maxsteps) end xaux = Array{Taylor1{U}}(undef, dof) @@ -206,9 +206,6 @@ for V in (:(Val{true}), :(Val{false})) x0 = deepcopy(q0) @inbounds tv[1] = t0 @inbounds xv[:,1] .= q0 - if $V == Val{true} - @inbounds polynV[:,1] .= deepcopy.(x) - end sign_tstep = copysign(1, tmax-t0) # Integration @@ -220,7 +217,7 @@ for V in (:(Val{true}), :(Val{false})) evaluate_threads!(x, δt, x0) # new initial condition if $V == Val{true} # Store the Taylor polynomial solution - @inbounds polynV[:,nsteps+1] .= deepcopy.(x) + @inbounds psol[:,nsteps] .= deepcopy.(x) end @inbounds Threads.@threads for i in eachindex(x0) x[i][0] = x0[i] @@ -240,13 +237,13 @@ for V in (:(Val{true}), :(Val{false})) end if $V == Val{true} - return TaylorInterpolant(tv[1], view(tv.-tv[1],1:nsteps), view(transpose(view(polynV,:,2:nsteps)),1:nsteps-1,:)) + return TaylorInterpolant(tv[1], view(tv.-tv[1],1:nsteps), view(transpose(view(psol,:,1:nsteps-1)),1:nsteps-1,:)) elseif $V == Val{false} return view(tv,1:nsteps), view(transpose(view(xv,:,1:nsteps)),1:nsteps,:) end end - function _taylorinteg_threads!(f!, t::Taylor1{T}, x::Array{Taylor1{U}, 1}, dx::Array{Taylor1{U}, 1}, q0::Array{U, 1}, t0::T, + function _taylorinteg_threads!(f!, t::Taylor1{T}, x::Array{Taylor1{U}, 1}, dx::Array{Taylor1{U}, 1}, q0::Array{U, 1}, t0::T, tmax::T, abstol::T, rv::TaylorIntegration.RetAlloc{Taylor1{U}}, ::$V, params; maxsteps::Int = 500) where {T <: Real, U <: Number} # Initialize the vector of Taylor1 expansions @@ -256,7 +253,7 @@ for V in (:(Val{true}), :(Val{false})) tv = Array{T}(undef, maxsteps+1) xv = Array{U}(undef, dof, maxsteps+1) if $V == Val{true} - polynV = Array{Taylor1{U}}(undef, dof, maxsteps+1) + psol = Array{Taylor1{U}}(undef, dof, maxsteps) end # Initial conditions @@ -264,9 +261,6 @@ for V in (:(Val{true}), :(Val{false})) x0 = deepcopy(q0) @inbounds tv[1] = t0 @inbounds xv[:,1] .= q0 - if $V == Val{true} - @inbounds polynV[:,1] .= deepcopy.(x) - end sign_tstep = copysign(1, tmax-t0) # Integration @@ -278,7 +272,7 @@ for V in (:(Val{true}), :(Val{false})) evaluate_threads!(x, δt, x0) # new initial condition if $V == Val{true} # Store the Taylor polynomial solution - @inbounds polynV[:,nsteps+1] .= deepcopy.(x) + @inbounds psol[:,nsteps] .= deepcopy.(x) end Threads.@threads for i in eachindex(x0) @@ -299,7 +293,7 @@ for V in (:(Val{true}), :(Val{false})) end if $V == Val{true} - return TaylorInterpolant(tv[1], view(tv.-tv[1],1:nsteps), view(transpose(view(polynV,:,2:nsteps)),1:nsteps-1,:)) + return TaylorInterpolant(tv[1], view(tv.-tv[1],1:nsteps), view(transpose(view(psol,:,1:nsteps-1)),1:nsteps-1,:)) elseif $V == Val{false} return view(tv,1:nsteps), view(transpose(view(xv,:,1:nsteps)),1:nsteps,:) end diff --git a/src/propagation.jl b/src/propagation.jl index 6694d0f..652ab82 100644 --- a/src/propagation.jl +++ b/src/propagation.jl @@ -1,3 +1,55 @@ +@doc raw""" + loadeph(ss16asteph::TaylorInterpolant, μ::Vector{<:Real}) + +Taking Solar System ephemeris `ss16asteph` and their gravitational parameters `μ` as input, returns for all bodies the point-mass Newtonian acceleration and the Newtonian N body potential. + +# Arguments + +- `ss16asteph`: Solar System ephemeris. +- `μ::Vector{<:Real}`: vector of mass parameters. +""" +function loadeph(ss16asteph::TaylorInterpolant, μ::Vector{<:Real}) + + # Compute point-mass Newtonian accelerations from ephemeris + # accelerations of all bodies are needed to compute the post-Newtonian acceleration of e.g. Solar System minor bodies + # Number of bodies that contibute to the asteroid's acceleration + Nm1 = numberofbodies(ss16asteph) + # Initialize a TaylorInterpolant for the point-mass Newtonian accelerations + acc_eph = TaylorInterpolant(ss16asteph.t0, ss16asteph.t, Matrix{eltype(ss16asteph.x)}(undef, length(ss16asteph.t)-1, 3Nm1)) + # Initialize a TaylorInterpolant for the newtonian N body potential + pot_eph = TaylorInterpolant(ss16asteph.t0, ss16asteph.t, Matrix{eltype(ss16asteph.x)}(undef, length(ss16asteph.t)-1, Nm1)) + # Fill TaylorInterpolant.x with zero polynomials + fill!(acc_eph.x, zero(ss16asteph.x[1])) + fill!(pot_eph.x, zero(ss16asteph.x[1])) + + # Iterator over all bodies except asteroid + for j in 1:Nm1 + for i in 1:Nm1 + if i == j + # + else + # Difference between two positions (\mathbf{r}_i - \mathbf{r}_j) + X_ij = ss16asteph.x[:, 3i-2] .- ss16asteph.x[:, 3j-2] # X-axis component + Y_ij = ss16asteph.x[:, 3i-1] .- ss16asteph.x[:, 3j-1] # Y-axis component + Z_ij = ss16asteph.x[:, 3i ] .- ss16asteph.x[:, 3j ] # Z-axis component + # Distance between two bodies squared ||\mathbf{r}_i - \mathbf{r}_j||^2 + r_p2_ij = ( (X_ij.^2) .+ (Y_ij.^2) ) .+ (Z_ij.^2) + # Distance between two bodies ||\mathbf{r}_i - \mathbf{r}_j|| + r_ij = sqrt.(r_p2_ij) + # Newtonian potential + pot_eph.x[:, j] .+= (μ[i]./r_ij) + end + end + + # Fill acelerations by differentiating velocities + acc_eph.x[:, 3j-2] .= NEOs.ordpres_differentiate.(ss16asteph.x[:, 3(Nm1+j)-2]) # X-axis component + acc_eph.x[:, 3j-1] .= NEOs.ordpres_differentiate.(ss16asteph.x[:, 3(Nm1+j)-1]) # Y-axis component + acc_eph.x[:, 3j ] .= NEOs.ordpres_differentiate.(ss16asteph.x[:, 3(Nm1+j) ]) # Z-axis component + end + + return acc_eph, pot_eph +end + @doc raw""" selecteph2jld2(sseph::TaylorInterpolant, bodyind::T, tspan::S) where {T <: AbstractVector{Int}, S <: Number} diff --git a/test/Project.toml b/test/Project.toml index 94452d1..403e0c3 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -5,11 +5,9 @@ JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Quadmath = "be4d8f0f-7fa4-5f49-b795-2f01399ab2dd" SPICE = "5bab7191-041a-5c2e-a744-024b9c3a5062" -TaylorIntegration = "92b13dbe-c966-51a2-8445-caca9f8a7d42" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [compat] JLD2 = "0.4" Quadmath = "0.5" SPICE = "0.2" -TaylorIntegration = "0.13" diff --git a/test/propagation.jl b/test/propagation.jl index dfe6f85..43ca1a6 100644 --- a/test/propagation.jl +++ b/test/propagation.jl @@ -3,7 +3,6 @@ using PlanetaryEphemeris using Dates using Quadmath -using TaylorIntegration using JLD2 using Test