Update to PlanetaryEphemeris v0.7 (#25)

* Update Project.toml

* Update loadpeeph

* Update propagation methods and tests

* Update Artifacts.toml

* Update artifact generation script

* Update pha/apophis.jl

* Update Project.toml

* Update propagation.jl

* Remove remaining neosinteg method; use taylorinteg instead

* Fix in propagation.jl

* WIP: SatelliteToolboxTransformations

* Update Project.toml

* Perform ECEF->ECI transformations via SatelliteToolboxTransformations

* Update propagation tests

* Scattered updates

* Update test/propagation.jl

* Fix propagation tests

* Update propagation tests

---------

Co-authored-by: Luis Eduardo Ramírez Montoya <[email protected]>

Artifacts.toml

@@ -63,9 +63,9 @@ lazy = true
     url = "https://naif.jpl.nasa.gov/pub/naif/generic_kernels/lsk/naif0012.tls"
 
 [sseph_p100]
-git-tree-sha1 = "6fb49d939d8cc8675acbb8d64d9014e918c1569c"
+git-tree-sha1 = "d4d91f0045063719c9bf581bf20bbc212327b2d0"
 lazy = true
 
     [[sseph_p100.download]]
-    sha256 = "5051b43a30e33b3bde2388d7bf48dcdfefffd1dc3c139b37042b94b2f877baee"
+    sha256 = "25b525ffef5c05cc51ee554760ae83107b9bd163577439470b60c2d4c7153a21"
     url = "https://github.com/LuEdRaMo/sseph/raw/main/sseph343ast016_p100y_et.tar.gz"

Project.toml

@@ -1,7 +1,7 @@
 name = "NEOs"
 uuid = "b41c07a2-2abb-11e9-070a-c3c1b239e7df"
 authors = ["Jorge A. Pérez Hernández", "Luis Benet", "Luis Eduardo Ramírez Montoya"]
-version = "0.6.4"
+version = "0.7.0"
 
 [deps]
 Artifacts = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"
@@ -24,7 +24,7 @@ Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Quadmath = "be4d8f0f-7fa4-5f49-b795-2f01399ab2dd"
 Requires = "ae029012-a4dd-5104-9daa-d747884805df"
 SPICE = "5bab7191-041a-5c2e-a744-024b9c3a5062"
-SatelliteToolbox = "6ac157d9-b43d-51bb-8fab-48bf53814f4a"
+SatelliteToolboxTransformations = "6b019ec1-7a1e-4f04-96c7-a9db1ca5514d"
 Scratch = "6c6a2e73-6563-6170-7368-637461726353"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 TaylorIntegration = "92b13dbe-c966-51a2-8445-caca9f8a7d42"
@@ -48,15 +48,15 @@ Healpix = "4"
 IntervalArithmetic = "0.20"
 IntervalRootFinding = "0.5.11"
 JLD2 = "0.4"
-PlanetaryEphemeris = "0.6"
+PlanetaryEphemeris = "0.7"
 Quadmath = "0.5"
 Requires = "0.5.2, 1"
 SPICE = "0.2"
-SatelliteToolbox = "0.10"
+SatelliteToolboxTransformations = "0.1"
 Scratch = "1.2"
 StatsBase = "0.33, 0.34"
 Tables = "1.10"
-TaylorIntegration = "0.13"
+TaylorIntegration = "0.14"
 TaylorSeries = "0.15"
 julia = "1.6"
 

dev/generate_artifacts.jl

@@ -2,9 +2,10 @@
 # - debiasing tables from Chesley et al. (2010), Farnocchia et al. (2015) and Eggl et al. (2020)
 # - JPL lsk and spk
 # - PlanetaryEphemeris.jl Solar System ephemeris
-
-# TO DO: adapt script for julia1.6+
-# TO DO: adapt script for 2010 debiasing tables (requires extra unpacking)
+# TODO: adapt script for 2010 debiasing tables (requires extra unpacking)
+#
+# This script can be run from NEOs.jl root directory with:
+#      julia --project dev/generate_artifacts.jl
 
 using Pkg.Artifacts
 using Pkg.PlatformEngines: sha256, unpack

pha/apophis.jl

@@ -1,5 +1,4 @@
 using ArgParse, NEOs, PlanetaryEphemeris, Dates, TaylorIntegration, JLD2
-# using Revise, NEOs, PlanetaryEphemeris, Dates, TaylorIntegration, JLD2
 
 # Load JPL ephemeris
 loadjpleph()
@@ -45,10 +44,6 @@ function parse_commandline()
             help = "Whether to use the taylorized method of jetcoeffs or not"
             arg_type = Bool
             default = true
-        "--ss_eph_file"
-            help = "Path to local Solar System ephemeris file (loaded via artifact by default)"
-            arg_type = Union{String,Nothing}
-            default = nothing
     end
 
     s.epilog = """
@@ -65,15 +60,20 @@ function parse_commandline()
     return parse_args(s)
 end
 
-function print_header(header::String)
+function print_header(header::String, level::Int = 1)
     L = length(header)
-    println(repeat("-", L))
+    if level == 1
+        c = "="
+    else
+        c = "-"
+    end
+    println(repeat(c, L))
     println(header)
-    println(repeat("-", L))
+    println(repeat(c, L))
 end
 
-function main(dynamics::D, maxsteps::Int, jd0_datetime::DateTime, nyears_bwd::T, nyears_fwd::T,
-              ss16asteph_et::TaylorInterpolant, order::Int, varorder::Int, abstol::T, parse_eqs::Bool) where {T <: Real, D}
+function main(dynamics::D, maxsteps::Int, jd0_datetime::DateTime, nyears_bwd::T, nyears_fwd::T, order::Int, varorder::Int,
+              abstol::T, parse_eqs::Bool) where {T <: Real, D}
 
     # Initial conditions from Apophis JPL solution #197
     q00 = kmsec2auday(apophisposvel197(datetime2et(jd0_datetime)))
@@ -94,26 +94,29 @@ function main(dynamics::D, maxsteps::Int, jd0_datetime::DateTime, nyears_bwd::T,
     # Initial date (in julian days)
     jd0 = datetime2julian(jd0_datetime)
 
-    print_header("Integrator warmup")
-    # sol_bwd = NEOs.propagate(dynamics, 1, jd0, nyears_fwd, ss16asteph_et, q0, Val(true);
-    #                      order = order, abstol = abstol, parse_eqs = parse_eqs)
+    print_header("Integrator warmup", 2)
+    _ = NEOs.propagate(dynamics, 1, jd0, nyears_fwd, q0, Val(true);
+                       order = order, abstol = abstol, parse_eqs = parse_eqs)
+    println()
 
-    print_header("Main integration")
+    print_header("Main integration", 2)
     tmax = nyears_bwd*yr
     println("• Initial time of integration: ", string(jd0_datetime))
     println("• Final time of integration: ", julian2datetime(jd0 + tmax))
 
-    sol_bwd = NEOs.propagate(dynamics, maxsteps, jd0, nyears_bwd, ss16asteph_et, q0, Val(true);
-                         order = order, abstol = abstol, parse_eqs = parse_eqs)
-    PlanetaryEphemeris.save2jld2andcheck("Apophis_bwd.jld2", (asteph = sol_bwd,))
+    sol_bwd = NEOs.propagate(dynamics, maxsteps, jd0, nyears_bwd, q0, Val(true);
+                             order = order, abstol = abstol, parse_eqs = parse_eqs)
+    PE.save2jld2andcheck("Apophis_bwd.jld2", (asteph = sol_bwd,))
 
     tmax = nyears_fwd*yr
     println("• Initial time of integration: ", string(jd0_datetime))
     println("• Final time of integration: ", julian2datetime(jd0 + tmax))
 
-    sol_fwd = NEOs.propagate(dynamics, maxsteps, jd0, nyears_fwd, ss16asteph_et, q0, Val(true),
-                         order = order, abstol = abstol, parse_eqs = parse_eqs)
-    PlanetaryEphemeris.save2jld2andcheck("Apophis_fwd.jld2", (asteph = sol_fwd,))
+    sol_fwd = NEOs.propagate(dynamics, maxsteps, jd0, nyears_fwd, q0, Val(true);
+                             order = order, abstol = abstol, parse_eqs = parse_eqs)
+    PE.save2jld2andcheck("Apophis_fwd.jld2", (asteph = sol_fwd,))
+
+    println()
 
     # NEO
     # Change t, x, v units, resp., from days, au, au/day to sec, km, km/sec
@@ -122,12 +125,12 @@ function main(dynamics::D, maxsteps::Int, jd0_datetime::DateTime, nyears_bwd::T,
     xva(et) = bwdfwdeph(et, sol_bwd, sol_fwd)
     # Earth
     # Change x, v units, resp., from au, au/day to km, km/sec
-    eph_ea = selecteph(ss16asteph_et, ea)
-    xve(et) = auday2kmsec(eph_ea(et))
+    eph_ea = selecteph(NEOs.sseph, ea)
+    xve(et) = auday2kmsec(eph_ea(et/daysec))
     # Sun
     # Change x, v units, resp., from au, au/day to km, km/sec
-    eph_su = selecteph(ss16asteph_et, su)
-    xvs(et) = auday2kmsec(eph_su(et))
+    eph_su = selecteph(NEOs.sseph, su)
+    xvs(et) = auday2kmsec(eph_su(et/daysec))
 
 
     radec_2004_2020 = read_radec_mpc("/Users/Jorge/projects/NEOs/data/99942_2004_2020.dat")
@@ -140,7 +143,7 @@ function main(dynamics::D, maxsteps::Int, jd0_datetime::DateTime, nyears_bwd::T,
     # deldop = vcat(deldop_2005_2013,deldop_2021)
 
     # Compute optical residuals
-    print_header("Compute residuals")
+    print_header("Compute residuals", 2)
     res_radec, w_radec = residuals(radec, xvs = xvs, xve = xve, xva = xva)
 
     # Compute radar residuals
@@ -149,7 +152,7 @@ function main(dynamics::D, maxsteps::Int, jd0_datetime::DateTime, nyears_bwd::T,
     res_del_bwd, w_del_bwd, res_dop_bwd, w_dop_bwd = residuals(deldop_2005_2013, xvs = xvs, xve = xve, xva = xva_bwd, niter=5)
     res_del_fwd, w_del_fwd, res_dop_fwd, w_dop_fwd = residuals(deldop_2021, xvs = xvs, xve = xve, xva = xva_fwd, niter=5)
 
-    PlanetaryEphemeris.save2jld2andcheck("resw_radec.jld2", (;res_radec,w_radec,res_del_bwd,w_del_bwd,res_dop_bwd,w_dop_bwd,res_del_fwd,w_del_fwd,res_dop_fwd,w_dop_fwd))
+    PE.save2jld2andcheck("resw_radec.jld2", (;res_radec,w_radec,res_del_bwd,w_del_bwd,res_dop_bwd,w_dop_bwd,res_del_fwd,w_del_fwd,res_dop_fwd,w_dop_fwd))
 
     @show NEOs.cte(res) NEOs.cte(w)
 
@@ -163,18 +166,16 @@ function main()
     parsed_args = parse_commandline()
 
     print_header("Asteroid Apophis")
-    print_header("General parameters")
+    println()
+    print_header("General parameters", 2)
 
     # Number of threads
     N_threads = Threads.nthreads()
     println("• Number of threads: ", N_threads)
 
     # Dynamical function
-    if N_threads == 1
-        dynamics = RNp1BP_pN_A_J23E_J2S_ng_eph!
-    else
-        dynamics = RNp1BP_pN_A_J23E_J2S_ng_eph_threads!
-    end
+    dynamics = RNp1BP_pN_A_J23E_J2S_ng_eph_threads!
+
     println("• Dynamical function: ", dynamics)
 
     # Maximum number of steps
@@ -190,12 +191,6 @@ function main()
     # Number of years in forward integration
     nyears_fwd = parsed_args["nyears_fwd"]
 
-    # Solar system ephemeris
-    print("• Loading Solar System ephemeris... ")
-    ss_eph_file = parsed_args["ss_eph_file"]
-    ss16asteph_et = ss_eph_file === nothing ? NEOs.sseph : JLD2.load(ss_eph_file, "ss16asteph_et")
-    println("Done")
-
     # Order of Taylor polynomials
     order = parsed_args["order"]
     println("• Order of Taylor polynomials: ", order)
@@ -210,10 +205,9 @@ function main()
 
     # Whether to use @taylorize or not
     parse_eqs = parsed_args["parse_eqs"]
-    println("• Use @taylorize: ", parse_eqs)
+    println("• Use @taylorize: ", parse_eqs, "\n")
 
-    main(dynamics, maxsteps, jd0_datetime,  nyears_bwd, nyears_fwd, ss16asteph_et,
-         order, varorder, abstol, parse_eqs)
+    main(dynamics, maxsteps, jd0_datetime,  nyears_bwd, nyears_fwd, order, varorder, abstol, parse_eqs)
 end
 
 main()

scripts/distributed.jl

@@ -19,7 +19,7 @@ end
     const dense = true#false
     const apophisjlpath = pkgdir(Apophis)
     const radarobsfile = joinpath(apophisjlpath, "Apophis_JPL_data_2012_2013.dat")
-    const dynamics = RNp1BP_pN_A_J23E_J2S_ng_eph!
+    const dynamics = RNp1BP_pN_A_J23E_J2S_ng_eph_threads!
     const t0 = datetime2julian(DateTime(2008,9,24,0,0,0)) #starting time of integration
     const tmax = t0+365.25nyears #final time of integration
     @show t0 == 2454733.5

scripts/main.jl

@@ -19,7 +19,6 @@ const dense = false #true
 const quadmath = false # use quadruple precision
 const debias_table = "2018" # "2014", "hires2018"
 const apophisjlpath = pkgdir(NEOs)
-# const dynamics = RNp1BP_pN_A_J23E_J2S_ng_eph!
 const dynamics = RNp1BP_pN_A_J23E_J2S_ng_eph_threads!
 const jd0 = datetime2julian(DateTime(2008,9,24,0,0,0)) #Julian date of integration initial time
 @show jd0 == 2454733.5

src/NEOs.jl

@@ -19,10 +19,9 @@ using PlanetaryEphemeris: daysec, su, ea, α_p_sun, δ_p_sun, t2c_jpl_de430, pol
       nbodyind, ordpres_differentiate, numberofbodies, kmsec2auday, auday2kmsec, meanmotion,
       meananomaly
 using Healpix: ang2pixRing, Resolution
-using SatelliteToolbox: get_iers_eop_iau_2000A, EOPData_IAU1980, EOPData_IAU2000A, JD_J2000,
-      orbsv, sv_ecef_to_eci, get_Δat, nutation_fk5, r_ecef_to_eci, T_ECIs, T_ECIs_IAU_2006,
-      we, OrbitStateVector, r_ecef_to_eci, DCM
-import SatelliteToolbox.sv_ecef_to_eci
+using SatelliteToolboxTransformations
+import SatelliteToolboxTransformations.sv_ecef_to_eci
+import SatelliteToolboxTransformations.sv_ecef_to_ecef
 using Dates: format
 using Downloads: download
 import Downloads
@@ -52,15 +51,15 @@ export loadjpleph, sunposvel, earthposvel, moonposvel, apophisposvel197, apophis
 # Osculating
 export pv2kep, yarkp2adot
 # Topocentric
-export obs_pos_ECEF, obsposvelECI, t2c_rotation_iau_76_80
+export obs_pos_ECEF, obsposvelECI
 # Process radec
 export compute_radec, debiasing, w8sveres17, radec_astrometry, residuals
 # Process radar
 export compute_delay, radar_astrometry
 # Gauss method
 export gauss_method
 # Asteroid dynamical models
-export RNp1BP_pN_A_J23E_J2S_ng_eph!, RNp1BP_pN_A_J23E_J2S_ng_eph_threads!, RNp1BP_pN_A_J23E_J2S_eph_threads!
+export RNp1BP_pN_A_J23E_J2S_ng_eph_threads!, RNp1BP_pN_A_J23E_J2S_eph_threads!
 # Propagate
 export propagate, propagate_lyap, propagate_root, save2jldandcheck
 

src/observations/jpl_eph.jl

@@ -115,32 +115,31 @@ See also [`getposvel`](@ref).
 """
 dtt_tdb(et) = getposvel(1000000001, 1000000000, cte(et))[4] # units: seconds/seconds
 
-# Load Solar System 2000-2100 ephemeris
+# Load Solar System, accelerations, newtonian potentials and TT-TDB 2000-2100 ephemeris
 const sseph_artifact_path = joinpath(artifact"sseph_p100", "sseph343ast016_p100y_et.jld2")
 const sseph::TaylorInterpolant{Float64, Float64, 2} = JLD2.load(sseph_artifact_path, "ss16ast_eph")
+const acceph::TaylorInterpolant{Float64, Float64, 2} = JLD2.load(sseph_artifact_path, "acc_eph")
+const poteph::TaylorInterpolant{Float64, Float64, 2} = JLD2.load(sseph_artifact_path, "pot_eph")
 const ttmtdb::TaylorInterpolant{Float64, Float64, 1} = TaylorInterpolant(sseph.t0, sseph.t, sseph.x[:,end])
 
 @doc raw"""
-    loadpeeph()
-    loadpeeph(et::Real)
-    loadpeeph(et_0::Real, et_f::Real)
+    loadpeeph(eph::TaylorInterpolant{Float64, Float64, 2} = sseph, t_0::Real = 0.0, t_f::Real = 36525.0)
 
-Load Solar System ephemeris produced by `PlanetaryEphemeris.jl` in timerange `[0, et]` (`[et_0, et_f]`) where `et` must have units
-of ephemeris seconds since J2000. If no `et` is given, return the full (100 years) integration.
+Load ephemeris produced by `PlanetaryEphemeris.jl` in timerange `[t_0, t_f] ⊆ [0.0, 36525.0]` where `t` must have units of TDB 
+days since J2000. The available options for `eph` are:
 
-**Caution**: running this function for the first time will download the `sseph_p100` artifact (∼556 MB) which can take several minutes.
-"""
-loadpeeph()::TaylorInterpolant{Float64, Float64, 2} = sseph
-
-function loadpeeph(et::Real)
-    i = searchsortedfirst(sseph.t, et)
-    return TaylorInterpolant(sseph.t0, sseph.t[1:i], sseph.x[1:i-1, :])
-end
+- `NEOs.sseph`: Solar system ephemeris. 
+- `NEOs.acceph`: accelerations ephemeris. 
+- `NEOs.poteph`: newtonian potentials ephemeris.
 
-function loadpeeph(et_0::Real, et_f::Real)
-    i_0 = searchsortedlast(sseph.t, et_0)
-    i_f = searchsortedfirst(sseph.t, et_f)
-    return TaylorInterpolant(sseph.t0, sseph.t[i_0:i_f], sseph.x[i_0:i_f-1, :])
+!!! warning
+    Running this function for the first time will download the `sseph_p100` artifact (885 MB) which can take several minutes.
+"""
+function loadpeeph(eph::TaylorInterpolant{Float64, Float64, 2} = sseph, t_0::Real = sseph.t0, t_f::Real = sseph.t0 + sseph.t[end])
+    @assert 0.0 ≤ t_0 ≤ t_f ≤ 36525.0
+    i_0 = searchsortedlast(eph.t, t_0)
+    i_f = searchsortedfirst(eph.t, t_f)
+    return TaylorInterpolant(eph.t0, eph.t[i_0:i_f], eph.x[i_0:i_f-1, :])
 end
 
 @doc raw"""

src/observations/process_radec.jl

@@ -294,7 +294,7 @@ function debiasing(obs::RadecMPC{T}, mpc_catalogue_codes_201X::Vector{String}, t
         # Seconds since J2000 (TDB)
         et_secs_i = datetime2et(obs.date)
         # Seconds sinde J2000 (TT)
-        tt_secs_i = et_secs_i - ttmtdb(et_secs_i)
+        tt_secs_i = et_secs_i - ttmtdb(et_secs_i/daysec)
         # Years since J2000
         yrs_J2000_tt = tt_secs_i/(daysec*yr)
         # Total debiasing correction in right ascension (arcsec)