From 081e19de8c1858ab80c19288d28aebe458a60903 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Luis=20Eduardo=20Ram=C3=ADrez=20Montoya?=
 <luiseduardoramirezmontoya@ciencias.unam.mx>
Date: Sun, 23 Jun 2024 14:07:26 -0600
Subject: [PATCH] Add bplane tests

---
 test/bplane.jl   | 53 ++++++++++++++++++++++++++++++++++++++++++++++++
 test/runtests.jl |  1 +
 2 files changed, 54 insertions(+)
 create mode 100644 test/bplane.jl

diff --git a/test/bplane.jl b/test/bplane.jl
new file mode 100644
index 00000000..3b53dbaf
--- /dev/null
+++ b/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
\ No newline at end of file
diff --git a/test/runtests.jl b/test/runtests.jl
index 509eb0f3..deced224 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -5,6 +5,7 @@ testfiles = (
     "observations.jl",
     "propagation.jl",
     "orbitdetermination.jl",
+    "bplane.jl",
     "dataframes.jl",
     "aqua.jl",
     )