Skip to content

Commit

Permalink
Update to TaylorIntegration v0.14; add loadeph (#23)
Browse files Browse the repository at this point in the history
  • Loading branch information
@PerezHz
PerezHz authored Jun 30, 2023
1 parent 79c44b6 commit bb49d32
Show file tree
Hide file tree
Showing 5 changed files with 69 additions and 26 deletions.
4 changes: 2 additions & 2 deletions Project.toml
View file
Original file line number Diff line number Diff line change
@@ -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"
Expand All @@ -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"
36 changes: 15 additions & 21 deletions src/plephinteg.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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}
Expand All @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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)

Expand All @@ -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
Expand All @@ -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]
Expand All @@ -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
Expand All @@ -256,17 +253,14 @@ 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
@inbounds t[0] = t0
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
Expand All @@ -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)
Expand All @@ -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
Expand Down
52 changes: 52 additions & 0 deletions src/propagation.jl
View file
Original file line number Diff line number Diff line change
@@ -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}
Expand Down
2 changes: 0 additions & 2 deletions test/Project.toml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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"
1 change: 0 additions & 1 deletion test/propagation.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,6 @@
using PlanetaryEphemeris
using Dates
using Quadmath
using TaylorIntegration
using JLD2
using Test

Expand Down

2 comments on commit bb49d32

@PerezHz
Copy link
Owner Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@JuliaRegistrator register

@JuliaRegistrator
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Registration pull request created: JuliaRegistries/General/86554

After the above pull request is merged, it is recommended that a tag is created on this repository for the registered package version.

This will be done automatically if the Julia TagBot GitHub Action is installed, or can be done manually through the github interface, or via:

git tag -a v0.7.2 -m "<description of version>" bb49d3212e1e14cdf49a6d65798743488e5eb3d7
git push origin v0.7.2

Please sign in to comment.