diff --git a/src/interpolation/PlanetaryEphemerisSerialization.jl b/src/interpolation/PlanetaryEphemerisSerialization.jl index 67dd4bc..cc5df78 100644 --- a/src/interpolation/PlanetaryEphemerisSerialization.jl +++ b/src/interpolation/PlanetaryEphemerisSerialization.jl @@ -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 @@ -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 @@ -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 \ No newline at end of file diff --git a/src/interpolation/TaylorInterpolantNSerialization.jl b/src/interpolation/TaylorInterpolantNSerialization.jl index 88c5924..33c6c2d 100644 --- a/src/interpolation/TaylorInterpolantNSerialization.jl +++ b/src/interpolation/TaylorInterpolantNSerialization.jl @@ -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 @@ -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) @@ -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 @@ -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 \ No newline at end of file diff --git a/src/propagation.jl b/src/propagation.jl index 73ef53c..87a5a16 100644 --- a/src/propagation.jl +++ b/src/propagation.jl @@ -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 diff --git a/test/propagation.jl b/test/propagation.jl index a964612..6da51aa 100644 --- a/test/propagation.jl +++ b/test/propagation.jl @@ -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