Skip to content

Commit

Permalink
Update TSA tests
Browse files Browse the repository at this point in the history
  • Loading branch information
@LuEdRaMo
LuEdRaMo committed Dec 22, 2023
1 parent fe9d054 commit 7237c66
Showing 1 changed file with 100 additions and 48 deletions.
148 changes: 100 additions & 48 deletions test/orbit_determination.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,59 +6,111 @@ using LinearAlgebra
using Test

using NEOs: NEOSolution, numberofdays
@testset "Gauss Method" begin
# Load observations
radec = read_radec_mpc(joinpath("data", "RADEC_2023_DW.dat"))
# Parameters
params = Parameters(abstol = 1e-20, order = 25, parse_eqs = true)

# Orbit Determination
sol = orbitdetermination(radec, params)
@testset "Orbit Determination" begin
@testset "Gauss Method" begin
# Load observations
radec = read_radec_mpc(joinpath("data", "RADEC_2023_DW.dat"))
# Parameters
params = Parameters(abstol = 1e-20, order = 25, parse_eqs = true)

# Orbit Determination
sol = orbitdetermination(radec, params)

# Values by December 21, 2023

# Vector of observations
@test length(radec) == 123
@test numberofdays(radec) < 21.0
# Orbit solution
@test isa(sol, NEOSolution{Float64, Float64})
# Observation nights
@test length(sol.nights) == 44
@test sol.nights[1].radec[1] == radec[1]
@test sol.nights[end].radec[end] == radec[end]
@test issorted(sol.nights)
# Backward integration
@test datetime2days(date(radec[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 datetime2days(date(radec[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) == 123
@test iszero(count(outlier.(sol.res)))
# Least squares fit
@test sol.fit.success
@test all( sigmas(sol) .< 1e-4 )
@test nrms(sol) < 0.38
# Scalig factors
@test all(sol.scalings .== 1e-6)
end

# Values by December 21, 2023
@testset "Too Short Arc" begin
# Optical astrometry file
filename = joinpath("data", "2008_EK68.txt")
# Download observations
get_radec_mpc("designation" => "2008 EK68", filename)
# Load observations
radec = read_radec_mpc(filename)
# Delete astrometry file
rm(filename)
# Parameters
params = Parameters(abstol = 1e-20, order = 25, parse_eqs = true)

# Vector of observations
@test length(radec) == 123
@test numberofdays(radec) < 21.0
# Orbit solution
@test isa(sol, NEOSolution{Float64, Float64})
# Observation nights
@test length(sol.nights) == 44
@test sol.nights[1].radec[1] == radec[1]
@test sol.nights[end].radec[end] == radec[end]
# Backward integration
@test datetime2days(date(radec[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 datetime2days(date(radec[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) == 123
@test iszero(count(outlier.(sol.res)))
# Least squares fit
@test sol.fit.success
@test all( sigmas(sol) .< 1e-4 )
@test nrms(sol) < 0.38
# Scalig factors
@test all(sol.scalings .== 1e-6)
# Orbit Determination
sol = orbitdetermination(radec, params)

# Gauss refinement (with outlier rejection)
# Values by December 22, 2023

# Vector of observations
@test length(radec) == 10
@test numberofdays(radec) < 0.05
# Orbit solution
@test isa(sol, NEOSolution{Float64, Float64})
# Observation nights
@test length(sol.nights) == 1
@test sol.nights[1].radec[1] == radec[1]
@test sol.nights[end].radec[end] == radec[end]
@test issorted(sol.nights)
# Backward integration
@test datetime2days(date(radec[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 datetime2days(date(radec[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) == 10
@test iszero(count(outlier.(sol.res)))
# Least squares fit
@test sol.fit.success
@test all( sigmas(sol) .< 6e-3 )
@test nrms(sol) < 0.85
# Scalig factors
@test all(sol.scalings .< 1e-5)
end
#=
sol = gauss_refinement(radec, sol, params)
@test "Outlier Rejection" begin
# Gauss refinement (with outlier rejection)
@test isa(sol, NEOSolution{Float64, Float64})
@test datetime2days(date(radec[1])) > sol.bwd.t0 + sol.bwd.t[end]
@test datetime2days(date(radec[end])) < sol.fwd.t0 + sol.fwd.t[end]
@test all( norm.(sol.bwd.x, Inf) .< 2 )
@test all( norm.(sol.fwd.x, Inf) .< 2 )
@test isempty(sol.t_bwd) && isempty(sol.x_bwd) && isempty(sol.g_bwd)
@test isempty(sol.t_fwd) && isempty(sol.x_fwd) && isempty(sol.g_fwd)
@test length(sol.res) == 123
@test sol.fit.success
@test all( sqrt.(diag(sol.fit.Γ)) .< 100 )
@test nrms(sol) < 0.4
@test count(outlier.(sol.res)) == 0
sol = gauss_refinement(radec, sol, params)
@test isa(sol, NEOSolution{Float64, Float64})
@test datetime2days(date(radec[1])) > sol.bwd.t0 + sol.bwd.t[end]
@test datetime2days(date(radec[end])) < sol.fwd.t0 + sol.fwd.t[end]
@test all( norm.(sol.bwd.x, Inf) .< 2 )
@test all( norm.(sol.fwd.x, Inf) .< 2 )
@test isempty(sol.t_bwd) && isempty(sol.x_bwd) && isempty(sol.g_bwd)
@test isempty(sol.t_fwd) && isempty(sol.x_fwd) && isempty(sol.g_fwd)
@test length(sol.res) == 123
@test sol.fit.success
@test all( sqrt.(diag(sol.fit.Γ)) .< 100 )
@test nrms(sol) < 0.4
@test count(outlier.(sol.res)) == 0
end
=#
end

0 comments on commit 7237c66

Please sign in to comment.