TaylorIntegration v0.16 update patch

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.8.5"
+version = "0.8.6"
 
 [deps]
 ArgParse = "c7e460c6-2fb9-53a9-8c5b-16f535851c63"
@@ -23,6 +23,6 @@ DelimitedFiles = "1"
 JLD2 = "0.4"
 PrecompileTools = "1.1"
 Quadmath = "0.5"
-TaylorIntegration = "0.15.3"
+TaylorIntegration = "0.16"
 TaylorSeries = "0.18"
 julia = "1.6"

scripts/integrate_ephemeris.jl

@@ -3,17 +3,17 @@ using ArgParse, PlanetaryEphemeris, Dates
 function parse_commandline()
 
     s = ArgParseSettings()
-    
+
     # Program name (for usage & help screen)
-    s.prog = "integrate_ephemeris.jl"  
+    s.prog = "integrate_ephemeris.jl"
     # Desciption (for help screen)
-    s.description = "Integrates JPL DE430 Ephemeris" 
+    s.description = "Integrates JPL DE430 Ephemeris"
 
     @add_arg_table! s begin
         "--maxsteps"
             help = "Maximum number of steps during integration"
             arg_type = Int
-            default = 1_000_000 
+            default = 1_000_000
         "--jd0"
             help = "Starting time of integration; options are: \"2000-01-01T12:00:00\" or \"1969-06-28T00:00:00\""
             arg_type = DateTime
@@ -41,7 +41,7 @@ function parse_commandline()
         "--parse_eqs"
             help = "Whether to use the taylorized method of jetcoeffs (a-priori faster) or not"
             arg_type = Bool
-            default = true 
+            default = true
         "--bodyind"
             help = "Body indices in output"
             arg_type = UnitRange{Int}
@@ -68,43 +68,43 @@ function print_header(header::String)
     println(repeat("-", L))
     println(header)
     println(repeat("-", L))
-end 
+end
 
 function main(maxsteps::Int, jd0_datetime::DateTime, nyears::Float64, dynamics::Function, nast::Int,
               bodyind::UnitRange{Int}, order::Int, abstol::Float64, parse_eqs::Bool)
-    
-    jd0 = datetime2julian(jd0_datetime) 
+
+    jd0 = datetime2julian(jd0_datetime)
     print_header("Integrator warmup")
-    _ = propagate(1, jd0, nyears, Val(true), dynamics = dynamics, nast = nast, order = order, abstol = abstol, 
+    _ = propagate(1, jd0, nyears, dynamics = dynamics, nast = nast, order = order, abstol = abstol,
                   parse_eqs = parse_eqs)
-    
+
     print_header("Full integration")
     println( "Initial time of integration: ", string(jd0_datetime) )
     println( "Final time of integration: ", string(julian2datetime(jd0 + nyears*yr)) )
-    sol = propagate(maxsteps, jd0, nyears, Val(true), dynamics = dynamics, nast = nast, order = order, abstol = abstol, 
+    sol = propagate(maxsteps, jd0, nyears, dynamics = dynamics, nast = nast, order = order, abstol = abstol,
                     parse_eqs = parse_eqs)
 
     selecteph2jld2(sol, bodyind, nyears)
 
-    nothing 
-end 
+    nothing
+end
 
 function main()
 
-    # Parse arguments from commandline 
+    # Parse arguments from commandline
     parsed_args = parse_commandline()
 
     print_header("Integrate Ephemeris")
     print_header("General parameters")
 
-    # Number of threads 
+    # Number of threads
     println("• Number of threads: ", Threads.nthreads())
 
-    # Dynamical function 
+    # Dynamical function
     dynamics = parsed_args["dynamics"]
     println("• Dynamical function: ", dynamics)
 
-    # Maximum number of steps 
+    # Maximum number of steps
     maxsteps = parsed_args["maxsteps"]
     println("• Maximum number of steps: ", maxsteps)
 
@@ -116,21 +116,21 @@ function main()
     abstol = parsed_args["abstol"]
     println("• Absolute tolerance: ", abstol)
 
-    # Wheter to use @taylorize or not 
+    # Wheter to use @taylorize or not
     parse_eqs = parsed_args["parse_eqs"]
     println("• Use @taylorize: ", parse_eqs)
 
-    # Initial date o fintegration 
+    # Initial date o fintegration
     jd0_datetime = parsed_args["jd0"]
-    # Number of years to be integrated 
+    # Number of years to be integrated
     nyears = parsed_args["nyears"]
-    # Number of asteroids to be saved 
+    # Number of asteroids to be saved
     nast = parsed_args["nast"]
-    # Indexes of bodies to be saved 
+    # Indexes of bodies to be saved
     bodyind = parsed_args["bodyind"]
 
     main(maxsteps, jd0_datetime, nyears, dynamics, nast, bodyind, order, abstol, parse_eqs)
-    
+
 end
 
 main()

src/PlanetaryEphemeris.jl

@@ -33,6 +33,7 @@ include("propagation.jl")
 include("osculating.jl")
 include("barycenter.jl")
 #include("precompile.jl")
+include("deprecated.jl")
 
 include("dynamics/trivial.jl")
 include("dynamics/dynamical_model.jl")

src/deprecated.jl

@@ -0,0 +1,5 @@
+Base.@deprecate propagate(maxsteps, jd0, tspan, ::Val{true}; dynamics = NBP_pN_A_J23E_J23M_J2S_threads!,
+                   nast = 343, order = order, abstol = abstol, parse_eqs = true) propagate(maxsteps, jd0, tspan; dynamics = NBP_pN_A_J23E_J23M_J2S_threads!,
+                   nast = 343, order = order, abstol = abstol, parse_eqs = true)
+Base.@deprecate propagate(maxsteps, jd0, tspan, ::Val{false}; dynamics = NBP_pN_A_J23E_J23M_J2S_threads!,
+                   nast = 343, order = order, abstol = abstol, parse_eqs = true) error("PlanetaryEphemeris.propagate may only be called with `Val(true)` in the fourth argument.")

src/precompile.jl

@@ -7,13 +7,11 @@ using PrecompileTools
     const jd0 = datetime2julian(DateTime(2000,1,1,12))
     # Number of years
     const nyears = 2031.0 - year(julian2datetime(jd0))
-    # Dense output
-    const dense = Val(true)
     # Dynamical function
     const dynamics = DE430!
     @compile_workload begin
         # All calls in this block will be precompiled, regardless of whether
         # they belong to your package or not (on Julia 1.8 and higher)
-        propagate(maxsteps, jd0, nyears, dense; dynamics = dynamics, order = order, abstol = abstol)
+        propagate(maxsteps, jd0, nyears; dynamics = dynamics, order = order, abstol = abstol)
     end
 end

src/propagation.jl

@@ -151,7 +151,7 @@ function save2jld2andcheck(outfilename::String, sol)
 end
 
 @doc raw"""
-    propagate(maxsteps::Int, jd0::T, tspan::T, ::Val{false/true}; dynamics::Function = NBP_pN_A_J23E_J23M_J2S_threads!,
+    propagate(maxsteps::Int, jd0::T, tspan::T; dynamics::Function = NBP_pN_A_J23E_J23M_J2S_threads!,
               nast::Int = 343, order::Int = order, abstol::T = abstol, parse_eqs::Bool = true) where {T <: Real}
 
 Integrate the Solar System via the Taylor method.
@@ -161,56 +161,44 @@ Integrate the Solar System via the Taylor method.
 - `maxsteps::Int`: maximum number of steps for the integration.
 - `jd0::T`: initial Julian date.
 - `tspan::T`: time span of the integration (in Julian days).
-- `::Val{false/true}`: whether to save the Taylor polynomials at each step (`true`) or not (`false`).
 - `dynamics::Function`: dynamical model function.
 - `nast::Int`: number of asteroids to be considered in the integration.
 - `order::Int`: order of the Taylor expansions to be used in the integration.
 - `abstol::T`: absolute tolerance.
 - `parse_eqs::Bool`: whether to use the specialized method of `jetcoeffs!` (`true`) created with `@taylorize` or not.
 """ propagate
 
-for V_dense in (:(Val{true}), :(Val{false}))
-    @eval begin
+function propagate(maxsteps::Int, jd0::T, tspan::T; dynamics::Function = NBP_pN_A_J23E_J23M_J2S_threads!,
+                    nast::Int = 343, order::Int = order, abstol::T = abstol, parse_eqs::Bool = true) where {T <: Real}
 
-        function propagate(maxsteps::Int, jd0::T, tspan::T, ::$V_dense; dynamics::Function = NBP_pN_A_J23E_J23M_J2S_threads!,
-                           nast::Int = 343, order::Int = order, abstol::T = abstol, parse_eqs::Bool = true) where {T <: Real}
+    # Total number of bodies (Sun + 8 planets + Moon + Pluto + Asteroid)
+    N = 11 + nast
 
-            # Total number of bodies (Sun + 8 planets + Moon + Pluto + Asteroid)
-            N = 11 + nast
+    # Get 6N + 13 initial conditions (3N positions + 3N velocities + 6 lunar mantle angles + 6 lunar core angles + TT-TDB)
+    q0 = initialcond(N, jd0)
 
-            # Get 6N + 13 initial conditions (3N positions + 3N velocities + 6 lunar mantle angles + 6 lunar core angles + TT-TDB)
-            q0 = initialcond(N, jd0)
+    # Set initial time equal to zero (improves accuracy in data reductions)
+    t0 = zero(T)
 
-            # Set initial time equal to zero (improves accuracy in data reductions)
-            t0 = zero(T)
+    # Parameters for dynamical function
+    params = (N, jd0)
 
-            # Parameters for dynamical function
-            params = (N, jd0)
+    # Final time of integration (days)
+    tmax = t0 + tspan*yr
 
-            # Final time of integration (days)
-            tmax = t0 + tspan*yr
+    # Integration
+    sol = @time taylorinteg(dynamics, q0, t0, tmax, order, abstol, params;
+                            maxsteps, parse_eqs)
 
-            # Integration
-            sol = @time taylorinteg(dynamics, q0, t0, tmax, order, abstol, $V_dense(), params;
-                                    maxsteps, parse_eqs)
-
-            if $V_dense == Val{true}
-                return TaylorInterpolant{T, T, 2}(jd0 - J2000, sol[1], sol[3])
-            else
-                return sol
-            end
-
-        end
-
-        function propagate(maxsteps::Int, jd0::T1, tspan::T2, ::$V_dense; dynamics::Function = NBP_pN_A_J23E_J23M_J2S_threads!,
-                           nast::Int = 343, order::Int = order, abstol::T3 = abstol, parse_eqs::Bool = true) where {T1, T2, T3 <: Real}
+    return TaylorInterpolant{T, T, 2}(jd0 - J2000, sol.t, sol.p)
+end
 
-            _jd0, _tspan, _abstol = promote(jd0, tspan, abstol)
+function propagate(maxsteps::Int, jd0::T1, tspan::T2; dynamics::Function = NBP_pN_A_J23E_J23M_J2S_threads!,
+                    nast::Int = 343, order::Int = order, abstol::T3 = abstol, parse_eqs::Bool = true) where {T1, T2, T3 <: Real}
 
-            return propagate(maxsteps, _jd0, _tspan, $V_dense(); dynamics = dynamics, nast = nast, order = order,
-                             abstol = _abstol, parse_eqs = parse_eqs)
+    _jd0, _tspan, _abstol = promote(jd0, tspan, abstol)
 
-        end
+    return propagate(maxsteps, _jd0, _tspan; dynamics = dynamics, nast = nast, order = order,
+                        abstol = _abstol, parse_eqs = parse_eqs)
 
-    end
 end

test/propagation.jl

@@ -26,8 +26,6 @@ using LinearAlgebra: norm
     local jd0 = datetime2julian(DateTime(2000,1,1,12))
     # Number of years
     local nyears = 2031.0 - year(julian2datetime(jd0))
-    # Dense output
-    local dense = Val(true)
     # Dynamical function
     local dynamics = DE430!
 
@@ -68,7 +66,7 @@ using LinearAlgebra: norm
         @test iszero(zero(TaylorInterpolant{T, T, 2, SubArray{T, 1}, SubArray{Taylor1{T}, 2}}))
         @test iszero(zero(TaylorInterpolant{T, U, 2, SubArray{T, 1}, SubArray{Taylor1{U}, 2}}))
         # Test integration
-        sol = propagate(5, jd0, nyears, dense; dynamics=PlanetaryEphemeris.freeparticle!, order, abstol)
+        sol = propagate(5, jd0, nyears; dynamics=PlanetaryEphemeris.freeparticle!, order, abstol)
         @test sol isa TaylorInterpolant{Float64, Float64, 2}
         q0 = initialcond(N, jd0)
         @test sol(sol.t0) == q0
@@ -106,7 +104,7 @@ using LinearAlgebra: norm
         # Float64
 
         # Test integration
-        sol64 = propagate(100, jd0, nyears, dense; dynamics, order, abstol)
+        sol64 = propagate(100, jd0, nyears; dynamics, order, abstol)
         # Save solution
         filename = selecteph2jld2(sol64, bodyind, nyears)
         # Recovered solution
@@ -183,7 +181,7 @@ using LinearAlgebra: norm
         # Float 128
         #=
         # Test integration
-        sol128 = propagate(1, Float128(jd0), nyears, dense; dynamics = dynamics, order = order, abstol = abstol)
+        sol128 = propagate(1, Float128(jd0), nyears; dynamics = dynamics, order = order, abstol = abstol)
         # Save solution
         filename = selecteph2jld2(sol128, bodyind, nyears)
         # Recovered solution