Modified Target Plane (#70)

* Add mtp(::Vector{U})

* Add bplane tests

* Bump patch version

* Add Roots to test/Project.toml

Author: LuEdRaMo

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

src/NEOs.jl

@@ -75,7 +75,7 @@ export gauss_method, gaussinitcond
 # Orbit determination
 export curvature, issinglearc, isgauss, orbitdetermination
 # B plane
-export valsecchi_circle, bopik
+export valsecchi_circle, bopik, mtp
 # Outlier rejection
 export outlier_rejection
 

src/postprocessing/b_plane.jl src/postprocessing/bplane.jlsrc/postprocessing/b_plane.jl renamed to src/postprocessing/bplane.jl

@@ -206,4 +206,35 @@ function valsecchi_circle(U_y, U_norm, k, h; m_pl=3.003489614915764e-6, a_pl=1.0
     D0 = c*sinθ/(cosθ0p-cosθ)
 
     return R0, D0
+end
+
+@doc raw"""
+    mtp(xae::Vector{U}) where {U <: Number}
+
+Return the cartesian coordinates [Earth radii] on the modified
+target plane of an asteroid's close approach with the Earth.
+
+# Arguments
+
+- `xae::Vector{U}`: asteroid's geocentric position/velocity
+   vector at close approach [au, au/day].
+
+!!! reference
+    See equations (43)-(44) in page 15 of https://doi.org/10.1007/s10569-019-9914-4.
+"""
+function mtp(xae::Vector{U}) where {U <: Number}
+    # Unfold geocentric position and velocity
+    r, v = xae[1:3], xae[4:6]
+    # Reference direction
+    V = [0, 0, -1]
+    # Unit vectors
+    ez = v ./ euclid3D(v)
+    _ey_ = cross(V, ez)
+    ey = _ey_ ./ euclid3D(_ey_)
+    ex = cross(ey, ez)
+    # Cartesian coordinates [Earth radii]
+    X = dot3D(r, ex)*au/RE
+    Y = dot3D(r, ey)*au/RE
+
+    return X, Y
 end

src/postprocessing/outlier_rejection.jl

@@ -1,4 +1,4 @@
-include("b_plane.jl")
+include("bplane.jl")
 
 @doc raw"""
     residual_norm(x::OpticalResidual{T, T}) where {T <: Real}

test/Project.toml

@@ -11,6 +11,7 @@ Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 PlanetaryEphemeris = "d83715d0-7e5f-11e9-1a59-4137b20d8363"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Query = "1a8c2f83-1ff3-5112-b086-8aa67b057ba1"
+Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 TaylorIntegration = "92b13dbe-c966-51a2-8445-caca9f8a7d42"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

test/bplane.jl

@@ -0,0 +1,53 @@
+# This file is part of the NEOs.jl package; MIT licensed
+
+using NEOs
+using PlanetaryEphemeris
+using LinearAlgebra
+using Roots
+using Test
+
+using NEOs: propres
+
+@testset "2018 LA" begin
+
+    # Fetch optical astrometry
+    radec = fetch_radec_mpc("designation" => "2018 LA")
+    # Parameters
+    params = NEOParameters(coeffstol = Inf, bwdoffset = 0.007, fwdoffset = 0.007)
+
+    # Orbit Determination
+    sol = orbitdetermination(radec[1:8], params)
+    params = NEOParameters(params; jtlsorder = 6)
+    sol = orbitdetermination(radec, sol, params)
+
+    # Radial velocity with respect to the Earth.
+    function rvelea(t, fwd, params)
+        # Geocentric state vector
+        rv = fwd(t) - params.eph_ea(t)
+        # Derivative of geocentric distance
+        return dot(rv[1:3], rv[4:6])
+    end
+
+    # Values by June 23, 2024
+
+    # Time of close approach
+    params = NEOParameters(params; fwdoffset = 0.3)
+    bwd, fwd, res = propres(radec, epoch(sol) + J2000, sol(), params)
+    t_CA = find_zeros(t -> rvelea(t, fwd, params), fwd.t0, fwd.t0 + fwd.t[end-1])[1]
+    # Asteroid's geocentric state vector
+    xae = fwd(t_CA) - params.eph_ea(t_CA)
+    # Earth's heliocentric state vector
+    xes = params.eph_ea(t_CA) - params.eph_su(t_CA)
+
+    # Öpik's coordinates
+    B = bopik(xae, xes)
+    # Modified Target Plane
+    X, Y = mtp(xae)
+
+    # Impact parameter
+    @test B.b >= 1.0
+    # Impact condition
+    @test hypot(B.ξ, B.ζ) <= B.b
+    @test hypot(X, Y) <= 1
+
+end

test/runtests.jl

@@ -5,6 +5,7 @@ testfiles = (
     "observations.jl",
     "propagation.jl",
     "orbitdetermination.jl",
+    "bplane.jl",
     "dataframes.jl",
     "aqua.jl",
     )