Add preliminary docstrings

ext/DataFramesExt.jl

@@ -5,7 +5,7 @@ using TaylorSeries: get_numvars
 using PlanetaryEphemeris: J2000, selecteph, su, ea, yr, daysec, auday2kmsec
 using NEOs: RadecMPC, date, gauss_triplets, propagate, RNp1BP_pN_A_J23E_J2S_eph_threads!, order, abstol, sseph,
             scaled_variables, gauss_method, residuals, bwdfwdeph, newtonls, diffcorr, nrms, hascoord, tryls,
-            TimeOfDay
+            TimeOfDay, OpticalResidual, heliocentric_energy
 
 import Base: convert
 import NEOs: AbstractAstrometry, extrapolation, reduce_nights, gaussinitcond, relax_factor
@@ -87,6 +87,12 @@ function reduce_nights(radec::Vector{RadecMPC{T}}) where {T <: AbstractFloat}
     return cdf.observatory, cdf.date, cdf.α, cdf.δ
 end 
 
+@doc raw"""
+    relax_factor(radec::Vector{RadecMPC{T}}, ξs::Vector{OpticalResidual{T}}) where {T <: Real}
+
+Modify the weights in `ξs` by a relaxing factor quantifying the correlation between observations taken on the 
+same night by the same observatory.
+"""
 function relax_factor(radec::Vector{RadecMPC{T}}) where {T <: AbstractFloat}
     # Convert to DataFrame 
     df = DataFrame(radec)
@@ -101,8 +107,14 @@ function relax_factor(radec::Vector{RadecMPC{T}}) where {T <: AbstractFloat}
     return map(x -> x > 4.0 ? x/4.0 : 1.0, Nv)
 end 
 
+function relax_factor(radec::Vector{RadecMPC{T}}, ξs::Vector{OpticalResidual{T, U}}) where {T <: Real, U <: Number}
+    rex = relax_factor(radec)
+    return OpticalResidual.(getfield.(ξs, :ξ_α), getfield.(ξs, :ξ_δ), getfield.(ξs, :w_α), getfield.(ξs, :w_δ), rex,
+                            getfield.(ξs, :outlier))
+end
+
 @doc raw"""
-    gaussinitcond(radec::Vector{RadecMPC{T}}; Δ::Period = Day(1), Q_max::T = 0.75, niter::Int = 5, maxsteps::Int = 100, 
+    gaussinitcond(radec::Vector{RadecMPC{T}}; max_triplets::Int = 10, Q_max::T = 100., niter::Int = 5, maxsteps::Int = 100, 
                   varorder::Int = 5, order::Int = order, abstol::T = abstol, parse_eqs::Bool = true) where {T <: AbstractFloat}
 
 Return initial conditions via Gauss method. 
@@ -111,7 +123,7 @@ See also [`gauss_method`](@ref).
 
 # Arguments
 - `radec::Vector{RadecMPC{T}}`: vector of observations.
-- `Δ::Period`: see [`gauss_triplets`](@ref).
+- `max_triplets::Int`: maximum number of triplets.
 - `Q_max::T`: The maximum nrms that is considered a good enough orbit.
 - `niter::Int`: number of iterations for Newton's method.
 - `maxsteps::Int`: maximum number of steps for propagation.
@@ -169,6 +181,9 @@ function gaussinitcond(radec::Vector{RadecMPC{T}}; max_triplets::Int = 10, Q_max
         # Gauss method solution 
         sol = gauss_method(observatories[triplet], dates[triplet], α[triplet] .+ dq[1:3], δ[triplet] .+ dq[4:6]; niter = niter)
 
+        # Filter Gauss solutions by heliocentric energy
+        filter!(x -> heliocentric_energy(x.statevect) <= 0, sol)
+
         # Iterate over Gauss solutions
         for i in eachindex(sol)
 
@@ -180,23 +195,20 @@ function gaussinitcond(radec::Vector{RadecMPC{T}}; max_triplets::Int = 10, Q_max
                                 order = order, abstol = abstol, parse_eqs = parse_eqs)
 
             # O-C residuals
-            res, w = residuals(radec; xvs = et -> auday2kmsec(eph_su(et/daysec)), xve = et -> auday2kmsec(eph_ea(et/daysec)), 
-                               xva = et -> bwdfwdeph(et, bwd, fwd))
+            res = residuals(radec; xvs = et -> auday2kmsec(eph_su(et/daysec)), xve = et -> auday2kmsec(eph_ea(et/daysec)), 
+                            xva = et -> bwdfwdeph(et, bwd, fwd))
 
             # Subset of radec for orbit fit
-            i_O = findall(x -> dates[triplet[1]] <= date(x) <= dates[triplet[3]], radec)
-            i_ξ = vcat(i_O, i_O .+ length(radec))
+            idxs = findall(x -> dates[triplet[1]] <= date(x) <= dates[triplet[3]], radec)
 
             # Relaxation factor
-            rex = relax_factor(radec)
-            rex = vcat(rex, rex)
-            w = w ./ rex
-
+            res = relax_factor(radec, res)
+
             # Orbit fit
-            fit = tryls(res[i_ξ], w[i_ξ], zeros(get_numvars()), niter)
+            fit = tryls(res[idxs], zeros(get_numvars()), niter)
 
             # NRMS
-            Q = nrms(res(fit.x), w)
+            Q = nrms(res, fit)
 
             # TO DO: check cases where newton converges but diffcorr no
 

src/observations/process_radec.jl

@@ -1,3 +1,70 @@
+@doc raw"""
+    OpticalResidual{T <: Real, U <: Number}
+
+An astrometric optical observed minus computed residual. 
+
+# Fields
+
+- `ξ_α::U`: right ascension residual.
+- `ξ_δ::U`: declination residual.
+- `w_α::T`: right ascension weight.
+- `w_δ::T`: declination weight.
+- `relax_factor::T`: relaxation factor.
+- `outlier::Bool`: whether the residual is an outlier or not.
+"""
+@auto_hash_equals struct OpticalResidual{T <: Real, U <: Number}
+    ξ_α::U
+    ξ_δ::U
+    w_α::T
+    w_δ::T
+    relax_factor::T
+    outlier::Bool
+    function OpticalResidual{T, U}(ξ_α::U, ξ_δ::U, w_α::T, w_δ::T, relax_factor::T = one(T),
+                                   outlier::Bool = false) where {T <: Real, U <:  Number}
+        new{T, U}(ξ_α, ξ_δ, w_α, w_δ, relax_factor, outlier)
+    end
+end
+
+function OpticalResidual(ξ_α::U, ξ_δ::U, w_α::T, w_δ::T, relax_factor::T = one(T),
+                         outlier::Bool = false) where {T <: Real, U <: Number}
+    return OpticalResidual{T, U}(ξ_α, ξ_δ, w_α, w_δ, relax_factor, outlier)
+end
+
+# Evaluate methods
+function evaluate(res::OpticalResidual{T, TaylorN{T}}, x::Vector{T}) where {T <: Real}
+    return OpticalResidual(res.ξ_α(x), res.ξ_δ(x), res.w_α, res.w_δ, res.relax_factor, res.outlier)
+end
+(res::OpticalResidual{T, TaylorN{T}})(x::Vector{T}) where {T <: Real} = evaluate(res, x)
+
+function evaluate(res::Vector{OpticalResidual{T, TaylorN{T}}}, x::Vector{T}) where {T <: Real}
+    res_new = Vector{OpticalResidual{T, T}}(undef, length(res))
+    for i in eachindex(res)
+        res_new[i] = evaluate(res[i], x)
+    end
+    return res_new
+end
+(res::AbstractVector{OpticalResidual{T, TaylorN{T}}})(x::Vector{T}) where {T <: Real} = evaluate(res, x)
+
+# Print method for OpticalResidual
+# Examples:
+# α: -138.7980136773549 δ: -89.8002527255012
+# α: -134.79449984291568 δ: -91.42509376643284
+function show(io::IO, x::OpticalResidual{T, U}) where {T <: Real, U <: Number}
+    print(io, "α: ", cte(x.ξ_α), " δ: ", cte(x.ξ_δ))
+end
+
+@doc raw"""
+    unfold(ξs::Vector{OpticalResidual{T, U}}) where {T <: Real, U <: Number}
+
+Concatenate right ascension and declination residuals for an orbit fit.
+"""
+function unfold(ξs::Vector{OpticalResidual{T, U}}) where {T <: Real, U <: Number}
+    res = vcat(getfield.(ξs, :ξ_α), getfield.(ξs, :ξ_δ))
+    relax_factor = getfield.(ξs, :relax_factor)
+    w = vcat(getfield.(ξs, :w_α) ./ relax_factor, getfield.(ξs, :w_δ) ./ relax_factor)
+    return res, w
+end
+
 @doc raw"""
     compute_radec(observatory::ObservatoryMPC{T}, t_r_utc::DateTime; niter::Int = 5,
     xve::EarthEph = earthposvel, xvs::SunEph = sunposvel, xva::AstEph) where {T <: AbstractFloat, EarthEph, SunEph, AstEph}
@@ -493,17 +560,12 @@ See also [`compute_radec`](@ref), [`debiasing`](@ref), [`w8sveres17`](@ref) and
 - `xva::AstEph`: asteroid ephemeris [et seconds since J2000] -> [barycentric position in km and velocity in km/sec].
 """
 function residuals(obs::Vector{RadecMPC{T}}; kwargs...) where {T <: AbstractFloat}
-
     # Optical astrometry (dates + observed + computed + debiasing + weights)
     x_jt = radec_astrometry(obs; kwargs...)
     # Right ascension residuals
     res_α = x_jt[2] .- x_jt[6] .- x_jt[4]
     # Declination residuals
     res_δ = x_jt[3] .- x_jt[7] .- x_jt[5]
     # Total residuals
-    res = vcat(res_α, res_δ)
-    # Weights
-    w = repeat(1 ./ x_jt[end].^2, 2)
-
-    return res, w
-end
+    return OpticalResidual.(res_α, res_δ, 1 ./ x_jt[end].^2, 1 ./ x_jt[end].^2, one(T), false)
+end

src/observations/topocentric.jl

@@ -1,6 +1,17 @@
 # Earth orientation parameters (eop) 2000
 const eop_IAU2000A::EopIau2000A = fetch_iers_eop(Val(:IAU2000A))
 
+@doc raw"""
+    TimeOfDay
+
+Day/night at a particular timezone.
+
+# Fields 
+- `light::Symbol`: `:day` or `:night`.
+- `start::Date`.
+- `stop::Date`.
+- `utc::Int`: hours from UTC.
+"""
 @auto_hash_equals struct TimeOfDay
     light::Symbol
     start::Date
@@ -27,23 +38,35 @@ const eop_IAU2000A::EopIau2000A = fetch_iers_eop(Val(:IAU2000A))
     end
 end
 
+TimeOfDay(radec::RadecMPC{T}) where {T <: AbstractFloat} = TimeOfDay(date(radec), observatory(radec))
+
 isday(x::TimeOfDay) = x.light == :day
 isnight(x::TimeOfDay) = x.light == :night
 
-# Print method for RadecMPC
+# Print method for TimeOfDay
 # Examples:
-# N00hp15 α: 608995.65 δ: -25653.3 t: 2020-12-04T10:41:43.209 obs: 703
-# 99942 α: 422475.3 δ: 97289.49 t: 2021-05-12T06:28:35.904 obs: F51
+# Night from 2023-06-29 to 2023-06-29 at UTC-7
+# Night from 2023-06-29 to 2023-06-30 at UTC+3
 function show(io::IO, m::TimeOfDay)
     print(io, uppercasefirst(string(m.light)), " from ", m.start, " to ", m.stop, " at UTC", @sprintf("%+d", m.utc))
 end
 
-hours_from_UTC(lon::T) where {T <: Number} = ceil(Int, 12*lon/π - 0.5)
+@doc raw"""
+    hours_from_UTC(lon::T) where {T <: Real}
+
+Return the naive hour difference between longitude `lon` [rad] and UTC.
+"""
+hours_from_UTC(lon::T) where {T <: Real} = ceil(Int, 12*lon/π - 0.5)
 function hours_from_UTC(observatory::ObservatoryMPC{T}) where  {T <: AbstractFloat} 
     lon, _ = lonlat(observatory)
     return hours_from_UTC(lon)
 end
 
+@doc raw"""
+    lonlat(observatory::ObservatoryMPC{T}) where {T <: AbstractFloat}
+
+Return longitude and latitude (both in rad) of an observatory.
+"""
 function lonlat(observatory::ObservatoryMPC{T}) where {T <: AbstractFloat}
     # ECEF [km]
     p_ECEF = obs_pos_ECEF(observatory)
@@ -55,6 +78,14 @@ function lonlat(observatory::ObservatoryMPC{T}) where {T <: AbstractFloat}
     return lon, lat_geocentric
 end 
 
+@doc raw"""
+    sunriseset(date::DateTime, observatory::ObservatoryMPC{T}) where {T <: AbstractFloat}
+
+Return `DateTime` of sunrise and sunset at a particular date and location.
+
+!!! reference
+    See "General Solar Position Calculations" by NOAA at https://gml.noaa.gov/grad/solcalc/solareqns.PDF.
+"""
 function sunriseset(date::DateTime, observatory::ObservatoryMPC{T}) where {T <: AbstractFloat}
     # Fractional year [rad]
     γ = 2π * (dayofyear(date) - 1 + (hour(date)-12)/24) / daysinyear(date)
@@ -81,6 +112,8 @@ function sunriseset(date::DateTime, observatory::ObservatoryMPC{T}) where {T <:
     return sunrise, sunset
 end
 
+sunriseset(radec::RadecMPC{T}) where {T <: AbstractFloat} = sunriseset(date(radec), observatory(radec))
+
 @doc raw"""
     obs_pos_ECEF(observatory::ObservatoryMPC{T}) where {T <: AbstractFloat}
     obs_pos_ECEF(x::RadecMPC{T}) where {T <: AbstractFloat}

src/orbit_determination/gauss_method.jl

@@ -3,8 +3,7 @@ include("osculating.jl")
 @doc raw"""
     GaussSolution{T <: Real, U <: Number}
 
-A preliminary orbit obtained from Gauss method of orbit determination. See Algorithm 5.5 in page 274 of
-https://doi.org/10.1016/C2016-0-02107-1.
+A preliminary orbit obtained from Gauss method of orbit determination.
 
 See also [`gauss_method`](@ref).
 
@@ -18,6 +17,9 @@ See also [`gauss_method`](@ref).
 - `τ_3::T`: time between third and second observations.
 - `f_1, g_1, f_3, g_3::U`: Lagrange coefficients.
 - `status::Symbol`: the status of the solution (`:empty`, `:unkown` or `:unique`)
+
+!!! reference
+    See Algorithm 5.5 in page 274 of https://doi.org/10.1016/C2016-0-02107-1.
 """
 @auto_hash_equals struct GaussSolution{T <: Real, U <: Number}
     statevect::Vector{U}
@@ -48,7 +50,7 @@ end
 # Examples:
 # unique Gauss solution (r = 1.0800950907383229 AU)
 function show(io::IO, g::GaussSolution{T, U}) where {T <: Real, U <: Number}
-    print(io, g.status, " Gauss solution (r = ", norm(cte.(g.statevect[1:3])), " AU)")
+    print(io, uppercasefirst(string(g.status)), " Gauss solution (r = ", norm(cte.(g.statevect[1:3])), " AU)")
 end
 
 @doc raw"""
@@ -178,7 +180,7 @@ end
     gauss_method(observatories::Vector{ObservatoryMPC{T}}, dates::Vector{DateTime}, α::Vector{U}, δ::Vector{U};
                  niter::Int = 5) where {T <: Real, U <: Number}
 
-Core Gauss method of Initial Orbit determination (IOD). See Algorithm 5.5 in page 274 https://doi.org/10.1016/C2016-0-02107-1.
+Core Gauss method of Initial Orbit determination (IOD).
 
 # Arguments
 
@@ -188,6 +190,9 @@ Core Gauss method of Initial Orbit determination (IOD). See Algorithm 5.5 in pag
 - `α::Vector{U}`: right ascension [rad].
 - `δ::Vector{U}`: declination [rad].
 - `niter::Int`: Number of iterations for Newton's method.
+
+!!! reference
+    See Algorithm 5.5 in page 274 https://doi.org/10.1016/C2016-0-02107-1.
 """
 function gauss_method(obs::Vector{RadecMPC{T}}; niter::Int = 5) where {T <: AbstractFloat}
 
@@ -317,9 +322,7 @@ function gauss_method(observatories::Vector{ObservatoryMPC{T}}, dates::Vector{Da
 
             sol_gauss[i] = GaussSolution{T, U}(vcat(r_vec[2, :], v_2_vec), D, R_vec, ρ_vec, τ_1, τ_3, f_1, g_1, f_3, g_3, status[i])
         end
-
-        filter!(x -> heliocentric_energy(x.statevect) <= 0, sol_gauss)
-
+        # Sort solutions by heliocentric range
         return sort!(sol_gauss, by = x -> norm(x.statevect[1:3]))
 
     end
@@ -351,14 +354,12 @@ for Gauss method. The function assumes `dates` is sorted.
 gauss_norm(dates::Vector{DateTime}) = abs( (dates[2] - dates[1]).value - (dates[3] - dates[2]).value )
 
 @doc raw"""
-    gauss_triplets(dates::Vector{DateTime}, Δ::Period = Day(1))
+    gauss_triplets(dates::Vector{DateTime}, max_triplets::Int = 10, max_iter::Int = 100)
 
-Return a vector of triplets to be used within [`gaussinitcond`](@ref) to select the best observations for Gauss method. 
+Return a vector of `max_triplets` triplets to be used within [`gaussinitcond`](@ref) to select the best observations for Gauss method. 
 The triplets are sorted by [`gauss_norm`](@ref).
-     
-See also [`closest_index_sorted`](@ref).
 """
-function gauss_triplets(dates::Vector{DateTime}, Δ_min::Period, Δ_max::Period, avoid::Vector{Vector{Int}}, max_triplets)
+function gauss_triplets(dates::Vector{DateTime}, Δ_min::Period, Δ_max::Period, avoid::Vector{Vector{Int}}, max_triplets::Int)
 
     triplets = Vector{Vector{Int}}(undef, 0)
     L = length(dates)