Skip to content

Commit

Permalink
Scattered fixes
Browse files Browse the repository at this point in the history
  • Loading branch information
@PerezHz
PerezHz committed Apr 1, 2024
1 parent 39aba75 commit 136be6e
Show file tree
Hide file tree
Showing 4 changed files with 115 additions and 104 deletions.
12 changes: 9 additions & 3 deletions src/interpolation/PlanetaryEphemerisSerialization.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,10 +21,12 @@ struct PlanetaryEphemerisSerialization{T}
end

# Tell JLD2 to save TaylorInterpolant{T, T, 2} as PlanetaryEphemerisSerialization{T}
writeas(::Type{TaylorInterpolant{T, T, 2}}) where {T<:Real} = PlanetaryEphemerisSerialization{T}
function writeas(::Type{<:TaylorInterpolant{T, T, 2, Vector{T}, Matrix{Taylor1{T}}}}) where {T<:Real}
return PlanetaryEphemerisSerialization{T}
end

# Convert method to write .jld2 files
function convert(::Type{PlanetaryEphemerisSerialization{T}}, eph::TaylorInterpolant{T, T, 2}) where {T <: Real}
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
Expand All @@ -44,7 +46,7 @@ function convert(::Type{PlanetaryEphemerisSerialization{T}}, eph::TaylorInterpol
end

# Convert method to read .jld2 files
function convert(::Type{TaylorInterpolant{T, T, 2}}, eph::PlanetaryEphemerisSerialization{T}) where {T<:Real}
function convert(::Type{TaylorInterpolant{T, T, 2, Vector{T}, Matrix{Taylor1{T}}}}, eph::PlanetaryEphemerisSerialization{T}) where {T<:Real}
# Taylor polynomials order
order = eph.order
# Number of coefficients in each polynomial
Expand All @@ -62,3 +64,7 @@ function convert(::Type{TaylorInterpolant{T, T, 2}}, eph::PlanetaryEphemerisSeri

return TaylorInterpolant{T, T, 2}(eph.t0, eph.t, x)
end

function convert(::Type{TaylorInterpolant{T, T, 2}}, eph::PlanetaryEphemerisSerialization{T}) where {T<:Real}
return convert(TaylorInterpolant{T, T, 2, Vector{T}, Matrix{Taylor1{T}}}, eph)
end
15 changes: 10 additions & 5 deletions src/interpolation/TaylorInterpolantNSerialization.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,11 +22,13 @@ struct TaylorInterpolantNSerialization{T}
x::Vector{T}
end

# Tell JLD2 to save TaylorInterpolant{T, TaylorN{T}, 2} as TaylorInterpolantNSerialization{T}
writeas(::Type{TaylorInterpolant{T, TaylorN{T}, 2}}) where {T} = TaylorInterpolantNSerialization{T}
# Tell JLD2 to save <:TaylorInterpolant{T, TaylorN{T}, 2} as TaylorInterpolantNSerialization{T}
function writeas(::Type{<:TaylorInterpolant{T, TaylorN{T}, 2, Vector{T}, Matrix{Taylor1{TaylorN{T}}}}}) where {T<:Real}
return TaylorInterpolantNSerialization{T}
end

# Convert method to write .jld2 files
function convert(::Type{TaylorInterpolantNSerialization{T}}, eph::TaylorInterpolant{T, TaylorN{T}, 2}) where {T}
function convert(::Type{TaylorInterpolantNSerialization{T}}, eph::TaylorInterpolant{T, TaylorN{T}, 2, Vector{T}, Matrix{Taylor1{TaylorN{T}}}}) where {T}
# Variables
vars = TS.get_variable_names()
# Number of variables
Expand All @@ -45,7 +47,6 @@ function convert(::Type{TaylorInterpolantNSerialization{T}}, eph::TaylorInterpol
L = varorder + 1
# Number of coefficients in each HomogeneousPolynomial
M = binomial(n + varorder, varorder)
# M = sum(binomial(n + i_3 - 1, i_3) for i_3 in 0:varorder)

# Vector of coefficients
x = Vector{T}(undef, N * k * M)
Expand All @@ -71,7 +72,7 @@ function convert(::Type{TaylorInterpolantNSerialization{T}}, eph::TaylorInterpol
end

# Convert method to read .jld2 files
function convert(::Type{TaylorInterpolant{T, TaylorN{T}, 2}}, eph::TaylorInterpolantNSerialization{T}) where {T}
function convert(::Type{TaylorInterpolant{T, TaylorN{T}, 2, Vector{T}, Matrix{Taylor1{TaylorN{T}}}}}, eph::TaylorInterpolantNSerialization{T}) where {T}
# Variables
vars = eph.vars
# Number of variables
Expand Down Expand Up @@ -120,4 +121,8 @@ function convert(::Type{TaylorInterpolant{T, TaylorN{T}, 2}}, eph::TaylorInterpo
end

return TaylorInterpolant{T, TaylorN{T}, 2}(eph.t0, eph.t, x)
end

function convert(::Type{TaylorInterpolant{T, TaylorN{T}, 2}}, eph::PlanetaryEphemerisSerialization{T}) where {T<:Real}
return convert(TaylorInterpolant{T, TaylorN{T}, 2, Vector{T}, Matrix{Taylor1{TaylorN{T}}}}, eph)
end
2 changes: 1 addition & 1 deletion src/propagation.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -195,7 +195,7 @@ for V_dense in (:(Val{true}), :(Val{false}))
maxsteps, parse_eqs)

if $V_dense == Val{true}
return TaylorInterpolant{T, T, 2, SubArray{T,1}, SubArray{Taylor1{T},2}}(jd0 - J2000, sol[1], sol[3])
return TaylorInterpolant{T, T, 2}(jd0 - J2000, sol[1], sol[3])
else
return sol
end
Expand Down
190 changes: 95 additions & 95 deletions test/propagation.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -97,100 +97,100 @@ using LinearAlgebra: norm
rm("test.jld2")
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

# Test selecteph
t0 = sol64.t0 + sol64.t[end]/3
tf = sol64.t0 + 2*sol64.t[end]/3
idxs = vcat(nbodyind(N, [su, ea, mo]), 6N+1:6N+13)
i_0 = searchsortedlast(sol64.t, t0)
i_f = searchsortedfirst(sol64.t, tf)

subsol = selecteph(sol64, [su, ea, mo], t0, tf; euler = true, ttmtdb = true)

@test subsol.t0 == sol64.t0
@test subsol.t0 + subsol.t[1] ≤ t0
@test subsol.t0 + subsol.t[end] ≥ tf
@test size(subsol.x) == (i_f - i_0, length(idxs))
@test subsol.x == sol64.x[i_0:i_f-1, idxs]
@test subsol(t0) == sol64(t0)[idxs]
@test subsol(tf) == sol64(tf)[idxs]

# Kernels URLs
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

# # Test selecteph
# t0 = sol64.t0 + sol64.t[end]/3
# tf = sol64.t0 + 2*sol64.t[end]/3
# idxs = vcat(nbodyind(N, [su, ea, mo]), 6N+1:6N+13)
# i_0 = searchsortedlast(sol64.t, t0)
# i_f = searchsortedfirst(sol64.t, tf)

# subsol = selecteph(sol64, [su, ea, mo], t0, tf; euler = true, ttmtdb = true)

# @test subsol.t0 == sol64.t0
# @test subsol.t0 + subsol.t[1] ≤ t0
# @test subsol.t0 + subsol.t[end] ≥ tf
# @test size(subsol.x) == (i_f - i_0, length(idxs))
# @test subsol.x == sol64.x[i_0:i_f-1, idxs]
# @test subsol(t0) == sol64(t0)[idxs]
# @test subsol(tf) == sol64(tf)[idxs]

# # Kernels URLs
# 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

0 comments on commit 136be6e

Please sign in to comment.