Minor scattered updates (#33)

* propagate_root with dense output now returns a TaylorInterpolant; update tests

* Update Apophis script

* Update apophis.jl

* Update pha/Project.toml

* Fix bug in Valsecchi circles computation

* Update jetcoeffs.jl

* Bump patch version

* Fix obs_mpc_url

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.7.1"
+version = "0.7.2"
 
 [deps]
 Artifacts = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

pha/Project.toml

@@ -1,5 +1,4 @@
 [deps]
-
 ArgParse = "c7e460c6-2fb9-53a9-8c5b-16f535851c63"
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
@@ -15,5 +14,4 @@ StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 TaylorIntegration = "92b13dbe-c966-51a2-8445-caca9f8a7d42"
 
 [compat]
-
-NEOs = "0.7"
+NEOs = "0.7"

pha/apophis.jl

@@ -125,7 +125,9 @@ function main(dynamics::D, maxsteps::Int, jd0_datetime::DateTime, nyears_bwd::T,
 
     print_header("Integrator warmup", 2)
     _ = NEOs.propagate(dynamics, 1, jd0, nyears_fwd, q0, Val(true);
-                       order = order, abstol = abstol, parse_eqs = parse_eqs)
+                       order, abstol, parse_eqs)
+    _ = NEOs.propagate_root(dynamics, 1, jd0, nyears_fwd, q0, Val(true);
+                       order, abstol, parse_eqs)
     println()
 
     print_header("Main integration", 2)
@@ -136,16 +138,16 @@ function main(dynamics::D, maxsteps::Int, jd0_datetime::DateTime, nyears_bwd::T,
     sol_bwd = NEOs.propagate(dynamics, maxsteps, jd0, nyears_bwd, q0, Val(true);
                              order, abstol, parse_eqs)
     jldsave("Apophis_bwd.jld2"; sol_bwd)
-    # sol_bwd = JLD2.load("Apophis_bwd.jld2", "asteph")
+    # sol_bwd = JLD2.load("Apophis_bwd.jld2", "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, q0, Val(true);
+    sol_fwd, tvS, xvS, gvS = NEOs.propagate_root(dynamics, maxsteps, jd0, nyears_fwd, q0, Val(true);
                              order, abstol, parse_eqs)
-    jldsave("Apophis_fwd.jld2"; sol_fwd)
-    # sol_fwd = JLD2.load("Apophis_fwd.jld2", "asteph")
+    jldsave("Apophis_fwd.jld2"; sol_fwd, tvS, xvS, gvS)
+    # sol_fwd = JLD2.load("Apophis_fwd.jld2", "sol_bwd")
     println()
 
     # load Solar System ephemeris
@@ -181,7 +183,7 @@ function main(dynamics::D, maxsteps::Int, jd0_datetime::DateTime, nyears_bwd::T,
 
     # Compute radar residuals
     res_del, w_del, res_dop, w_dop = NEOs.residuals(deldop; xvs, xve, xva, niter=10, tord=10)
-    jldsave("Apophis_res_w_radec.jld2"; res_del, w_del, res_dop, w_dop)
+    jldsave("Apophis_res_w_deldop.jld2"; res_del, w_del, res_dop, w_dop)
     # JLD2.@load "Apophis_res_w_deldop.jld2"
 
     ### Process optical astrometry (filter, weight, debias)

src/constants.jl

@@ -1,20 +1,20 @@
-# Internal paths 
+# Internal paths
 
-# Path to NEOs src directory 
+# Path to NEOs src directory
 const src_path = dirname(pathof(NEOs))
-# Path to scratch space 
+# Path to scratch space
 const scratch_path = Ref{String}("")
 
 # URLs
 
-# MPC catalogues file url 
+# MPC catalogues file url
 const mpc_catalogues_url = "https://minorplanetcenter.net/iau/info/CatalogueCodes.html"
-# MPC observatories file url 
+# MPC observatories file url
 const mpc_observatories_url = "https://minorplanetcenter.net/iau/lists/ObsCodes.html"
-# MPC database search url 
+# MPC database search url
 const search_mpc_url = "https://www.minorplanetcenter.net/db_search/show_object?utf8=%E2%9C%93&object_id="
-# MPC observations url 
-const obs_mpc_url = "https://www.minorplanetcenter.net/tmp/"
+# MPC observations url
+const obs_mpc_url = "https://www.minorplanetcenter.net/tmp2/"
 
 # Abbreviations
 const cte = constant_term
@@ -28,31 +28,31 @@ const μ_DE430 = PE.μ
 const μ_B16_DE430 = μ_DE430[12:27]     # DE430 GM's of 16 most massive asteroids
 const μ_ast343_DE430 = μ_DE430[12:end] # DE430 GM's of 343 main belt asteroids included in DE430 integration
 # Gravitational parameter of the Sun [au^2 / day^3]
-const μ_S = PE.GMS  
+const μ_S = PE.GMS
 
 # Standard value of nominal mean angular velocity of Earth (rad/sec)
 # See Explanatory Supplement to the Astronomical Almanac 2014 Sec 7.4.3.3 p. 296
 const ω = 7.2921151467e-5 # 7.292115e-5 rad/sec
 
 # Solar corona parameters
-# See Table 8.5 in page 329 of Explanatory Supplement to the Astronomical Almanac 2014 
+# See Table 8.5 in page 329 of Explanatory Supplement to the Astronomical Almanac 2014
 const A_sun = 1.06e8  # [cm^-3]
 const a_sun = 4.89e5  # [cm^-3]
 const b_sun = 3.91e5  # [cm^-3]
 
 # Sun radiated power intensity at photosphere surface, Watt/meter^2
-const S0_sun = 63.15E6 
+const S0_sun = 63.15E6
 # Conversion factor from m^2/sec^3 to au^2/day^3
-const m2_s3_to_au2_day3 = 1e-6daysec^3/au^2 
+const m2_s3_to_au2_day3 = 1e-6daysec^3/au^2
 
 # Vector of J_2*R^2 values
 # J_2: second zonal harmonic coefficient
-# R: radius of the body 
+# R: radius of the body
 const Λ2 = zeros(11)
 Λ2[ea] = 1.9679542578489185e-12  # Earth
 # Vector of J_3*R^3 values
 # J_3: third zonal harmonic coefficient
-# R: radius of the body 
+# R: radius of the body
 const Λ3 = zeros(11)
 Λ3[ea] = -1.962633335678878e-19  # Earth
 
@@ -67,6 +67,6 @@ const d_EM_km = 384_400
 # Earth-Moon distance in [au]
 const d_EM_au = 384_400 / au
 
-# Zeroth order obliquity of the ecliptic in degrees 
+# Zeroth order obliquity of the ecliptic in degrees
 # See equation (5-153) in page 5-61 of https://doi.org/10.1002/0471728470.
 const ϵ0_deg = 84381.448/3_600

src/postprocessing/b_plane.jl

@@ -5,14 +5,14 @@ Returns the "critical" ``B``, derived from conservation of energy and angular mo
 ```math
 B = \sqrt{1 + \frac{2\mu_P}{R_P v_\infty^2}},
 ```
-i.e., what impact parameter ``B`` corresponds to a grazing impact in hyperbolic close 
-encounter. ``\mu_P`` is the planet's gravitational parameter, ``R_P`` is the planet's 
-radius and ``v_\infty`` is the asymptotic inbound velocity. If actual ``B`` is equal or 
-less to this, then impact happens. Output is in planet radii. 
+i.e., what impact parameter ``B`` corresponds to a grazing impact in hyperbolic close
+encounter. ``\mu_P`` is the planet's gravitational parameter, ``R_P`` is the planet's
+radius and ``v_\infty`` is the asymptotic inbound velocity. If actual ``B`` is equal or
+less to this, then impact happens. Output is in planet radii.
 
 See equations (13)-(14) in pages 4-5 of https://doi.org/10.1007/s10569-019-9914-4.
 
-# Arguments 
+# Arguments
 
 - `μ_P`: planetary gravitational parameter (au^3/day^2).
 - `R_P`: planetary radius (au).
@@ -36,7 +36,7 @@ close encounter. Returns a named tuple with the following fields:
 
 - `b` = ``b_E``, where ``b_E`` is the Earth impact cross section ("critical B"). See [`crosssection`](@ref).
 
-# Arguments 
+# Arguments
 
 - `xae`: asteroid's geocentric position/velocity vector at closest approach in au, au/day.
 - `xes`: planet's heliocentric position/velocity vector at asteroid's closest approach in au, au/day.
@@ -86,7 +86,7 @@ function bopik(xae, xes)
 
     # Computation of Öpik's coordinates of impact parameter vector \mathbf{B}
 
-    # B-vector: "vector from the planet center to the intersection between the B-plane 
+    # B-vector: "vector from the planet center to the intersection between the B-plane
     # and the asymptote".
     # See equation (4) in page 3 of https://doi.org/10.1007/s10569-019-9914-4
     Bvec = cross(S_v, hvec)./v_infty
@@ -95,23 +95,16 @@ function bopik(xae, xes)
     B_dot_ζ = dot(Bvec, ζ_v)
 
     # Computation of planetocentric velocity (U) vector
-    # res = xes[1:3] # Earth's heliocentric position, au
-    # res_norm = norm(res) # Earth's heliocentric range, au
-    # res_unit = res/res_norm # Earth's heliocentric radius unit vector
-    # @show res_norm
-    ves = v_pl*au/daysec # Earth's velocity, km/s
-    ves_norm = norm(ves) # Earth's speed, km/s
+    ves = v_pl # Earth's velocity, au/day
+    ves_norm = sqrt(ves[1]^2 + ves[2]^2 + ves[3]^2) # Earth's speed, au/day
     ves_unit = ves/ves_norm # Earth's velocity unit vector
-    # X-Y angle (degrees)
-    # @show rad2deg(acos(dot(ves_unit, res_unit)))
     # angle between Y-axis and \vec{U}
     cosθ = dot(S_v, ves_unit)
-    # @show cosθ()
     v_infty_kms = v_infty*au/daysec # Asteroid unperturbed speed, km/sec
     # @show v_infty_kms()
     # The norm of \vec{U} in appropriate units
-    U_unit = (2pi*au)/(yr*daysec) # 1 U in km/s
-    U_norm = v_infty_kms/U_unit
+    U_unit = ves_norm # 1 U = v_ast/v_pl
+    U_norm = v_infty/U_unit
     # U_y
     U_y = U_norm*cosθ
     # @show U_y U_norm
@@ -124,18 +117,18 @@ function bopik(xae, xes)
 end
 
 @doc raw"""
-    valsecchi_circle(a, e, i, k, h; m_pl=3.003489614915764e-6)  
+    valsecchi_circle(a, e, i, k, h; m_pl=3.003489614915764e-6)
 
 Computes Valsecchi circle associated to a mean motion resonance. Returns radius ``R`` (au) and
 ``\zeta``-axis coordinate ``D`` (au). This function first computes the Y-component and norm
 of the planetocentric velocity vector ``\mathbf{U}``
 ```math
-U_y = \sqrt{a(1-e^2)}\cos i - 1 \quad \text{and} \quad 
+U_y = \sqrt{a(1-e^2)}\cos i - 1 \quad \text{and} \quad
 U = ||\mathbf{U}|| = \sqrt{3 - \frac{1}{a} - 2\sqrt{a(1-e^2)}\cos i},
 ```
 and then substitutes into `valsecchi_circle(U_y, U_norm, k, h; m_pl=3.003489614915764e-6, a_pl=1.0)`.
-`a `, `e` and `i` are the asteroid heliocentric semimajor axis (au), eccentricity and 
-inclination (rad) respectively. 
+`a `, `e` and `i` are the asteroid heliocentric semimajor axis (au), eccentricity and
+inclination (rad) respectively.
 
 See section 2.1 in page 1181 of https://doi.org/10.1051/0004-6361:20031039.
 
@@ -169,13 +162,13 @@ Computes Valsecchi circle associated to a mean motion resonance. Returns radius
 R = \left|\frac{c\sin\theta_0'}{\cos\theta_0' - \cos\theta}\right| \quad \text{and} \quad
 D = \frac{c\sin\theta}{\cos\theta_0' - \cos\theta},
 ```
-where ``c = m/U^2`` with ``m`` the mass of the planet and ``U = ||\mathbf{U}||`` the norm of 
+where ``c = m/U^2`` with ``m`` the mass of the planet and ``U = ||\mathbf{U}||`` the norm of
 the planetocentric velocity vector; and ``\theta``, ``\theta_0'`` are the angles between Y-axis
-and ``\mathbf{U}`` pre and post encounter respectively. 
+and ``\mathbf{U}`` pre and post encounter respectively.
 
 See pages 1181, 1182 and 1187 of https://doi.org/10.1051/0004-6361:20031039.
 
-# Arguments 
+# Arguments
 - `U_y`: Y-component of unperturbed planetocentric velocity (Y-axis coincides with the direction of motion of the planet).
 - `U_norm`: Euclidean norm of unperturbed planetocentric velocity. Both `U_y`, `U_norm` are in units such that the heliocentric velocity of the planet is 1.
 - `k/h`: `h` heliocentric revolutions of asteroid per `k` heliocentric revolutions of Earth.
@@ -197,8 +190,8 @@ function valsecchi_circle(U_y, U_norm, k, h; m_pl=3.003489614915764e-6, a_pl=1.0
     # See page 1187 of https://doi.org/10.1051/0004-6361:20031039
     cosθ0p = (1-(U_norm^2)-(1/a0p))/(2U_norm)
     sinθ0p = sin(acos(cosθ0p))
-    # c = m/U^2 
-    # See first sentence below equation (3) in section 2.3, page 1182 of 
+    # c = m/U^2
+    # See first sentence below equation (3) in section 2.3, page 1182 of
     # https://doi.org/10.1051/0004-6361:20031039
     c = m_pl/(U_norm^2)
     # Radius of the Valsecchi circle R and its ζ-axis component D

src/propagation/jetcoeffs.jl

@@ -138,7 +138,7 @@ function TaylorIntegration._allocate_jetcoeffs!(::Val{RNp1BP_pN_A_J23E_J2S_ng_ep
     accY = Taylor1(identity(constant_term(zero_q_1)), order)
     accZ = Taylor1(identity(constant_term(zero_q_1)), order)
     local M_ = Array{S}(undef, 3, 3, N)
-    local M_[:, :, ea] = t2c_jpl_de430(dsj2k)
+    local M_[:, :, ea] = t2c_jpl_de430(dsj2k) .+ zero_q_1
     dq[1] = Taylor1(identity(constant_term(q[4])), order)
     dq[2] = Taylor1(identity(constant_term(q[5])), order)
     dq[3] = Taylor1(identity(constant_term(q[6])), order)
@@ -810,7 +810,7 @@ function TaylorIntegration.jetcoeffs!(::Val{RNp1BP_pN_A_J23E_J2S_ng_eph_threads!
     accZ.coeffs[1] = identity(constant_term(zero_q_1))
     accZ.coeffs[2:order + 1] .= zero(accZ.coeffs[1])
     local M_ = Array{S}(undef, 3, 3, N)
-    local M_[:, :, ea] = t2c_jpl_de430(dsj2k)
+    local M_[:, :, ea] = t2c_jpl_de430(dsj2k) .+ zero_q_1
     (dq[1]).coeffs[1] = identity(constant_term(q[4]))
     (dq[1]).coeffs[2:order + 1] .= zero((dq[1]).coeffs[1])
     (dq[2]).coeffs[1] = identity(constant_term(q[5]))
@@ -1749,7 +1749,7 @@ function TaylorIntegration._allocate_jetcoeffs!(::Val{RNp1BP_pN_A_J23E_J2S_eph_t
     accY = Taylor1(identity(constant_term(zero_q_1)), order)
     accZ = Taylor1(identity(constant_term(zero_q_1)), order)
     local M_ = Array{S}(undef, 3, 3, N)
-    local M_[:, :, ea] = t2c_jpl_de430(dsj2k)
+    local M_[:, :, ea] = t2c_jpl_de430(dsj2k) .+ zero_q_1
     dq[1] = Taylor1(identity(constant_term(q[4])), order)
     dq[2] = Taylor1(identity(constant_term(q[5])), order)
     dq[3] = Taylor1(identity(constant_term(q[6])), order)
@@ -2324,7 +2324,7 @@ function TaylorIntegration.jetcoeffs!(::Val{RNp1BP_pN_A_J23E_J2S_eph_threads!},
     accZ.coeffs[1] = identity(constant_term(zero_q_1))
     accZ.coeffs[2:order + 1] .= zero(accZ.coeffs[1])
     local M_ = Array{S}(undef, 3, 3, N)
-    local M_[:, :, ea] = t2c_jpl_de430(dsj2k)
+    local M_[:, :, ea] = t2c_jpl_de430(dsj2k) .+ zero_q_1
     (dq[1]).coeffs[1] = identity(constant_term(q[4]))
     (dq[1]).coeffs[2:order + 1] .= zero((dq[1]).coeffs[1])
     (dq[2]).coeffs[1] = identity(constant_term(q[5]))

src/propagation/propagation.jl

@@ -162,12 +162,12 @@ for V_dense in V_true_false
                            parse_eqs::Bool = true) where {T <: Real, U <: Number, D}
 
             # Parameters for taylorinteg
-            _q0, _t0, _tmax, _params = propagate_params(jd0, tspan, q0; μ_ast = μ_ast, order = order, abstol = abstol)
+            _q0, _t0, _tmax, _params = propagate_params(jd0, tspan, q0; μ_ast, order, abstol)
 
             # Propagate orbit
 
             @time sol = taylorinteg(dynamics, _q0, _t0, _tmax, order, abstol, $V_dense(), _params;
-                                    maxsteps = maxsteps, parse_eqs = parse_eqs)
+                                    maxsteps, parse_eqs)
 
             # Dense output (save Taylor polynomials in each step)
             if $V_dense == Val{true}
@@ -186,20 +186,20 @@ for V_dense in V_true_false
 
             _jd0, _tspan, _abstol = promote(jd0, tspan, abstol)
 
-            return propagate(dynamics, maxsteps, _jd0, _tspan, q0, $V_dense(); μ_ast = μ_ast, order = order,
-                             abstol = _abstol, parse_eqs = parse_eqs)
+            return propagate(dynamics, maxsteps, _jd0, _tspan, q0, $V_dense(); μ_ast, order,
+                             abstol = _abstol, parse_eqs)
         end
 
         function propagate(dynamics::D, maxsteps::Int, jd0::T, nyears_bwd::T, nyears_fwd::T, q0::Vector{U}, ::$V_dense;
                            μ_ast::Vector = μ_ast343_DE430[1:end], order::Int = order, abstol::T = abstol,
                            parse_eqs::Bool = true) where {T <: Real, U <: Number, D}
 
             # Backward integration
-            bwd = propagate(dynamics, maxsteps, jd0, nyears_bwd, q0, $V_dense(); μ_ast = μ_ast, order = order,
-                            abstol = abstol, parse_eqs = parse_eqs)
+            bwd = propagate(dynamics, maxsteps, jd0, nyears_bwd, q0, $V_dense(); μ_ast, order,
+                            abstol, parse_eqs)
             # Forward integration
-            fwd = propagate(dynamics, maxsteps, jd0, nyears_fwd, q0, $V_dense(); μ_ast = μ_ast, order = order,
-                            abstol = abstol, parse_eqs = parse_eqs)
+            fwd = propagate(dynamics, maxsteps, jd0, nyears_fwd, q0, $V_dense(); μ_ast, order,
+                            abstol, parse_eqs)
 
             return bwd, fwd
 
@@ -210,13 +210,20 @@ for V_dense in V_true_false
                                 order::Int = order, abstol::T = abstol) where {T <: Real, U <: Number, D}
 
             # Parameters for neosinteg
-            _q0, _t0, _tmax, _params = propagate_params(jd0, tspan, q0; μ_ast = μ_ast, order = order, abstol = abstol)
+            _q0, _t0, _tmax, _params = propagate_params(jd0, tspan, q0; μ_ast, order, abstol)
 
             # Propagate orbit
             @time sol = taylorinteg(dynamics, rvelea, _q0, _t0, _tmax, order, abstol, $V_dense(), _params;
                                   maxsteps, parse_eqs, eventorder, newtoniter, nrabstol)
 
-            return sol
+            # Dense output (save Taylor polynomials in each step)
+            if $V_dense == Val{true}
+                tv, xv, psol, tvS, xvS, gvS = sol
+                return TaylorInterpolant(jd0 - JD_J2000, tv .- tv[1], psol), tvS, xvS, gvS
+            # Point output
+            elseif $V_dense == Val{false}
+                return sol
+            end
 
         end
 

test/propagation.jl

@@ -281,6 +281,7 @@ using InteractiveUtils: methodswith
         rm("test.jld2")
 
         sol, tvS, xvS, gvS = NEOs.propagate_root(dynamics, 1, jd0, 0.02, q0, Val(true); order, abstol, parse_eqs)
+        @test sol isa TaylorInterpolant
 
         jldsave("test.jld2"; sol, tvS, xvS, gvS)
         recovered_sol = JLD2.load("test.jld2", "sol")