Adapt changes in TI # 156 (refact taylorize) (#8)

* Refactor taylorize * Explicit taylorize code * Type-stability for propagate * Simplify propagate * Add tests * ArgParse for integrate_ephemeris.jl * Add constants * Update TS and TI, custom serializatoin, simplify propagate and precompilation * Fix an error in DE430! * Minor fixes * More minor fixes * Custom print for TaylorInterpolant * Minor fixes * Relative state vector callability methods * Bump minor version * Specify what we mean by backward integration * Minor fix * Update TaylorSeries and Abstract TaylorInterpolant fields * Fix join method for TaylorInterpolant * Compat version for JLD2 * Add ArgParse as a dependency * Remove number of threads from main
PerezHz · Apr 2, 2023 · 0af6c0d · 0af6c0d · PerezHz · Apr 2, 2023
1 parent c8f8b9c
commit 0af6c0d
Show file tree

Hide file tree

Showing 18 changed files with 14,276 additions and 5,589 deletions.
diff --git a/.gitignore b/.gitignore
@@ -1 +1,2 @@
-*.jld
+*.jld
+*.jld2
diff --git a/Project.toml b/Project.toml
@@ -1,25 +1,28 @@
 name = "PlanetaryEphemeris"
 uuid = "d83715d0-7e5f-11e9-1a59-4137b20d8363"
 authors = ["Jorge A. Pérez Hernández", "Luis Benet"]
-version = "0.2.1"
+version = "0.3.0"
 
 [deps]
+ArgParse = "c7e460c6-2fb9-53a9-8c5b-16f535851c63"
 AutoHashEquals = "15f4f7f2-30c1-5605-9d31-71845cf9641f"
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
 DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
-JLD = "4138dd39-2aa7-5051-a626-17a0bb65d9c8"
+JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Quadmath = "be4d8f0f-7fa4-5f49-b795-2f01399ab2dd"
+SnoopPrecompile = "66db9d55-30c0-4569-8b51-7e840670fc0c"
 TaylorIntegration = "92b13dbe-c966-51a2-8445-caca9f8a7d42"
 TaylorSeries = "6aa5eb33-94cf-58f4-a9d0-e4b2c4fc25ea"
 
 [compat]
+ArgParse = "1.1"
 AutoHashEquals = "0.2"
-JLD = "0.12"
+JLD2 = "0.4"
 Quadmath = "0.5"
-TaylorIntegration = "0.8"
-TaylorSeries = "0.12"
+TaylorIntegration = "0.11"
+TaylorSeries = "0.14"
 julia = "1.6"
 
 [extras]

diff --git a/README.md b/README.md
@@ -15,23 +15,22 @@ Instituto de Ciencias Físicas, Universidad Nacional Autónoma de México (UNAM)
 
 ## Installation
 
-The current development version of this package may be installed in Julia via:
+The current version of this package may be installed in Julia pkg manager via:
 ```
-import Pkg
-Pkg.add(Pkg.PackageSpec(url="https://github.com/PerezHz/PlanetaryEphemeris.jl.git", rev="main"))
+] add PlanetaryEphemeris
 ```
 
 ## Usage
 
-`PlanetaryEphemeris.propagate` is a high-level function which performs the
+`PlanetaryEphemeris.propagate_dense` is a high-level function which performs the
 numerical integration. The file `integrate_ephemeris.jl` in the `scripts` directory
 contains an example script. This script may be called as:
 
-`julia --project=@. integrate_ephemeris.jl`
+`julia --project integrate_ephemeris.jl --help`
 
 `PlanetaryEphemeris.propagate` also supports multi-threading:
 
-`JULIA_NUM_THREADS=<number-of-threads> julia --project=@. integrate_ephemeris.jl`
+`julia -t <number-of-threads> --project integrate_ephemeris.jl --help`
 
 ## Acknowledgments
 

diff --git a/scripts/integrate_ephemeris.jl b/scripts/integrate_ephemeris.jl
@@ -1,35 +1,113 @@
-#Multi-threaded:
-#JULIA_NUM_THREADS=<number-of-threads> julia --project=@. integrate_ephemeris.jl
-#Single-threaded:
-#julia --project=@. integrate_ephemeris.jl
-
-@show Threads.nthreads()
-
-using PlanetaryEphemeris
-using Dates
-
-#script parameters (TODO: use ArgParse.jl instead)
-const maxsteps = 1000000
-# jd0 = datetime2julian(DateTime(1969,6,28,0,0,0)) #starting time of integration
-const jd0 = datetime2julian(DateTime(2000,1,1,12)) #starting time of integration
-const nyears = 2031.0 - year(julian2datetime(jd0))
-@show jd0, J2000, jd0-J2000
-const dense = true #false
-const dynamics = DE430!
-@show dynamics
-const nast = 343 #16 # number of asteroid perturbers
-const quadmath = false #true # use quadruple precision
-###bodyind = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 17, 18, 19, 20, 21, 22, 23, 25, 27, 28, 30, 40, 41, 46, 55, 62, 73, 113, 115, 277, 322] # SS + 25 ast perturbers
-const bodyind = 1:(11+16) #1:(11+nast) # body indices in output
-
-# integration parameters
-const order = 25
-const abstol = 1.0E-20
-
-#integrator warmup
-PlanetaryEphemeris.propagate(1, jd0, nyears, output=false, dense=dense, dynamics=dynamics, nast=nast, quadmath=quadmath, bodyind=bodyind, order=order, abstol=abstol)
-println("*** Finished warmup")
-
-# perform full integration
-PlanetaryEphemeris.propagate(maxsteps, jd0, nyears, dense=dense, dynamics=dynamics, nast=nast, quadmath=quadmath, bodyind=bodyind, order=order, abstol=abstol)
-println("*** Finished full integration")
+using ArgParse, PlanetaryEphemeris, Dates
+
+function parse_commandline()
+
+    s = ArgParseSettings()
+
+    # Program name (for usage & help screen)
+    s.prog = "integrate_ephemeris.jl"  
+    # Desciption (for help screen)
+    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 
+        "--jd0"
+            help = "Starting time of integration; options are: \"2000-01-01T12:00:00\" or \"1969-06-28T00:00:00\""
+            arg_type = DateTime
+            default = DateTime(2000, 1, 1, 12) # DateTime(1969, 6, 28, 0, 0, 0)
+        "--nyears"
+            help = "Number of years"
+            arg_type = Float64
+            default = 31.0
+        "--dynamics"
+            help = "Dynamical model function"
+            arg_type = Function
+            default = DE430!
+        "--nast"
+            help = "Number of asteroid perturbers"
+            arg_type = Int
+            default = 343 # 16
+        "--order"
+            help = "Order of Taylor polynomials expansions during integration"
+            arg_type = Int
+            default = 25
+        "--abstol"
+            help = "Absolute tolerance"
+            arg_type = Float64
+            default = 1.0E-20
+        "--parse_eqs"
+            help = "Whether to use the taylorized method of jetcoeffs (a-priori faster) or not"
+            arg_type = Bool
+            default = true 
+        "--bodyind"
+            help = "Body indices in output"
+            arg_type = UnitRange{Int}
+            default = 1:(11+16) # 1:(11+nast)
+            # bodyind = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 17, 18, 19, 20, 21, 22, 23, 25, 27, 28, 30, 40, 41, 46, 55, 62, 73, 113, 115, 277, 322] # SS + 25 ast perturbers
+    end
+
+    s.epilog = """
+        examples:\n
+        \n
+        # Multi-threaded\n
+        julia -t 4 --project integrate_ephemeris.jl --maxsteps 100 --jd0 "2000-01-01T12:00:00"\n
+        \n
+        # Single-threaded\n
+        julia --project integrate_ephemeris.jl --maxsteps 100 --jd0 "2000-01-01T12:00:00"\n
+        \n
+    """
+
+    return parse_args(s)
+end
+
+function main(maxsteps::Int, jd0_datetime::DateTime, nyears::Float64, dynamics::Function, nast::Int,
+              bodyind::UnitRange{Int}, order::Int, abstol::Float64, parse_eqs::Bool)
+
+    println("*** Integrate Ephemeris ***")
+    println("Number of threads: ", Threads.nthreads())
+    println("Dynamical function: ", dynamics) 
+
+    jd0 = datetime2julian(jd0_datetime) 
+    println( "Initial time of integration: ", string(jd0_datetime) )
+    println( "Final time of integration: ", string(julian2datetime(jd0 + nyears*yr)) )
+
+    println("*** Integrator warmup ***")
+    _ = propagate(1, jd0, nyears, Val(true), dynamics = dynamics, nast = nast, order = order, abstol = abstol, 
+                  parse_eqs = parse_eqs)
+    println("*** Finished warmup ***")
+
+    println("*** Full integration ***")
+    sol = propagate(maxsteps, jd0, nyears, Val(true), dynamics = dynamics, nast = nast, order = order, abstol = abstol, 
+                    parse_eqs = parse_eqs)
+    println("*** Finished full integration ***")
+
+    # Total number of bodies (Sun + 8 Planets + Moon + Pluto + Asteroids)
+    N = 11 + nast
+
+    selecteph2jld2(sol, bodyind, nyears, N)
+
+    nothing 
+end 
+
+function main()
+
+    parsed_args = parse_commandline()
+
+    maxsteps = parsed_args["maxsteps"] :: Int
+    jd0_datetime = parsed_args["jd0"] :: DateTime
+    nyears = parsed_args["nyears"] :: Float64
+    dynamics = parsed_args["dynamics"] :: Function
+    nast = parsed_args["nast"] :: Int 
+    bodyind = parsed_args["bodyind"] :: UnitRange{Int}
+    order = parsed_args["order"] :: Int
+    abstol = parsed_args["abstol"] :: Float64
+    parse_eqs = parsed_args["parse_eqs"] :: Bool
+
+    main(maxsteps, jd0_datetime, nyears, dynamics, nast, bodyind, order, abstol, parse_eqs)
+
+end
+
+main()
diff --git a/src/PlanetaryEphemeris.jl b/src/PlanetaryEphemeris.jl
@@ -2,51 +2,25 @@ module PlanetaryEphemeris
 
 # __precompile__(false)
 
-export PE, au, yr, sundofs, earthdofs,
-    c_au_per_day, μ, NBP_pN_A_J23E_J23M_J2S!,
-    NBP_pN_A_J23E_J23M_J2S_threads!, DE430!,
-    semimajoraxis, eccentricity, inclination,
-    longascnode, argperi, longperi,
-    trueanomaly, ecanomaly, meananomaly,
-    timeperipass, lrlvec, eccentricanomaly,
-    meanan2truean, meanmotion, time2truean,
-    su, ea, mo, au, yr, daysec, clightkms,
-    c_au_per_day, c_au_per_sec, c_cm_per_sec,
-    J2000, R_sun, α_p_sun, δ_p_sun, au,
-    UJ_interaction, de430_343ast_ids, Rx, Ry, Rz,
-    ITM_und, ITM1, ITM2, R_moon, τ_M, k_2M,
-    JSEM, CM, SM, n1SEM, n2M, J2E, J2EDOT, RE,
-    k_20E, k_21E, k_22E, τ_0p, τ_1p, τ_2p, τ_0, τ_1, τ_2, ω_E, EMRAT,
-    TaylorInterpolant, selecteph2jld, ssb_posvel_pN, nbodyind
+export PE, au, yr, sundofs, earthdofs, c_au_per_day, μ, NBP_pN_A_J23E_J23M_J2S!, NBP_pN_A_J23E_J23M_J2S_threads!, DE430!,
+       semimajoraxis, eccentricity, inclination, longascnode, argperi, longperi, trueanomaly, ecanomaly, meananomaly,
+       timeperipass, lrlvec, eccentricanomaly, meanan2truean, meanmotion, time2truean, su, ea, mo, au, yr, daysec, clightkms,
+       c_au_per_day, c_au_per_sec, c_cm_per_sec, J2000, R_sun, α_p_sun, δ_p_sun, au, UJ_interaction, de430_343ast_ids, Rx, Ry, 
+       Rz, ITM_und, ITM1, ITM2, R_moon, τ_M, k_2M, JSEM, CM, SM, n1SEM, n2M, J2E, J2EDOT, RE, k_20E, k_21E, k_22E, τ_0p, τ_1p, 
+       τ_2p, τ_0, τ_1, τ_2, ω_E, EMRAT, TaylorInterpolant, selecteph2jld, ssb_posvel_pN, nbodyind, propagate, t2c_jpl_de430, 
+       c2t_jpl_de430, pole_rotation, selecteph2jld2, save2jld2andcheck, numberofbodies, kmsec2auday, auday2kmsec
 
 using AutoHashEquals
 using TaylorIntegration, LinearAlgebra
 using Printf
-using Dates: DateTime, julian2datetime, datetime2julian, year
+using Dates: DateTime, datetime2julian, year
 using DelimitedFiles
-using JLD
+using JLD2
 using Quadmath
 
-import Base.reverse
-
-@doc raw"""
-    nbodyind(N::Int, i::Int)
-    nbodyind(N::Int, ivec::AbstractVector{Int})
-
-Returns the indexes of the positions and velocities of the `i`-th body (or the 
-`ivec`-th bodies) in a vector with `N` bodies. The function assumes that the vector has 
-the form: `3N` positions + `3N` velocities (+ Lunar physical librations + TT-TDB). 
-"""
-nbodyind(N::Int, i::Int) = union(3i-2:3i, 3*(N+i)-2:3*(N+i))
-
-function nbodyind(N::Int, ivec::AbstractVector{Int})
-    a = Int[]
-    for i in ivec
-        i > N && continue
-        a = union(a, nbodyind(N,i))
-    end
-    return sort(a)
-end
+import Base: convert, reverse, show, join 
+import Dates: julian2datetime
+import JLD2: writeas
 
 include("constants.jl")
 include("jpl-de-430-431-earth-orientation-model.jl")
@@ -58,5 +32,6 @@ include("plephinteg.jl")
 include("propagation.jl")
 include("osculating.jl")
 include("barycenter.jl")
+#include("precompile.jl")
 
 end # module
diff --git a/src/barycenter.jl b/src/barycenter.jl
@@ -1,7 +1,7 @@
 @doc raw"""
     μ_star_fun(μ, q, i)
 
-Returns the mass parameter of the `i`-th body with relativistic corrections up to order ``1/c^2``
+Return the mass parameter of the `i`-th body with relativistic corrections up to order ``1/c^2``
 ```math
 \mu_i^* = \mu_i\left(1 + \frac{v_i^2}{2c^2} - \frac{1}{2c^2}\sum_{j\neq i}\frac{\mu_j}{r_{ij}}\right),
 ```
@@ -51,7 +51,7 @@ end
 @doc raw"""
     ssb_posvel_pN(μ, q)
 
-Returns the two sums needed to determine the Solar System Barycenter (SSB)
+Return the two sums needed to determine the Solar System Barycenter (SSB)
 ```math
 \left\lbrace
 \begin{array}{l}
@@ -102,7 +102,7 @@ end
 @doc raw"""
     sun_posvel_pN(μ, q)
 
-Solves
+Solve
 ```math
 \left\lbrace
 \begin{array}{l}

diff --git a/src/constants.jl b/src/constants.jl
@@ -1,11 +1,14 @@
 # PlanetaryEphemeris abbreviation
 const PE = PlanetaryEphemeris
 
+# Path to PlanetaryEphemeris src directory 
+const src_path = dirname(pathof(PlanetaryEphemeris))
+
 # Integration parameters
 const order = 30
 const abstol = 1.0E-20
 
-# Important bodies indexes 
+# Important bodies indices 
 const su = 1 # Sun's index
 const ea = 4 # Earth's index
 const mo = 5 # Moon's index
@@ -439,6 +442,30 @@ const c_cm_per_sec = 100_000*clightkms         # Speed of light in cm per sec
 const c_p2 = 29979.063823897606                # Speed of light^2 in au^2/day^2
 const c_m2 = 3.3356611996764786e-5             # Speed of light^-2 in day^2/au^2
 
+@doc raw"""
+    nbodyind(N::Int, i::Int)
+    nbodyind(N::Int, ivec::AbstractVector{Int})
+
+Return the indices of the positions and velocities of the `i`-th body (or the 
+`ivec`-th bodies) in a vector with `N` bodies. The function assumes that the vector has 
+the form: `3N` positions + `3N` velocities (+ Lunar physical librations + TT-TDB). 
+"""
+nbodyind(N::Int, i::Int) = union(3i-2:3i, 3*(N+i)-2:3*(N+i))
+
+function nbodyind(N::Int, ivec::T) where {T <: AbstractVector{Int}}
+    a = Vector{Int}(undef, 0)
+    for i in ivec
+        i > N && continue
+        a = union(a, nbodyind(N, i))
+    end
+
+    return sort(a)
+end
+
+numberofbodies(L::Int) = (L - 13) ÷ 6
+numberofbodies(v::Vector{T}) where {T} = numberofbodies(length(v))
+numberofbodies(m::Matrix{T}) where {T} = numberofbodies(size(m, 2))
+
 const sundofs = nbodyind(length(μ), su)        # Sun's position and velocity indices
 const earthdofs = nbodyind(length(μ), ea)      # Earth's position and velocity indices
 
@@ -598,6 +625,11 @@ const lnm5 = [2n-1 for n in 1:6]                         # (2n - 1)
 const lnm6 = [-(n+1) for n in 1:6]                       # -(n + 1)
 const lnm7 = [m for m in 1:6]                            # m  
 
+# Number of bodies in extended-body accelerations
+const N_ext = 11     
+# Number of bodies used to compute time-delayed tidal interactions
+const N_bwd = 11                   
+
 # Diagonal elements of undistorted lunar mantle moment of inertia
 # See equation (37) in page 16 of https://ui.adsabs.harvard.edu/abs/2014IPNPR.196C...1F%2F/abstract
 const A_T = (2(1-β_L*γ_L)/(2β_L-γ_L+β_L*γ_L))*J2_M_und
@@ -650,6 +682,10 @@ const τ_0 = 0.0                    # Rotational time-lag for long-period deform
 const τ_1 = 7.3632190228041890E-03 # Rotational time-lag for diurnal deformation, days
 const τ_2 = 2.5352978633388720E-03 # Rotational time-lag for semi-diurnal deformation, days
 
+# Tidal acceleration auxiliaries
+const μ_mo_div_μ_ea = μ[mo]/μ[ea]                # Ratio of Moon and Earth mass parameters
+const tid_num_coeff = 1.5*(1.0 + μ_mo_div_μ_ea)  # Overall numerical factor in equation (32)
+
 # Standard value of nominal mean angular velocity of Earth (rad/day), 
 # See Explanatory Supplement to the Astronomical Almanac 2014, section 7.4.3.3,
 # page 296: 7.2921151467e-5 rad/second