Add 2008 TC3 test

src/orbit_determination/gauss_method.jl

@@ -467,9 +467,9 @@ function gaussinitcond(radec::Vector{RadecMPC{T}}, tracklets::Vector{Tracklet{T}
         # Julian day when to start propagation
         jd0 = datetime2julian(dates[triplet[2]])
         # Number of years in forward integration
-        nyears_fwd = (tf - jd0 + 2) / yr
+        nyears_fwd = (tf - jd0 + params.fwdoffset) / yr
         # Number of years in backward integration
-        nyears_bwd = -(jd0 - t0 + 2) / yr
+        nyears_bwd = -(jd0 - t0 + params.bwdoffset) / yr
         # Gauss method solution 
         sol = gauss_method(observatories[triplet], dates[triplet], α[triplet] .+ dq[1:3],
                            δ[triplet] .+ dq[4:6]; niter = params.niter)
@@ -485,15 +485,15 @@ function gaussinitcond(radec::Vector{RadecMPC{T}}, tracklets::Vector{Tracklet{T}
             bwd::TaylorInterpolant{T, TaylorN{T}, 2} = propagate(
                 RNp1BP_pN_A_J23E_J2S_eph_threads!, jd0, nyears_bwd, q0, params
             ) 
-            if !issuccessfulprop(bwd, t0 - jd0)
+            if !issuccessfulprop(bwd, t0 - jd0; tol = params.coeffstol)
                 continue
             end
 
             # Forward propagation
             fwd::TaylorInterpolant{T, TaylorN{T}, 2} = propagate(
                 RNp1BP_pN_A_J23E_J2S_eph_threads!, jd0, nyears_fwd, q0, params
             )
-            if !issuccessfulprop(fwd, tf - jd0)
+            if !issuccessfulprop(fwd, tf - jd0; tol = params.coeffstol)
                 continue
             end
 

src/orbit_determination/outlier_rejection.jl

@@ -35,9 +35,9 @@ function outlier_rejection(radec::Vector{RadecMPC{T}}, sol::NEOSolution{T, T},
     # Julian day of first (last) observation
     t0, tf = datetime2julian(date(radec[1])), datetime2julian(date(radec[end]))
     # Number of years in forward integration 
-    nyears_fwd = (tf - jd0 + 2) / yr
+    nyears_fwd = (tf - jd0 + params.fwdoffset) / yr
     # Number of years in backward integration
-    nyears_bwd = -(jd0 - t0 + 2) / yr
+    nyears_bwd = -(jd0 - t0 + params.bwdoffset) / yr
     # Dynamical function 
     dynamics = RNp1BP_pN_A_J23E_J2S_eph_threads!
     # Initial conditions (T)
@@ -52,14 +52,14 @@ function outlier_rejection(radec::Vector{RadecMPC{T}}, sol::NEOSolution{T, T},
     # Backward integration
     bwd = propagate(dynamics, jd0, nyears_bwd, q, params)
 
-    if !issuccessfulprop(bwd, t0 - jd0)
+    if !issuccessfulprop(bwd, t0 - jd0; tol = params.coeffstol)
         return zero(NEOSolution{T, T})
     end
 
     # Forward integration
     fwd = propagate(dynamics, jd0, nyears_fwd, q, params)
 
-    if !issuccessfulprop(fwd, tf - jd0)
+    if !issuccessfulprop(fwd, tf - jd0; tol = params.coeffstol)
         return zero(NEOSolution{T, T})
     end
 

src/orbit_determination/tooshortarc.jl

@@ -321,12 +321,13 @@ function propres(radec::Vector{RadecMPC{T}}, jd0::T, q0::Vector{U},
     # Time of first (last) observation
     t0, tf = datetime2julian(date(radec[1])), datetime2julian(date(radec[end]))
     # Years in backward (forward) integration
-    nyears_bwd = -(jd0 - t0 + 0.5) / yr
-    nyears_fwd = (tf - jd0 + 0.5) / yr
+    nyears_bwd = -(jd0 - t0 + params.bwdoffset) / yr
+    nyears_fwd = (tf - jd0 + params.fwdoffset) / yr
     # Backward (forward) integration
     bwd = propagate(RNp1BP_pN_A_J23E_J2S_eph_threads!, jd0, nyears_bwd, q0, params)
     fwd = propagate(RNp1BP_pN_A_J23E_J2S_eph_threads!, jd0, nyears_fwd, q0, params)
-    if !issuccessfulprop(bwd, t0 - jd0) || !issuccessfulprop(fwd, tf - jd0)
+    if !issuccessfulprop(bwd, t0 - jd0; tol = params.coeffstol) || 
+       !issuccessfulprop(fwd, tf - jd0; tol = params.coeffstol)
         return bwd, fwd, Vector{OpticalResidual{T, U}}(undef, 0)
     end
     # Sun (Earth) ephemeris

src/propagation/parameters.jl

@@ -5,6 +5,9 @@ struct NEOParameters{T <: AbstractFloat}
     order::Int
     abstol::T
     parse_eqs::Bool
+    bwdoffset::T
+    fwdoffset::T
+    coeffstol::T
     # Residuals parameters
     mpc_catalogue_codes_201X::Vector{String}
     truth::String
@@ -37,6 +40,8 @@ Parameters for all orbit determination functions.
 - `order::Int`: order of Taylor expansions wrt time (default: 25).
 - `abstol::T`: absolute tolerance used to compute propagation timestep (default: `1e-20`).
 - `parse_eqs::Bool`: whether to use the specialized method of `jetcoeffs` or not (default: `true`).
+- `bwdoffset/fwdoffset::T`: days to propagate beyond first (bwd) / last (fwd) observation (default: `0.5`).
+- `coeffstol::T`: maximum size of the coefficients (default: `10.0`).
 
 # Residuals Parameters
 
@@ -64,16 +69,17 @@ Parameters for all orbit determination functions.
 - `max_per::T`: maximum allowed rejection percentage (default: `18.0`).
 """
 function NEOParameters(; maxsteps::Int = 500, μ_ast::Vector{T} = μ_ast343_DE430[1:end],
-                      order::Int = 25, abstol::T = 1e-20, parse_eqs::Bool = true, 
-                      debias_table::String = "2018", niter::Int = 5, max_triplets::Int = 10,
-                      varorder::Int = 5, Q_max::T = 5.0, maxiter::Int = 100,
-                      max_per::T = 18.0) where {T <: AbstractFloat}
+                      order::Int = 25, abstol::T = 1e-20, parse_eqs::Bool = true, bwdoffset::T = 0.5,
+                      fwdoffset::T = 0.5, coeffstol::T = 10.0, debias_table::String = "2018",
+                      niter::Int = 5, max_triplets::Int = 10, varorder::Int = 5, Q_max::T = 5.0,
+                      maxiter::Int = 100, max_per::T = 18.0) where {T <: AbstractFloat}
     # Unfold debiasing matrix
     mpc_catalogue_codes_201X, truth, resol, bias_matrix = select_debiasing_table(debias_table)
     # Assemble NEOParameters
-    return NEOParameters{T}(maxsteps, μ_ast, order, abstol, parse_eqs, mpc_catalogue_codes_201X,
-                         truth, resol, bias_matrix, niter, max_triplets, varorder, Q_max,
-                         maxiter, max_per)
+    return NEOParameters{T}(maxsteps, μ_ast, order, abstol, parse_eqs, bwdoffset,
+                            fwdoffset, coeffstol, mpc_catalogue_codes_201X, truth,
+                            resol, bias_matrix, niter, max_triplets, varorder, Q_max,
+                            maxiter, max_per)
 end                            
 
 function NEOParameters(params::NEOParameters{T}; kwargs...) where {T <: AbstractFloat}

test/orbit_determination.jl

@@ -1,6 +1,7 @@
 # This file is part of the NEOs.jl package; MIT licensed
 
 using NEOs
+using Dates
 using PlanetaryEphemeris
 using LinearAlgebra
 using Test
@@ -141,11 +142,12 @@ using NEOs: NEOSolution, numberofdays
     end
 
     @testset "Interesting NEOs" begin
+        using NEOs: reduce_tracklets, selecteph, sseph, propres, evalfit
 
-        # 2014 AA hit the Earth around January 2, 2014 02:49 UTC
+        # 2014 AA hit the Earth around January 2, 2014, 02:49 UTC
 
         # Optical astrometry file
-        filename = joinpath("data", "2011_AA.txt")
+        filename = joinpath("data", "2014_AA.txt")
         # Download observations
         get_radec_mpc("designation" => "2014 AA", filename)
         # Load observations 
@@ -187,6 +189,104 @@ using NEOs: NEOSolution, numberofdays
         @test nrms(sol) < 0.13
         # Scalig factors
         @test all(sol.scalings .< 1e-5)
+
+        # 2008 TC3 entered the Earth's atmosphere around October 7, 2008, 02:46 UTC
+
+        # Optical astrometry file
+        filename = joinpath("data", "2008_TC3.txt")
+        # Download observations
+        get_radec_mpc("designation" => "2008 TC3", filename)
+        # Load observations 
+        radec = read_radec_mpc(filename)
+        # Delete astrometry file
+        rm(filename)
+        # Parameters
+        params = NEOParameters(abstol = 1e-20, order = 25, parse_eqs = true,
+                               coeffstol = Inf)
+
+        # Observations before  October 6, 2008, 10:00 UTC
+        idxs = findall(x -> date(x) < DateTime(2008, 10, 06, 10), radec)
+        # Restricted Orbit Determination
+        sol = orbitdetermination(radec[idxs], params)
+
+        # Tracklets 
+        tracklets = reduce_tracklets(radec)
+        # Time of first (last) observation [julian days]
+        t0, tf = datetime2julian(date(radec[1])), datetime2julian(date(radec[end]))
+        # Sun (earth's) ephemeris
+        eph_su = selecteph(NEOs.sseph, su)
+        eph_ea = selecteph(sseph, ea)
+
+        # Initial time of propagation [julian days]
+        jd0 = sol.bwd.t0 + J2000
+        # Years in backward (forward) propagation
+        nyears_bwd = -(jd0 - t0 + 0.5) / yr
+        nyears_fwd = (tf - jd0 + 0.5) / yr
+        # Initial conditions
+        q0 = sol()
+        scalings = abs.(q0) ./ 10^6
+        dq = NEOs.scaled_variables("dx", scalings, order = 5)
+        q = q0 + dq
+        # Propagation and residuals
+        bwd = NEOs.propagate(RNp1BP_pN_A_J23E_J2S_eph_threads!, jd0, nyears_bwd, q, params)
+        fwd = NEOs.propagate(RNp1BP_pN_A_J23E_J2S_eph_threads!, jd0, nyears_fwd, q, params)
+        res = residuals(radec, params;
+                        xvs = et -> auday2kmsec(eph_su(et/daysec)),
+                        xve = et -> auday2kmsec(eph_ea(et/daysec)),
+                        xva = et -> bwdfwdeph(et, bwd, fwd))
+        # Least squares fit
+        fit = tryls(res, zeros(6), params.niter)
+        # Orbit with all the observations
+        sol = evalfit(NEOSolution(tracklets, bwd, fwd, res, fit, scalings))
+
+        # Initial time of propagation [julian days]
+        jd0 = sol.bwd.t0 + J2000
+        # Initial conditions
+        q0 = sol()
+        scalings = abs.(q0) ./ 10^6
+        dq = NEOs.scaled_variables("dx", scalings, order = 5)
+        q = q0 + dq
+        # Propagation and residuals
+        bwd = NEOs.propagate(RNp1BP_pN_A_J23E_J2S_eph_threads!, jd0, nyears_bwd, q, params)
+        fwd = NEOs.propagate(RNp1BP_pN_A_J23E_J2S_eph_threads!, jd0, nyears_fwd, q, params)
+        res = residuals(radec, params;
+                        xvs = et -> auday2kmsec(eph_su(et/daysec)),
+                        xve = et -> auday2kmsec(eph_ea(et/daysec)),
+                        xva = et -> bwdfwdeph(et, bwd, fwd))
+        # Least squares fit
+        fit = tryls(res, zeros(6), params.niter)
+        # Orbit refinement
+        sol = evalfit(NEOSolution(tracklets, bwd, fwd, res, fit, scalings))
+
+        # Values by January 26, 2024
+
+        # Vector of observations
+        @test length(radec) == 883
+        @test numberofdays(radec) < 0.80
+        # Orbit solution
+        @test isa(sol, NEOSolution{Float64, Float64})
+        # Tracklets
+        @test length(sol.tracklets) == 29
+        @test sol.tracklets[1].radec[1] == radec[1]
+        @test sol.tracklets[end].radec[end] == radec[end]
+        @test issorted(sol.tracklets)
+        # 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) .< 1e55 )
+        @test isempty(sol.t_fwd) && isempty(sol.x_fwd) && isempty(sol.g_fwd)
+        # Vector of residuals
+        @test length(sol.res) == 883
+        @test iszero(count(outlier.(sol.res)))
+        # Least squares fit
+        @test sol.fit.success
+        @test all( sigmas(sol) .< 1e-7 )
+        @test nrms(sol) < 0.61
+        # Scalig factors
+        @test all(sol.scalings .< 1e-5)
     end
 
 end