Skip to content

Commit

Permalink
Uncomment tests
Browse files Browse the repository at this point in the history
  • Loading branch information
@PerezHz
PerezHz committed Apr 1, 2024
1 parent 136be6e commit fe9240b
Showing 1 changed file with 95 additions and 95 deletions.
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 fe9240b

Please sign in to comment.