New methods of IOD (#85)

* Add adamQtol to NEOParameters

* Introduce addradec! to jtls

* Bimodal addradec!

* Introduce iod(::Vector{RadecMPC}, ::NEOParameters)

* Update scripts/names.jl

* Update scripts/orbitdetermination.jl

* Iterate tracklets in iod

* Use on_error in scripts/orbitdetermination.jl

* Add initcond as iod kwarg

* Update iod kwargs signature

* Fix iod

* Introduce updatesol to iod

* Manual weights for F52, U68 and L01

* Fill rin to have >= 3 obs in jtls

* ADAM requires a minimum of 2 observations

* Introduce _gaussiod and _tsaiod

* Fix _gausstriplets2!

* Update params in scripts/orbitdetermination.jl

* Update scripts/names.jl

* Add discovery filter to scripts/

* Choose best orbit in scripts/orbitdetermination.jl

* Introduce isodcompatible to scripts/names.jl

* Fix filter in scripts/

* Use MPC OBS API to fetch optical astrometry

* Add date > 2000 filter to scripts/

* Introduce _initialtracklets to jtls

* Fix fetch_radec_mpc signature

* Update test/orbitdetermination.jl

* Add iod tests

* Fix iod test

* Bump patch version

* Deprecate gaussinitcond and tooshortarc

* Use dtutc2jdtdb instead of datetime2julian in OD

* Change datetime for dtutc in remaining units.jl functions

* Fix tests

* Fix NEOCP test

* Use MPC NEOCP OBS API to fetch NEOCP astrometry

Project.toml

@@ -1,7 +1,7 @@
 name = "NEOs"
 uuid = "b41c07a2-2abb-11e9-070a-c3c1b239e7df"
 authors = ["Jorge A. Pérez Hernández", "Luis Benet", "Luis Eduardo Ramírez Montoya"]
-version = "0.9.0"
+version = "0.9.1"
 
 [deps]
 Artifacts = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

examples/recipes.jl

@@ -3,7 +3,7 @@
 using NEOs, Plots
 
 # Download optical astrometry of asteroid 2023 DW
-radec = fetch_radec_mpc("designation" => "2023 DW")
+radec = fetch_radec_mpc("2023 DW")
 
 # Orbit determination
 params = NEOParameters(bwdoffset = 0.007, fwdoffset = 0.007)

pha/apophis.jl

@@ -41,7 +41,7 @@ function parse_commandline()
 
     @add_arg_table! s begin
         "--jd0"
-            help = "Initial epoch calendar date (TDB)"
+            help = "Initial epoch calendar date (UTC)"
             arg_type = DateTime
             default = DateTime(2020, 12, 17)
         "--varorder"
@@ -121,7 +121,7 @@ function main(dynamics::D, maxsteps::Int, jd0_datetime::DateTime, nyears_bwd::T,
     q0 = vcat(q00, 0.0, 0.0) .+ dq
 
     # Initial date (in Julian days)
-    jd0 = datetime2julian(jd0_datetime)
+    jd0 = dtutc2jdtdb(jd0_datetime)
 
     print_header("Integrator warmup", 2)
     params = NEOParameters(;maxsteps=1, order, abstol, parse_eqs)
@@ -132,7 +132,7 @@ function main(dynamics::D, maxsteps::Int, jd0_datetime::DateTime, nyears_bwd::T,
     print_header("Main integration", 2)
     tmax = nyears_bwd*yr
     println("• Initial time of integration: ", string(jd0_datetime))
-    println("• Final time of integration: ", julian2datetime(jd0 + tmax))
+    println("• Final time of integration: ", jdtdb2utc(jd0 + tmax))
 
     params = NEOParameters(;maxsteps, order, abstol, parse_eqs)
     sol_bwd = NEOs.propagate(dynamics, jd0, nyears_bwd, q0, params)
@@ -141,7 +141,7 @@ function main(dynamics::D, maxsteps::Int, jd0_datetime::DateTime, nyears_bwd::T,
 
     tmax = nyears_fwd*yr
     println("• Initial time of integration: ", string(jd0_datetime))
-    println("• Final time of integration: ", julian2datetime(jd0 + tmax))
+    println("• Final time of integration: ", jdtdb2utc(jd0 + tmax))
 
     sol_fwd, tvS, xvS, gvS = NEOs.propagate_root(dynamics, jd0, nyears_fwd, q0, params)
     jldsave("Apophis_fwd.jld2"; sol_fwd, tvS, xvS, gvS)

scripts/adam.jl

@@ -50,7 +50,7 @@ end
         μ::T = 0.75, ν::T = 0.9, ϵ::T = 1e-8, Qtol::T = 0.001, adamorder::Int = 2,
         dynamics::D = newtonian!) where {T <: Real, D}
         # Initial time of integration [julian days]
-        jd0 = datetime2julian(A.date)
+        jd0 = dtutc2jdtdb(A.date)
         # Maximum number of iterations
         maxiter = params.adamiter
         # Allocate memory
@@ -71,7 +71,7 @@ end
         set_variables(Float64, "dx"; order = adamorder, numvars = 6)
         dae = [scalings[i] * TaylorN(i, order = adamorder) for i in 1:6]
         # Propagation buffer
-        t0, tf = datetime2days(date(radec[1])), datetime2days(date(radec[end]))
+        t0, tf = dtutc2days(date(radec[1])), dtutc2days(date(radec[end]))
         tlim = (t0 - params.bwdoffset, tf + params.fwdoffset)
         buffer = PropagationBuffer(dynamics, jd0, tlim, aes[:, 1] .+ dae, params)
         # Vector of O-C residuals

scripts/mcmov.jl

@@ -61,9 +61,9 @@ end
         _, y_max = argoldensearch(A, 10^x_min, 10^x_max, :max, :outer, 1e-20)
         scalings[5:6] .= (x_max - x_min) / 100, (y_max - y_min) / 100
         # Propagation buffer
-        t0, tf = datetime2days(date(radec[1])), datetime2days(date(radec[end]))
+        t0, tf = dtutc2days(date(radec[1])), dtutc2days(date(radec[end]))
         tlim = (t0 - params.bwdoffset, tf + params.fwdoffset)
-        buffer = PropagationBuffer(newtonian!, datetime2julian(A.date), tlim,
+        buffer = PropagationBuffer(newtonian!, dtutc2jdtdb(A.date), tlim,
                                    ae .+ scalings .* dae, params)
         # Vector of residuals
         res = Vector{OpticalResidual{T, TaylorN{T}}}(undef, length(radec))
@@ -82,7 +82,7 @@ end
             # Barycentric initial conditions (JT)
             q = attr2bary(A, AE, params)
             # Reference epoch in julian days (corrrected for light-time)
-            jd0 = datetime2julian(A.date) - constant_term(AE[5]) / c_au_per_day
+            jd0 = dtutc2jdtdb(A.date) - constant_term(AE[5]) / c_au_per_day
             # Propagation and residuals
             propres!(res, radec, jd0, q, params; buffer)
             # Least squares fit

scripts/names.jl

@@ -1,5 +1,5 @@
-using ArgParse, Downloads, Random, HORIZONS, NEOs
-using NEOs: numberofdays
+using ArgParse, Downloads, Random, NEOs, Dates
+using NEOs: RadecMPC, numberofdays, issatellite, reduce_tracklets
 
 # Potentially Hazardous Asteroids MPC list
 const PHAS_URL = "https://cgi.minorplanetcenter.net/cgi-bin/textversion.cgi?f=lists/PHAs.html"
@@ -20,37 +20,61 @@ function parse_commandline()
 
     @add_arg_table! s begin
         "--N", "-N"
-            help = "Number of NEOs"
+            help = "number of NEOs"
             arg_type = Int
-            default = 100
+            default = 0
         "--output", "-o"
-            help = "Output file"
+            help = "output file"
             arg_type = String
             default = "names.txt"
     end
 
     s.epilog = """
-        examples:\n
-        \n
-        julia --project names.jl -N 250 -o names.txt\n
+        note: if izero(N), then the script will save all compatible NEOs\n
         \n
+        examples:\n
+        # Save 1000 NEOs to names.txt\n
+        julia --project names.jl -N 1000 -o names.txt\n
     """
 
     return parse_args(s)
 end
 
+function isodcompatible(radec::Vector{RadecMPC{T}}) where {T <: Real}
+    # Eliminate observations before oficial discovery
+    firstobs = findfirst(r -> !isempty(r.discovery), radec)
+    isnothing(firstobs) && return false
+    radec = radec[firstobs:end]
+    # Filter out incompatible observations
+    filter!(radec) do r
+        hascoord(r.observatory) && !issatellite(r.observatory) &&
+        date(r) >= Date(2000)
+    end
+    length(radec) < 3 && return false
+    # There is at least one set of 3 tracklets with a < 15 days timespan
+    tracklets = reduce_tracklets(radec)
+    for i in 1:length(tracklets)-2
+        numberofdays(tracklets[i:i+2]) > 15.0 && continue
+        tracklets = tracklets[i:i+2]
+        radec = reduce(vcat, getfield.(tracklets, :radec))
+        sort!(radec)
+        break
+    end
+    return numberofdays(radec) <= 15.0
+end
+
 function main()
 
     # Parse arguments from commandline
     parsed_args = parse_commandline()
-
-    println("NEOs selector for orbitdetermination.jl")
-
     # Number of NEOs
     N::Int = parsed_args["N"]
-    println("• Number of NEOs: ", N)
     # Output file
-    outfile::String = parsed_args["output"]
+    output::String = parsed_args["output"]
+    # Print header
+    println("NEOs selector for orbitdetermination.jl")
+    println("• Number of NEOs: ", N)
+    println("• Output file: ", output)
 
     # Parse PHAs, Atens, Apollos and Amors
     provdesig = Vector{SubString{String}}(undef, 0)
@@ -71,34 +95,40 @@ function main()
     unique!(provdesig)
     # We can only process asteroids discovered between 2000 and 2024
     filter!(desig -> "2000" <= desig[1:4] <= "2024", provdesig)
-    # Shuffle the provisional designations list
-    shuffle!(provdesig)
+
+    if 0 < N < length(provdesig)
+        # Shuffle the provisional designations list
+        shuffle!(provdesig)
+    else
+        # Consider all compatible NEOs
+        N = length(provdesig)
+    end
 
     # Assemble a list with N NEOs
     names = Vector{String}(undef, N)
 
-    for i in eachindex(names)
-        while !isempty(provdesig)
-            neo = String(pop!(provdesig))
-            radec = fetch_radec_mpc("designation" => neo)
-            jplorbit = sbdb("des" => neo)["orbit"]
-            if (numberofdays(radec) <= 30.0) &&  (jplorbit["n_obs_used"] == length(radec)) &&
-                isnothing(jplorbit["n_del_obs_used"]) && isnothing(jplorbit["n_dop_obs_used"])
-                names[i] = neo
-                break
-            end
+    i = 1
+    while i <= N && !isempty(provdesig)
+        neo = String(popfirst!(provdesig))
+        radec = fetch_radec_mpc(neo)
+        if isodcompatible(radec)
+            names[i] = neo
+            i += 1
         end
     end
+    # Eliminate #undef elements
+    names = names[1:i-1]
+    # Sort provisional designations
     sort!(names)
 
     # Save names list
-    open(outfile, "w") do file
+    open(output, "w") do file
         for i in eachindex(names)
             write(file, names[i], i == N ? "" : "\n")
         end
     end
 
-    println("• Saved ", N, " names to ", outfile)
+    println("• Saved ", N, " names to ", output)
 
     nothing
 end