Initial Orbit Determination (#31)

* Add QueryExt.jl

* Minor fix in QueryExt.jl

* Gauss initial conditions

* Remove show

* Change Query for DataFrames approach

* Minor fix in gauss_idxs

* Change to Interpolations.jl

* Modify gauss_idxs

* Document new functions

* Change gauss_idxs for gauss_triplets

* Fix gauss_triplets

* Add Q_max kwarg to gaussinitcond

* Minor fix

* Fix documentation typo

* Fix propagate_root return

* Heliocentric energy

* Modify gauss_triplets

* Change min delta

* Minor fix

* Add tryls

* Adaptative time delta

* Minor correction

* Modify gauss_triplets

* Relaxation factor

* Fix infinite while

* Fix M_

* Change reduce from Date to TimeOfDay

* Add preliminary docstrings

* Try fit Gauss to full radec

* Upgrade DF to fulldep

* Minor fix

* Equality of OrbitFit

* Adaptative maxsteps

* HTTP/JSON for get_radec_mpc

* Fix get_radec_mpc

* Fix adaptative maxsteps

* Fix ra discontinuity

* Softer _gaussinitcond Val(2) condition

* Fix unfold and nrms

* Fix Val(2) _gaussinitcond

* New adaptative maxsteps

* Yet another adaptative maxsteps

* Avoid loading debiasing table more than once

* Fix signatures

* Fix more signatures

* Update CatalogueMPC to new observations API

* Update ObservatoryMPC to new observations API

* Minor fix

* Update RadecMPC to new observations API

* Minor fixes

* Try previous maxsteps

* Fix reduce_nights

* Yet another adaptative maxsteps

* Fix tests

* Fix obstech

* Fix reduce_nights

* Fix discontinuity in reduce_nights

* Fix residuals computation

* Try adaptative triplets

* Fix tests

* Add residual rejection

* Remove nyears as orbitdet args

* Retry residual rejection

* Fix tests

* Fix more tests

* Add preliminary gauss tests

* Fix Gauss tests

* Comment out residual rejection

* Outlier rejection by mean point

* Minor fixes

* Another minor fix

* Clustering outlier rejection criteria

* K-means clustering

* Minor fix

* Too Short Arc

* issinglearc

* Fixes to tooshortarc

* Automatic scaling

* Tiny object boundary

* Minor fixes

* ObservationNight

* Fixes to Gauss Method

* Minor fix

* Update Gauss Method tests

* Update Too Short Arc

* Update TSA tests

* Update Outlier Rejection

* Update Outlier Rejection tests

* Minor fix

* tsanightorder & naive maxsteps

* Fixes to TSA

* Fix istsa

* Fix boundary_projection

* Add init to kmeans

* Minor fixes

* Polynomial fit

* Heel anomaly in tsa

* Minor fix

* Another minor fix

* unsuccessful_propagation

* Try reset CI.yml to origin/main

* Reset .github/ and data/ to origin/main

* Reset dev/ and notebooks/ to origin/main

* Reset pha/ and scripts/ to origin/main

* Reset NEOs.jl/ to origin/main

* Reset test/ to origin/main

* Reset src/ to origin/main

* Remove .DS_Store

* OrbitFit -> LeastSquaresFit

* Parameters -> NEOParameters

* issuccessfulprop

* ObservationNight -> Tracklet

* Change all night refs to tracklet

* Add 2014 AA test

* Add 2008 TC3 test

* Minor fix

* Reduce fwdoffset in 2008 TC3 test

* Minor fix in constants.jl

* Update HORIZONS and TS

* Update TaylorInterpolant

* Minor fix

* Add eph_su(ea) to NEOParameters

* Float64 -> T in AdmissibleRegion

* Minor fixes

* Light-time correction

* Minor fix

* Incremental radec in TSA

* Minor fixes

* Add Newtonian model as default dynamics for preliminary orbit determination (#53)

* Add newtonian dynamical model

* Use Newtonian dynamical model throughout preliminary orbit determination methods

* Make newtonian! default dynamical model for preliminary orbit determination

---------

Co-authored-by: LuEdRaMo <[email protected]>

* fetch_radec_mpc

* H_max to NEOParameters

* nms(::NEOSolution)

* log adam

* Isolate outlier_rejection

* Qmax -> NEOParameters

* ecef_to_geocentric(::Vector{TaylorN}})

* Require julia v1.9+

* Bump minor version

* Minor fixes

* Add a_max to NEOParameters

* Minor fix

* Another minor fix

* Yet more minor fixes

* Levenberg-Marquardt

* Update OD tests

* Minor fix

* Some aesthetic fixes (#54)

* Change end of line sequence

* Add NEOParameters test

* NEOParameters fix

* Try mode change

* More mode changes

* More mode changes

* Add remaining mode changes

* Update least_squares.jl

* Revert undesired changes

* Update to TaylorSeries 0.17 (#58)

* Update to TaylorSeries 0.17

* Formatting

* Fix tests

* Custom print for AdmissibleRegion

* Change Gauss/TSA order in OD

* jplcompare

* Update Apophis script

* Custom print for NEOParameters

* Update .gitignore

* Update some comments and docstrings

* Fix initial guess in LS

* Allow jd0::U in tooshortarc

* Allow jdo::U in propagation

* Fix first guess in LS

* propagate_params(..., ::NEOParameters{T})

* Add some internal variables as kwargs in methods

* Adapt diffcorr for idxs

* Reduce integrations in adam

* Minor fix

* orbitdetermination(::Vector{RadecMPC}, ::NEOSolution,...)

* Minor fix

* Test 2008 TC3 sigmas(sol1)

* Spelling corrections in test/orbit_determination.jl

* test/extensions.jl -> test/dataframe.jl

* Rename test/dataframe.jl -> test/dataframes.jl

* Add commented test and TODO

* Remove dep (Tables) in tests environment

---------

Co-authored-by: Luis Eduardo Ramírez Montoya <[email protected]>
Co-authored-by: Jorge Pérez <[email protected]>
Co-authored-by: Jorge A. Pérez-Hernández <[email protected]>

.github/workflows/CI.yml

@@ -4,10 +4,6 @@ on:
     paths-ignore:
       - 'LICENSE.md'
       - 'README.md'
-  pull_request:
-    branches:
-      - main
-    tags: '*'
 concurrency:
   # Skip intermediate builds: always.
   # Cancel intermediate builds: only if it is a pull request build.
@@ -27,7 +23,6 @@ jobs:
           - 'nightly'
         os:
           - ubuntu-latest
-          - macos-latest
           - windows-latest
         arch:
           - x64
@@ -43,9 +38,10 @@ jobs:
         env:
           JULIA_NUM_THREADS: 2
       - uses: julia-actions/julia-processcoverage@v1
-      - uses: codecov/codecov-action@v3
+      - uses: codecov/codecov-action@v4
         with:
-          files: lcov.info
+          token: ${{ secrets.CODECOV_TOKEN }}
+          fail_ci_if_error: false
       - uses: coverallsapp/github-action@v2
         with:
           path-to-lcov: lcov.info

.gitignore

@@ -8,3 +8,5 @@ jldeph/
 notebooks/figs
 EOP_IAU*.TXT
 Manifest.toml
+*.DS_Store
+.vscode

Project.toml

@@ -1,11 +1,13 @@
 name = "NEOs"
 uuid = "b41c07a2-2abb-11e9-070a-c3c1b239e7df"
 authors = ["Jorge A. Pérez Hernández", "Luis Benet", "Luis Eduardo Ramírez Montoya"]
-version = "0.7.5"
+version = "0.8.0"
 
 [deps]
 Artifacts = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"
 AutoHashEquals = "15f4f7f2-30c1-5605-9d31-71845cf9641f"
+Clustering = "aaaa29a8-35af-508c-8bc3-b662a17a0fe5"
+DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
 DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
 Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b"
@@ -14,55 +16,44 @@ HORIZONS = "5a3ac768-beb4-554a-9c98-3342fe3377f5"
 HTTP = "cd3eb016-35fb-5094-929b-558a96fad6f3"
 Healpix = "9f4e344d-96bc-545a-84a3-ae6b9e1b672b"
 InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
-IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253"
-IntervalRootFinding = "d2bf35a9-74e0-55ec-b149-d360ff49b807"
 JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
+JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
 LazyArtifacts = "4af54fe1-eca0-43a8-85a7-787d91b784e3"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+LsqFit = "2fda8390-95c7-5789-9bda-21331edee243"
 PlanetaryEphemeris = "d83715d0-7e5f-11e9-1a59-4137b20d8363"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Quadmath = "be4d8f0f-7fa4-5f49-b795-2f01399ab2dd"
-Requires = "ae029012-a4dd-5104-9daa-d747884805df"
+Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
 SPICE = "5bab7191-041a-5c2e-a744-024b9c3a5062"
 SatelliteToolboxTransformations = "6b019ec1-7a1e-4f04-96c7-a9db1ca5514d"
 Scratch = "6c6a2e73-6563-6170-7368-637461726353"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
+Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
 TaylorIntegration = "92b13dbe-c966-51a2-8445-caca9f8a7d42"
 TaylorSeries = "6aa5eb33-94cf-58f4-a9d0-e4b2c4fc25ea"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
-[weakdeps]
-DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
-Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
-
-[extensions]
-DataFramesExt = ["DataFrames", "Tables"]
-
 [compat]
 Artifacts = "1"
 AutoHashEquals = "0.2"
+Clustering = "0.15"
 DataFrames = "1.5"
 DelimitedFiles = "1"
 Downloads = "1"
-HORIZONS = "0.3"
+HORIZONS = "0.4"
 HTTP = "1.9.5"
 Healpix = "4"
-IntervalArithmetic = "0.20, 0.22"
-IntervalRootFinding = "0.5.11"
 JLD2 = "0.4"
+JSON = "0.21"
 LazyArtifacts = "1"
-PlanetaryEphemeris = "0.7"
+PlanetaryEphemeris = "0.8"
 Quadmath = "0.5"
-Requires = "0.5.2, 1"
 SPICE = "0.2"
 SatelliteToolboxTransformations = "0.1"
 Scratch = "1.2"
 StatsBase = "0.33, 0.34"
 Tables = "1.10"
-TaylorIntegration = "0.14"
-TaylorSeries = "0.16"
-julia = "1.6"
-
-[extras]
-DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
-Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
+TaylorIntegration = "0.15"
+TaylorSeries = "0.17"
+julia = "1.9"

pha/Project.toml

@@ -14,4 +14,4 @@ StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 TaylorIntegration = "92b13dbe-c966-51a2-8445-caca9f8a7d42"
 
 [compat]
-NEOs = "0.7"
+NEOs = "0.8"

pha/apophis.jl

@@ -124,41 +124,37 @@ function main(dynamics::D, maxsteps::Int, jd0_datetime::DateTime, nyears_bwd::T,
     jd0 = datetime2julian(jd0_datetime)
 
     print_header("Integrator warmup", 2)
-    _ = NEOs.propagate(dynamics, 1, jd0, nyears_fwd, q0, Val(true);
-                       order, abstol, parse_eqs)
-    _ = NEOs.propagate_root(dynamics, 1, jd0, nyears_fwd, q0, Val(true);
-                       order, abstol, parse_eqs)
+    params = NEOParameters(;maxsteps=1, order, abstol, parse_eqs)
+    _ = NEOs.propagate(dynamics, jd0, nyears_fwd, q0, params)
+    _ = NEOs.propagate_root(dynamics, jd0, nyears_fwd, q0, params)
     println()
 
     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))
 
-    sol_bwd = NEOs.propagate(dynamics, maxsteps, jd0, nyears_bwd, q0, Val(true);
-                             order, abstol, parse_eqs)
+    params = NEOParameters(;maxsteps, order, abstol, parse_eqs)
+    sol_bwd = NEOs.propagate(dynamics, jd0, nyears_bwd, q0, params)
     jldsave("Apophis_bwd.jld2"; sol_bwd)
     # sol_bwd = JLD2.load("Apophis_bwd.jld2", "sol_bwd")
 
     tmax = nyears_fwd*yr
     println("• Initial time of integration: ", string(jd0_datetime))
     println("• Final time of integration: ", julian2datetime(jd0 + tmax))
 
-    sol_fwd, tvS, xvS, gvS = NEOs.propagate_root(dynamics, maxsteps, jd0, nyears_fwd, q0, Val(true);
-                             order, abstol, parse_eqs)
+    sol_fwd, tvS, xvS, gvS = NEOs.propagate_root(dynamics, jd0, nyears_fwd, q0, params)
     jldsave("Apophis_fwd.jld2"; sol_fwd, tvS, xvS, gvS)
     # sol_fwd = JLD2.load("Apophis_fwd.jld2", "sol_bwd")
     println()
 
-    # load Solar System ephemeris
-    sseph::TaylorInterpolant{Float64,Float64,2} = loadpeeph(NEOs.sseph, sol_bwd.t0+sol_bwd.t[end], sol_fwd.t0+sol_fwd.t[end])
-    eph_su::TaylorInterpolant{Float64,Float64,2} = selecteph(sseph, su)
-    eph_ea::TaylorInterpolant{Float64,Float64,2} = selecteph(sseph, ea)
+    # Load Sun ephemeris
+    eph_su = selecteph(NEOs.sseph, su, sol_bwd.t0+sol_bwd.t[end], sol_fwd.t0+sol_fwd.t[end])
+    # Load Earth ephemeris
+    eph_ea = selecteph(NEOs.sseph, ea, sol_bwd.t0+sol_bwd.t[end], sol_fwd.t0+sol_fwd.t[end])
 
-    # NEO
+    # Apophis
     # Change t, x, v units, resp., from days, au, au/day to sec, km, km/sec
-    xva_bwd(et) = auday2kmsec(sol_bwd(et/daysec)[1:6])
-    xva_fwd(et) = auday2kmsec(sol_fwd(et/daysec)[1:6])
     xva(et) = bwdfwdeph(et, sol_bwd, sol_fwd)
     # Earth
     # Change x, v units, resp., from au, au/day to km, km/sec
@@ -177,7 +173,8 @@ function main(dynamics::D, maxsteps::Int, jd0_datetime::DateTime, nyears_bwd::T,
     deldop = vcat(deldop_2005_2013,deldop_2021)
 
     # Compute optical residuals
-    res_radec_all, w_radec_all = NEOs.residuals(radec; xvs, xve, xva)
+    _res_radec_all_ = NEOs.residuals(radec; xvs, xve, xva)
+    res_radec_all, w_radec_all = NEOs.unfold(_res_radec_all_)
     jldsave("Apophis_res_w_radec.jld2"; res_radec_all, w_radec_all)
     # JLD2.@load "Apophis_res_w_radec.jld2"
 
@@ -233,9 +230,9 @@ function main(dynamics::D, maxsteps::Int, jd0_datetime::DateTime, nyears_bwd::T,
     res = vcat(res_radec, res_del, res_dop)
     w = vcat(w_radec, w_del, w_dop)
 
-    success, δx_OR8, Γ_OR8 = newtonls(res, w, zeros(get_numvars()), 10)
-    x_OR8 = sol_fwd(sol_fwd.t0)(δx_OR8)
-    σ_OR8 = sqrt.(diag(Γ_OR8)).*scalings
+    fit_OR8 = newtonls(res, w, zeros(get_numvars()), 10)
+    x_OR8 = sol_fwd(sol_fwd.t0)(fit_OR8.x)
+    σ_OR8 = sqrt.(diag(fit_OR8.Γ)).*scalings
 
     nradec = length(res_radec)
     res_ra = view(res_radec, 1:nradec÷2)
@@ -248,19 +245,19 @@ function main(dynamics::D, maxsteps::Int, jd0_datetime::DateTime, nyears_bwd::T,
     print_header("Orbital fit (8-DOF) and post-fit statistics", 2)
 
     # orbital fit
-    println("Success flag                                                      : ", success, "\n")
+    println("Success flag                                                      : ", fit_OR8.success, "\n")
     println("Nominal solution             [au,au,au,au/d,au/d,au/d,au/d²,au/d²]: ", x_OR8, "\n")
     println("1-sigma formal uncertainties [au,au,au,au/d,au/d,au/d,au/d²,au/d²]: ", σ_OR8, "\n")
 
     # post-fit statistics
-    println("Normalized RMS (optical-only)               [adimensional] : ", nrms(res_radec(δx_OR8),w_radec))
-    println("Normalized RMS (radar-only)                 [adimensional] : ", nrms(vcat(res_del,res_dop)(δx_OR8),vcat(w_del,w_dop)))
-    println("Normalized RMS (combined optical and radar) [adimensional] : ", nrms(res(δx_OR8),w), "\n")
-    println("Mean weighted right-ascension residual      [arcseconds]   : ", mean(res_ra(δx_OR8), weights(w_ra)))
-    println("Mean weighted declination residual          [arcseconds]   : ", mean(res_dec(δx_OR8), weights(w_dec)))
-    println("Mean weighted time-delay residual           [micro-seconds]: ", mean(res_del(δx_OR8), weights(w_del)))
-    println("Mean weighted Doppler-shift residual        [Hz]           : ", mean(res_dop(δx_OR8), weights(w_dop)), "\n")
-    println("Chi-squared statistic (χ²):                 [adimensional] : ", chi2(res(δx_OR8),w))
+    println("Normalized RMS (optical-only)               [adimensional] : ", nrms(res_radec(fit_OR8.x),w_radec))
+    println("Normalized RMS (radar-only)                 [adimensional] : ", nrms(vcat(res_del,res_dop)(fit_OR8.x),vcat(w_del,w_dop)))
+    println("Normalized RMS (combined optical and radar) [adimensional] : ", nrms(res(fit_OR8.x),w), "\n")
+    println("Mean weighted right-ascension residual      [arcseconds]   : ", mean(res_ra(fit_OR8.x), weights(w_ra)))
+    println("Mean weighted declination residual          [arcseconds]   : ", mean(res_dec(fit_OR8.x), weights(w_dec)))
+    println("Mean weighted time-delay residual           [micro-seconds]: ", mean(res_del(fit_OR8.x), weights(w_del)))
+    println("Mean weighted Doppler-shift residual        [Hz]           : ", mean(res_dop(fit_OR8.x), weights(w_dop)), "\n")
+    println("Chi-squared statistic (χ²):                 [adimensional] : ", chi2(res(fit_OR8.x),w))
 
     return sol_bwd, sol_fwd, res_radec, res_del, res_dop, w_radec, w_del, w_dop
 end
@@ -316,3 +313,36 @@ function main()
 end
 
 main()
+
+#=
+-0.18034827489412805
+0.9406910783153754
+0.3457360118643932
+-0.016265940057745887
+4.3915296805381036e-5
+-0.0003952032399008921
+-2.883842658925719e-14
+-1.6977850502837174e-11
+
+
+6.991681792970175e-9
+3.3210758500485412e-9
+8.763895947794195e-9
+5.485729074188412e-11
+9.846074910774499e-11
+1.3782767177665416e-10
+2.525644870489942e-16
+3.2193289354129716e-12
+
+Normalized RMS (optical-only)               [adimensional] : 0.4365429773707507
+Normalized RMS (radar-only)                 [adimensional] : 0.5421268387032326
+Normalized RMS (combined optical and radar) [adimensional] : 0.43685679642513925
+
+Mean weighted right-ascension residual      [arcseconds]   : -0.002787272519210798
+Mean weighted declination residual          [arcseconds]   : -0.003376099885054764
+Mean weighted time-delay residual           [micro-seconds]: -0.0007329351962194035
+Mean weighted Doppler-shift residual        [Hz]           : -0.03025601266756564
+
+Chi-squared statistic (χ²):                 [adimensional] : 3597.4067719864497
+
+=#