Update OD tests

PerezHz · Dec 12, 2024 · 0f57a5f · 0f57a5f
1 parent ed13e57
commit 0f57a5f
Showing 1 changed file with 73 additions and 17 deletions.
diff --git a/test/orbitdetermination.jl b/test/orbitdetermination.jl
@@ -7,7 +7,7 @@ using LinearAlgebra
 using Test
 
 using NEOs: NEOSolution, RadecMPC, reduce_tracklets,
-    indices, numberofdays, nout
+    indices, numberofdays, nout, scaled_variables
 
 function iodsubradec(radec::Vector{RadecMPC{T}}, N::Int = 3) where {T <: Real}
     tracklets = reduce_tracklets(radec)
@@ -17,7 +17,7 @@ function iodsubradec(radec::Vector{RadecMPC{T}}, N::Int = 3) where {T <: Real}
 end
 
 @testset "Orbit Determination" begin
-    @testset "Gauss Method (without ADAM)" begin
+    @testset "Straight Gauss Method" begin
         # Load observations
         radec = read_radec_mpc(joinpath(pkgdir(NEOs), "test", "data", "RADEC_2023_DW.dat"))
         # Subset of radec for IOD
@@ -112,30 +112,30 @@ end
         @test uncertaintyparameter(od, sol1, params) == 7
     end
 
-    @testset "Gauss Method (with ADAM)" begin
+    @testset "Unsafe Gauss Method" begin
         # Load observations
-        radec = fetch_radec_mpc("2016 TU93")
+        radec = fetch_radec_mpc("2005 TM173")
         # Subset of radec for IOD
         # subradec = iodsubradec(radec, 3)
         # In this case: radec == subradec
 
-        @test length(radec) == 9
-        @test numberofdays(radec) < 13.1
+        @test length(radec) == 6
+        @test numberofdays(radec) < 1.95
 
         # Parameters
-        params = NEOParameters(bwdoffset = 0.007, fwdoffset = 0.007)
+        params = NEOParameters(bwdoffset = 0.007, fwdoffset = 0.007, safegauss = false)
         # Orbit determination problem
         od = ODProblem(newtonian!, radec)
 
         # Initial Orbit Determination
         sol = orbitdetermination(od, params)
 
-        # Values by Nov 21, 2024
+        # Values by Dec 11, 2024
 
         # Orbit solution
         @test isa(sol, NEOSolution{Float64, Float64})
         # Tracklets
-        @test length(sol.tracklets) == 3
+        @test length(sol.tracklets) == 2
         @test sol.tracklets[1].radec[1] == radec[1]
         @test sol.tracklets[end].radec[end] == radec[end]
         @test issorted(sol.tracklets)
@@ -148,23 +148,79 @@ end
         @test all( norm.(sol.fwd.x, Inf) .< 2 )
         @test isempty(sol.t_fwd) && isempty(sol.x_fwd) && isempty(sol.g_fwd)
         # Vector of residuals
-        @test length(sol.res) == 9
+        @test length(sol.res) == 6
         @test iszero(nout(sol.res))
         # Least squares fit
         @test sol.fit.success
-        @test all( sigmas(sol) .< 8e-5 )
-        @test all( snr(sol) .> 36)
+        @test all( sigmas(sol) .< 5e-3 )
+        @test all( snr(sol) .> 21.5)
         @test nrms(sol) < 0.46
         # Jacobian
         @test size(sol.jacobian) == (6, 6)
         @test !isdiag(sol.jacobian)
-        @test maximum(sol.jacobian) < 4.1e-5
+        @test maximum(sol.jacobian) < 6.1e-4
+        # Compatibility with JPL
+        JPL = [1.0042569058151192, 0.2231639040146286, 0.11513854178693468,
+            -0.010824212819531798, 0.017428798232689943, 0.0071046780555307385]
+        @test all(abs.(sol() - JPL) ./ sigmas(sol) .< 8.1e-3)
+        # MPC Uncertainty Parameter
+        @test uncertaintyparameter(od, sol, params) == 10
+    end
+
+    @testset "Gauss Method with ADAM refinement" begin
+        # Load observations
+        radec = fetch_radec_mpc("2024 MK")
+        # Subset of radec for IOD
+        subradec = radec[10:21]
+
+        @test length(subradec) == 12
+        @test numberofdays(subradec) < 42.8
+
+        # Parameters
+        params = NEOParameters(bwdoffset = 0.007, fwdoffset = 0.007)
+        # Orbit determination problem
+        od = ODProblem(newtonian!, subradec)
+
+        # Initial Orbit Determination
+        varorder = max(params.tsaorder, params.gaussorder, params.jtlsorder)
+        scaled_variables("dx", ones(6); order = varorder)
+        sol = gaussiod(od, params)
+
+        # Values by Dec 11, 2024
+
+        # Orbit solution
+        @test isa(sol, NEOSolution{Float64, Float64})
+        # Tracklets
+        @test length(sol.tracklets) == 3
+        @test sol.tracklets[1].radec[1] == subradec[1]
+        @test sol.tracklets[end].radec[end] == subradec[end]
+        @test issorted(sol.tracklets)
+        # Backward integration
+        @test dtutc2days(date(subradec[1])) > sol.bwd.t0 + sol.bwd.t[end]
+        @test all( norm.(sol.bwd.x, Inf) .< 2 )
+        @test isempty(sol.t_bwd) && isempty(sol.x_bwd) && isempty(sol.g_bwd)
+        # Forward integration
+        @test dtutc2days(date(subradec[end])) < sol.fwd.t0 + sol.fwd.t[end]
+        @test all( norm.(sol.fwd.x, Inf) .< 2 )
+        @test isempty(sol.t_fwd) && isempty(sol.x_fwd) && isempty(sol.g_fwd)
+        # Vector of residuals
+        @test length(sol.res) == 12
+        @test iszero(nout(sol.res))
+        # Least squares fit
+        @test sol.fit.success
+        @test all( sigmas(sol) .< 6.6e-4 )
+        @test all( snr(sol) .> 38.8)
+        @test nrms(sol) < 0.32
+        # Jacobian
+        @test size(sol.jacobian) == (6, 6)
+        @test isdiag(sol.jacobian)
+        @test maximum(sol.jacobian) < 9.5e-6
         # Compatibility with JPL
-        JPL = [1.0102558767253402, 0.2935121552882981, 0.10468669797982912,
-            -0.0002639687633186843, 0.01837366168395344, 0.007208431369660604]
-        @test all(abs.(sol() - JPL) ./ sigmas(sol) .< 4.8e-2)
+        JPL = [-0.12722461679828806, -0.9466098076903212, -0.4526816007640767,
+            0.02048875631534963, -0.00022720097573790754, 0.00321302850930331]
+        @test all(abs.(sol() - JPL) ./ sigmas(sol) .< 0.16)
         # MPC Uncertainty Parameter
-        @test uncertaintyparameter(od, sol, params) == 8
+        @test uncertaintyparameter(od, sol, params) == 9
     end
 
     @testset "Admissible region" begin