diff --git a/test/orbitdetermination.jl b/test/orbitdetermination.jl index fa1d12b0..bae4f97c 100644 --- 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