diff --git a/src/orbit_determination/gauss_method.jl b/src/orbit_determination/gauss_method.jl index 429dd480..b296f8b1 100755 --- a/src/orbit_determination/gauss_method.jl +++ b/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,7 +485,7 @@ 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 @@ -493,7 +493,7 @@ function gaussinitcond(radec::Vector{RadecMPC{T}}, tracklets::Vector{Tracklet{T} 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 diff --git a/src/orbit_determination/outlier_rejection.jl b/src/orbit_determination/outlier_rejection.jl index d1ab00eb..baf94a0a 100644 --- a/src/orbit_determination/outlier_rejection.jl +++ b/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 diff --git a/src/orbit_determination/tooshortarc.jl b/src/orbit_determination/tooshortarc.jl index baf74e21..aa371735 100644 --- a/src/orbit_determination/tooshortarc.jl +++ b/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 diff --git a/src/propagation/parameters.jl b/src/propagation/parameters.jl index 5b69756c..a1503759 100644 --- a/src/propagation/parameters.jl +++ b/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} diff --git a/test/orbit_determination.jl b/test/orbit_determination.jl index 29bab1f3..ec0e3ceb 100755 --- a/test/orbit_determination.jl +++ b/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 \ No newline at end of file