From 45e1b6b20c384215064edb1ceb996f63a35eee95 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Luis=20Eduardo=20Ram=C3=ADrez=20Montoya?= Date: Sun, 13 Aug 2023 15:19:37 -0600 Subject: [PATCH] Add preliminary docstrings --- ext/DataFramesExt.jl | 38 ++-- src/observations/process_radec.jl | 76 ++++++- src/observations/topocentric.jl | 41 +++- src/orbit_determination/gauss_method.jl | 25 +-- src/postprocessing/least_squares.jl | 148 ++++++++++--- src/propagation/jetcoeffs.jl | 276 ++++++++++++------------ 6 files changed, 405 insertions(+), 199 deletions(-) diff --git a/ext/DataFramesExt.jl b/ext/DataFramesExt.jl index c35cd876..c63a6d5e 100644 --- a/ext/DataFramesExt.jl +++ b/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 diff --git a/src/observations/process_radec.jl b/src/observations/process_radec.jl index 657efa5e..109037f2 100644 --- a/src/observations/process_radec.jl +++ b/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,7 +560,6 @@ 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 @@ -501,9 +567,5 @@ function residuals(obs::Vector{RadecMPC{T}}; kwargs...) where {T <: AbstractFloa # 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 \ No newline at end of file diff --git a/src/observations/topocentric.jl b/src/observations/topocentric.jl index e3d193cd..c70257fa 100644 --- a/src/observations/topocentric.jl +++ b/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} diff --git a/src/orbit_determination/gauss_method.jl b/src/orbit_determination/gauss_method.jl index 3224bd3c..30542b7e 100644 --- a/src/orbit_determination/gauss_method.jl +++ b/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) diff --git a/src/postprocessing/least_squares.jl b/src/postprocessing/least_squares.jl index 5e70f170..8dbb9f3c 100644 --- a/src/postprocessing/least_squares.jl +++ b/src/postprocessing/least_squares.jl @@ -1,5 +1,16 @@ include("b_plane.jl") +@doc raw""" + OrbitFit{T <: Real} + +A least squares fit. + +# Fields +- `success::Bool`: wheter the routine converged or not. +- `x::Vector{T}`: deltas that minimize the objective function. +- `Γ::Matrix{T}`: covariance matrix. +- `routine::Symbol`: minimization routine (`:newton` or `:diffcorr`). +""" struct OrbitFit{T <: Real} success::Bool x::Vector{T} @@ -12,14 +23,85 @@ OrbitFit(success::Bool, x::Vector{T}, Γ::Matrix{T}, routine::Symbol) where {T < # Print method for OrbitFit # 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 +# Succesful Newton +# Succesful differential corrections function show(io::IO, fit::OrbitFit{T}) where {T <: Real} success_s = fit.success ? "Succesful" : "Unsuccesful" routine_s = fit.routine == :newton ? "Newton" : "differential corrections" print(io, success_s, " ", routine_s) end +@doc raw""" + carpino_smoothing(n::T) where {T <: Real} + +Fudge term for rejection condition in [`outlier_rejection`](@ref). + +!!! reference + See page 253 of https://doi.org/10.1016/S0019-1035(03)00051-4. +""" +carpino_smoothing(n::T) where {T <: Real} = 400*(1.2)^(-n) + +@doc raw""" + outlier_rejection(ξs::Vector{OpticalResidual{T}}, fit::OrbitFit{T}; χ2_rec::T = 7., χ2_rej::T = 8., + α::T = 0.25) where {T <: Real} + +Outlier rejection algorithm. + +# Arguments + +- `ξs::Vector{OpticalResidual{T}}`: vector of residuals. +- `fit::OrbitFit{T}`: least squares fit. +- `χ2_rec::T`: recovery conditions. +- `χ2_rej::T`: rejection condition. +- `α::T`: scaling factor for maximum chi. + +!!! reference + See https://doi.org/10.1016/S0019-1035(03)00051-4. +""" +function outlier_rejection(ξs::Vector{OpticalResidual{T}}, fit::OrbitFit{T}; χ2_rec::T = 7., χ2_rej::T = 8., + α::T = 0.25) where {T <: Real} + # Number of residuals + N = length(ξs) + # Vector of chi2s + χ2s = Vector{T}(undef, N) + + for i in eachindex(χ2s) + # Current observation covariance matrix + γ = diagm([ξs[i].w_α/ξs[i].relax_factor, ξs[i].w_δ/ξs[i].relax_factor]) + # Current model matrix + A = hcat(TS.gradient(ξs[i].ξ_α)(fit.x), TS.gradient(ξs[i].ξ_δ)(fit.x)) + # Outlier sign + outlier_sign = ξs[i].outlier*2-1 + # Current residual covariance matrix + γ_ξ = γ + outlier_sign*(A')*fit.Γ*A + # Current residual + ξ = [ξs[i].ξ_α(fit.x), ξs[i].ξ_δ(fit.x)] + # Current chi2 + χ2s[i] = ξ' * inv(γ_ξ) * ξ + end + + χ2_max = maximum(χ2s) + new_ξs = Vector{OpticalResidual{T}}(undef, N) + N_sel = count(x -> !x.outlier, ξs) + + for i in eachindex(χ2s) + if χ2s[i] > max(χ2_rej + carpino_smoothing(N_sel), α*χ2_max) + new_ξs[i] = OpticalResidual(ξs[i].ξ_α, ξs[i].ξ_δ, ξs[i].w_α, ξs[i].w_δ, ξs[i].relax_factor, true) + elseif χ2s[i] < χ2_rec + new_ξs[i] = OpticalResidual(ξs[i].ξ_α, ξs[i].ξ_δ, ξs[i].w_α, ξs[i].w_δ, ξs[i].relax_factor, false) + else + new_ξs[i] = ξs[i] + end + end + + return new_ξs +end + +@doc raw""" + project(y::Vector{TaylorN{T}}, fit::OrbitFit{T}) where {T <: Real} + +Project `fit`'s covariance matrix into `y`. +""" function project(y::Vector{TaylorN{T}}, fit::OrbitFit{T}) where {T <: Real} J = Matrix{T}(undef, get_numvars(), length(y)) for i in eachindex(y) @@ -29,8 +111,9 @@ function project(y::Vector{TaylorN{T}}, fit::OrbitFit{T}) where {T <: Real} end @doc raw""" - nrms(res, w) - + nrms(res::Vector{U}, w::Vector{T}) where {T <: Real, U <: Number} + nrms(res::Vector{OpticalResidual{T}}, fit::OrbitFit{T}) where {T <: Real} + Returns the normalized root mean square error ```math \texttt{NRMS} = \sqrt{\frac{\chi^2}{m}}, @@ -44,15 +127,21 @@ See also [`chi2`](@ref). - `res`: Vector of residuals. - `w`: Vector of weights. +- `fit`: least squares fit. """ -function nrms(res, w) +function nrms(res::Vector{U}, w::Vector{T}) where {T <: Real, U <: Number} # Have as many residuals as weights @assert length(res) == length(w) # Normalized root mean square error return sqrt( chi2(res, w)/length(res) ) end +function nrms(res::Vector{OpticalResidual{T, U}}, fit::OrbitFit{T}) where {T <: Real, U <: Number} + res_, w = unfold(res) + return nrms(res_(fit.x), w) +end + @doc raw""" chi2(res, w) @@ -159,7 +248,7 @@ function ξTH(w, res, H_mat, npar) end @doc raw""" - diffcorr(res, w, x0, niters=5) + diffcorr(ξs::Vector{OpticalResidual{T}}, x0::Vector{T}, niters::Int = 5) where {T <: Real} Differential corrections subroutine for least-squares fitting. Returns the `niters`-th correction @@ -181,14 +270,13 @@ See also [`BHC`](@ref). # Arguments -- `res`: Vector of residuals. -- `w`: Vector of weights. -- `x_0`: First guess for the initial conditions. -- `niters`: Number of iterations. +- `ξs::Vector{OpticalResidual{T, TaylorN{T}}}`: Vector of residuals. +- `x_0::Vector{T}`: First guess for the initial conditions. +- `niters::Int`: Number of iterations. """ -function diffcorr(res::Vector{TaylorN{T}}, w::Vector{T}, x0::Vector{T}, niters::Int = 5) where {T <: Real} - # Have as many residuals as weights - @assert length(res) == length(w) +function diffcorr(ξs::Vector{OpticalResidual{T, TaylorN{T}}}, x0::Vector{T}, niters::Int = 5) where {T <: Real} + # Unfold residuals and weights + res, w = unfold(ξs) # Degrees of freedom npar = length(x0) # Design matrix B, H array and normal matrix C @@ -260,14 +348,13 @@ See also [`chi2`](@ref). # Arguments -- `res`: Vector of residuals. -- `w`: Vector of weights. -- `x_0`: First guess for the initial conditions. -- `niters`: Number of iterations. +- `ξs::Vector{OpticalResidual{T, TaylorN{T}}}`: Vector of residuals. +- `x_0::Vector{T}`: First guess for the initial conditions. +- `niters::Int`: Number of iterations. """ -function newtonls(res::Vector{TaylorN{T}}, w::Vector{T}, x0::Vector{T}, niters::Int = 5) where {T <: Real} - # Have as many residuals as weights - @assert length(res) == length(w) +function newtonls(ξs::Vector{OpticalResidual{T, TaylorN{T}}}, x0::Vector{T}, niters::Int = 5) where {T <: Real} + # Unfold residuals and weights + res, w = unfold(ξs) # Number of observations nobs = length(res) # Degrees of freedom @@ -321,14 +408,25 @@ function newtonls(res::Vector{TaylorN{T}}, w::Vector{T}, x0::Vector{T}, niters:: return OrbitFit(true, x_new, Γ, :newton) end -function tryls(res::Vector{TaylorN{T}}, w::Vector{T}, x0::Vector{T}, niters::Int = 5) where {T <: Real} +@doc raw""" + tryls(ξs::Vector{OpticalResidual{T}}, x0::Vector{T}, niters::Int = 5) where {T <: Real} + +Returns the best least squares fit between two routines: [`newtonls`](@ref) and [`diffcorr`](@ref). + +# Arguments + +- `ξs::Vector{OpticalResidual{T, TaylorN{T}}}`: Vector of residuals. +- `x_0::Vector{T}`: First guess for the initial conditions. +- `niters::Int`: Number of iterations. +""" +function tryls(ξs::Vector{OpticalResidual{T, TaylorN{T}}}, x0::Vector{T}, niters::Int = 5) where {T <: Real} # Newton's method - fit_1 = newtonls(res, w, x0, niters) + fit_1 = newtonls(ξs, x0, niters) # Differential corrections - fit_2 = diffcorr(res, w, x0, niters) + fit_2 = diffcorr(ξs, x0, niters) if fit_1.success && fit_2.success - Q_1 = nrms(res(fit_1.x), w) - Q_2 = nrms(res(fit_2.x), w) + Q_1 = nrms(ξs, fit_1) + Q_2 = nrms(ξs, fit_2) if Q_1 <= Q_2 return fit_1 else diff --git a/src/propagation/jetcoeffs.jl b/src/propagation/jetcoeffs.jl index 9547bbdc..f4990048 100644 --- a/src/propagation/jetcoeffs.jl +++ b/src/propagation/jetcoeffs.jl @@ -145,87 +145,87 @@ function TaylorIntegration._allocate_jetcoeffs!(::Val{RNp1BP_pN_A_J23E_J2S_ng_ep dq[3] = Taylor1(identity(constant_term(q[6])), order) newtonianNb_Potential[N] = Taylor1(identity(constant_term(zero_q_1)), order) tmp687 = Array{Taylor1{_S}}(undef, size(dq)) - tmp687 .= Taylor1(zero(_S), order) + tmp687 .= Taylor1(zero(constant_term(q[1])), order) tmp689 = Array{Taylor1{_S}}(undef, size(ui)) - tmp689 .= Taylor1(zero(_S), order) + tmp689 .= Taylor1(zero(constant_term(q[1])), order) tmp692 = Array{Taylor1{_S}}(undef, size(dq)) - tmp692 .= Taylor1(zero(_S), order) + tmp692 .= Taylor1(zero(constant_term(q[1])), order) tmp694 = Array{Taylor1{_S}}(undef, size(vi)) - tmp694 .= Taylor1(zero(_S), order) + tmp694 .= Taylor1(zero(constant_term(q[1])), order) tmp697 = Array{Taylor1{_S}}(undef, size(dq)) - tmp697 .= Taylor1(zero(_S), order) + tmp697 .= Taylor1(zero(constant_term(q[1])), order) tmp699 = Array{Taylor1{_S}}(undef, size(wi)) - tmp699 .= Taylor1(zero(_S), order) + tmp699 .= Taylor1(zero(constant_term(q[1])), order) pn2x = Array{Taylor1{_S}}(undef, size(X)) - pn2x .= Taylor1(zero(_S), order) + pn2x .= Taylor1(zero(constant_term(q[1])), order) pn2y = Array{Taylor1{_S}}(undef, size(Y)) - pn2y .= Taylor1(zero(_S), order) + pn2y .= Taylor1(zero(constant_term(q[1])), order) pn2z = Array{Taylor1{_S}}(undef, size(Z)) - pn2z .= Taylor1(zero(_S), order) + pn2z .= Taylor1(zero(constant_term(q[1])), order) tmp707 = Array{Taylor1{_S}}(undef, size(UU)) - tmp707 .= Taylor1(zero(_S), order) + tmp707 .= Taylor1(zero(constant_term(q[1])), order) tmp710 = Array{Taylor1{_S}}(undef, size(X)) - tmp710 .= Taylor1(zero(_S), order) + tmp710 .= Taylor1(zero(constant_term(q[1])), order) tmp712 = Array{Taylor1{_S}}(undef, size(Y)) - tmp712 .= Taylor1(zero(_S), order) + tmp712 .= Taylor1(zero(constant_term(q[1])), order) tmp713 = Array{Taylor1{_S}}(undef, size(tmp710)) - tmp713 .= Taylor1(zero(_S), order) + tmp713 .= Taylor1(zero(constant_term(q[1])), order) tmp715 = Array{Taylor1{_S}}(undef, size(Z)) - tmp715 .= Taylor1(zero(_S), order) + tmp715 .= Taylor1(zero(constant_term(q[1])), order) tmp723 = Array{Taylor1{_S}}(undef, size(pn2x)) - tmp723 .= Taylor1(zero(_S), order) + tmp723 .= Taylor1(zero(constant_term(q[1])), order) tmp724 = Array{Taylor1{_S}}(undef, size(tmp723)) - tmp724 .= Taylor1(zero(_S), order) + tmp724 .= Taylor1(zero(constant_term(q[1])), order) tmp819 = Array{Taylor1{_S}}(undef, size(ui)) - tmp819 .= Taylor1(zero(_S), order) + tmp819 .= Taylor1(zero(constant_term(q[1])), order) tmp821 = Array{Taylor1{_S}}(undef, size(vi)) - tmp821 .= Taylor1(zero(_S), order) + tmp821 .= Taylor1(zero(constant_term(q[1])), order) tmp822 = Array{Taylor1{_S}}(undef, size(tmp819)) - tmp822 .= Taylor1(zero(_S), order) + tmp822 .= Taylor1(zero(constant_term(q[1])), order) tmp824 = Array{Taylor1{_S}}(undef, size(wi)) - tmp824 .= Taylor1(zero(_S), order) + tmp824 .= Taylor1(zero(constant_term(q[1])), order) tmp735 = Array{Taylor1{_S}}(undef, size(X)) - tmp735 .= Taylor1(zero(_S), order) + tmp735 .= Taylor1(zero(constant_term(q[1])), order) tmp737 = Array{Taylor1{_S}}(undef, size(Y)) - tmp737 .= Taylor1(zero(_S), order) + tmp737 .= Taylor1(zero(constant_term(q[1])), order) tmp739 = Array{Taylor1{_S}}(undef, size(Z)) - tmp739 .= Taylor1(zero(_S), order) + tmp739 .= Taylor1(zero(constant_term(q[1])), order) tmp741 = Array{Taylor1{_S}}(undef, size(t31)) - tmp741 .= Taylor1(zero(_S), order) + tmp741 .= Taylor1(zero(constant_term(q[1])), order) tmp948 = Array{Taylor1{_S}}(undef, size(sin_ϕ)) - tmp948 .= Taylor1(zero(_S), order) + tmp948 .= Taylor1(zero(constant_term(q[1])), order) tmp949 = Array{Taylor1{_S}}(undef, size(ϕ)) - tmp949 .= Taylor1(zero(_S), order) + tmp949 .= Taylor1(zero(constant_term(q[1])), order) tmp751 = Array{Taylor1{_S}}(undef, size(sin2_ϕ)) - tmp751 .= Taylor1(zero(_S), order) + tmp751 .= Taylor1(zero(constant_term(q[1])), order) tmp757 = Array{Taylor1{_S}}(undef, size(sin_ϕ)) - tmp757 .= Taylor1(zero(_S), order) + tmp757 .= Taylor1(zero(constant_term(q[1])), order) tmp759 = Array{Taylor1{_S}}(undef, size(sin3_ϕ)) - tmp759 .= Taylor1(zero(_S), order) + tmp759 .= Taylor1(zero(constant_term(q[1])), order) tmp763 = Array{Taylor1{_S}}(undef, size(sin2_ϕ)) - tmp763 .= Taylor1(zero(_S), order) + tmp763 .= Taylor1(zero(constant_term(q[1])), order) tmp766 = Array{Taylor1{_S}}(undef, size(r_p2)) - tmp766 .= Taylor1(zero(_S), order) + tmp766 .= Taylor1(zero(constant_term(q[1])), order) tmp767 = Array{Taylor1{_S}}(undef, size(Λ2)) - tmp767 .= Taylor1(zero(_S), order) + tmp767 .= Taylor1(zero(constant_term(q[1])), order) tmp770 = Array{Taylor1{_S}}(undef, size(r_p1d2)) - tmp770 .= Taylor1(zero(_S), order) + tmp770 .= Taylor1(zero(constant_term(q[1])), order) tmp771 = Array{Taylor1{_S}}(undef, size(Λ3)) - tmp771 .= Taylor1(zero(_S), order) + tmp771 .= Taylor1(zero(constant_term(q[1])), order) tmp773 = Array{Taylor1{_S}}(undef, size(cos_ϕ)) - tmp773 .= Taylor1(zero(_S), order) + tmp773 .= Taylor1(zero(constant_term(q[1])), order) tmp775 = Array{Taylor1{_S}}(undef, size(cos_ϕ)) - tmp775 .= Taylor1(zero(_S), order) + tmp775 .= Taylor1(zero(constant_term(q[1])), order) tmp778 = Array{Taylor1{_S}}(undef, size(Λ2j_div_r4)) - tmp778 .= Taylor1(zero(_S), order) + tmp778 .= Taylor1(zero(constant_term(q[1])), order) tmp782 = Array{Taylor1{_S}}(undef, size(Λ3j_div_r5)) - tmp782 .= Taylor1(zero(_S), order) + tmp782 .= Taylor1(zero(constant_term(q[1])), order) tmp785 = Array{Taylor1{_S}}(undef, size(X)) - tmp785 .= Taylor1(zero(_S), order) + tmp785 .= Taylor1(zero(constant_term(q[1])), order) tmp787 = Array{Taylor1{_S}}(undef, size(Y)) - tmp787 .= Taylor1(zero(_S), order) + tmp787 .= Taylor1(zero(constant_term(q[1])), order) tmp789 = Array{Taylor1{_S}}(undef, size(Z)) - tmp789 .= Taylor1(zero(_S), order) + tmp789 .= Taylor1(zero(constant_term(q[1])), order) #= C:\Users\luisi\.julia\dev\NEOs.jl\cosa.jl:245 =# Threads.@threads for i = 1:Nm1 ui[i] = Taylor1(identity(constant_term(ss16asteph_t[3 * ((N - 1) + i) - 2])), order) vi[i] = Taylor1(identity(constant_term(ss16asteph_t[3 * ((N - 1) + i) - 1])), order) @@ -361,13 +361,13 @@ function TaylorIntegration._allocate_jetcoeffs!(::Val{RNp1BP_pN_A_J23E_J2S_ng_ep tmp832 = Taylor1(constant_term(q[6]) ^ float(constant_term(2)), order) v2[N] = Taylor1(constant_term(tmp830) + constant_term(tmp832), order) temp_004 = Array{Taylor1{_S}}(undef, size(newtonian1b_Potential)) - temp_004 .= Taylor1(zero(_S), order) + temp_004 .= Taylor1(zero(constant_term(q[1])), order) tmp835 = Array{Taylor1{_S}}(undef, size(μ)) - tmp835 .= Taylor1(zero(_S), order) + tmp835 .= Taylor1(zero(constant_term(q[1])), order) tmp837 = Array{Taylor1{_S}}(undef, size(μ)) - tmp837 .= Taylor1(zero(_S), order) + tmp837 .= Taylor1(zero(constant_term(q[1])), order) tmp839 = Array{Taylor1{_S}}(undef, size(μ)) - tmp839 .= Taylor1(zero(_S), order) + tmp839 .= Taylor1(zero(constant_term(q[1])), order) for i = 1:Nm1 temp_004[i] = Taylor1(constant_term(newtonian1b_Potential[i]) + constant_term(newtonianNb_Potential[N]), order) newtonianNb_Potential[N] = Taylor1(identity(constant_term(temp_004[i])), order) @@ -385,35 +385,35 @@ function TaylorIntegration._allocate_jetcoeffs!(::Val{RNp1BP_pN_A_J23E_J2S_ng_ep end _4ϕj[N] = Taylor1(constant_term(4) * constant_term(newtonianNb_Potential[N]), order) tmp845 = Array{Taylor1{_S}}(undef, size(v2)) - tmp845 .= Taylor1(zero(_S), order) + tmp845 .= Taylor1(zero(constant_term(q[1])), order) tmp847 = Array{Taylor1{_S}}(undef, size(vi_dot_vj)) - tmp847 .= Taylor1(zero(_S), order) + tmp847 .= Taylor1(zero(constant_term(q[1])), order) tmp848 = Array{Taylor1{_S}}(undef, size(tmp845)) - tmp848 .= Taylor1(zero(_S), order) + tmp848 .= Taylor1(zero(constant_term(q[1])), order) Xij_t_Ui = Array{Taylor1{_S}}(undef, size(X)) - Xij_t_Ui .= Taylor1(zero(_S), order) + Xij_t_Ui .= Taylor1(zero(constant_term(q[1])), order) Yij_t_Vi = Array{Taylor1{_S}}(undef, size(Y)) - Yij_t_Vi .= Taylor1(zero(_S), order) + Yij_t_Vi .= Taylor1(zero(constant_term(q[1])), order) Zij_t_Wi = Array{Taylor1{_S}}(undef, size(Z)) - Zij_t_Wi .= Taylor1(zero(_S), order) + Zij_t_Wi .= Taylor1(zero(constant_term(q[1])), order) tmp854 = Array{Taylor1{_S}}(undef, size(Xij_t_Ui)) - tmp854 .= Taylor1(zero(_S), order) + tmp854 .= Taylor1(zero(constant_term(q[1])), order) Rij_dot_Vi = Array{Taylor1{_S}}(undef, size(tmp854)) - Rij_dot_Vi .= Taylor1(zero(_S), order) + Rij_dot_Vi .= Taylor1(zero(constant_term(q[1])), order) tmp857 = Array{Taylor1{_S}}(undef, size(Rij_dot_Vi)) - tmp857 .= Taylor1(zero(_S), order) + tmp857 .= Taylor1(zero(constant_term(q[1])), order) pn1t7 = Array{Taylor1{_S}}(undef, size(tmp857)) - pn1t7 .= Taylor1(zero(_S), order) + pn1t7 .= Taylor1(zero(constant_term(q[1])), order) tmp860 = Array{Taylor1{_S}}(undef, size(pn1t7)) - tmp860 .= Taylor1(zero(_S), order) + tmp860 .= Taylor1(zero(constant_term(q[1])), order) pn1t2_7 = Array{Taylor1{_S}}(undef, size(ϕs_and_vs)) - pn1t2_7 .= Taylor1(zero(_S), order) + pn1t2_7 .= Taylor1(zero(constant_term(q[1])), order) tmp867 = Array{Taylor1{_S}}(undef, size(pNX_t_X)) - tmp867 .= Taylor1(zero(_S), order) + tmp867 .= Taylor1(zero(constant_term(q[1])), order) tmp868 = Array{Taylor1{_S}}(undef, size(tmp867)) - tmp868 .= Taylor1(zero(_S), order) + tmp868 .= Taylor1(zero(constant_term(q[1])), order) tmp869 = Array{Taylor1{_S}}(undef, size(tmp868)) - tmp869 .= Taylor1(zero(_S), order) + tmp869 .= Taylor1(zero(constant_term(q[1])), order) #= C:\Users\luisi\.julia\dev\NEOs.jl\cosa.jl:447 =# Threads.@threads for i = 1:10 ϕi_plus_4ϕj[i] = Taylor1(constant_term(newtonianNb_Potential_t[i]) + constant_term(_4ϕj[N]), order) tmp845[i] = Taylor1(constant_term(2) * constant_term(v2[i]), order) @@ -446,23 +446,23 @@ function TaylorIntegration._allocate_jetcoeffs!(::Val{RNp1BP_pN_A_J23E_J2S_ng_ep pNZ_t_pn3[i] = Taylor1(constant_term(acceph_t[3i]) * constant_term(pn3[i]), order) end tmp877 = Array{Taylor1{_S}}(undef, size(U_t_pn2)) - tmp877 .= Taylor1(zero(_S), order) + tmp877 .= Taylor1(zero(constant_term(q[1])), order) termpnx = Array{Taylor1{_S}}(undef, size(X_t_pn1)) - termpnx .= Taylor1(zero(_S), order) + termpnx .= Taylor1(zero(constant_term(q[1])), order) sumpnx = Array{Taylor1{_S}}(undef, size(termpnx)) - sumpnx .= Taylor1(zero(_S), order) + sumpnx .= Taylor1(zero(constant_term(q[1])), order) tmp880 = Array{Taylor1{_S}}(undef, size(V_t_pn2)) - tmp880 .= Taylor1(zero(_S), order) + tmp880 .= Taylor1(zero(constant_term(q[1])), order) termpny = Array{Taylor1{_S}}(undef, size(Y_t_pn1)) - termpny .= Taylor1(zero(_S), order) + termpny .= Taylor1(zero(constant_term(q[1])), order) sumpny = Array{Taylor1{_S}}(undef, size(termpny)) - sumpny .= Taylor1(zero(_S), order) + sumpny .= Taylor1(zero(constant_term(q[1])), order) tmp883 = Array{Taylor1{_S}}(undef, size(W_t_pn2)) - tmp883 .= Taylor1(zero(_S), order) + tmp883 .= Taylor1(zero(constant_term(q[1])), order) termpnz = Array{Taylor1{_S}}(undef, size(Z_t_pn1)) - termpnz .= Taylor1(zero(_S), order) + termpnz .= Taylor1(zero(constant_term(q[1])), order) sumpnz = Array{Taylor1{_S}}(undef, size(termpnz)) - sumpnz .= Taylor1(zero(_S), order) + sumpnz .= Taylor1(zero(constant_term(q[1])), order) for i = 1:10 tmp877[i] = Taylor1(constant_term(U_t_pn2[i]) + constant_term(pNX_t_pn3[i]), order) termpnx[i] = Taylor1(constant_term(X_t_pn1[i]) + constant_term(tmp877[i]), order) @@ -1758,87 +1758,87 @@ function TaylorIntegration._allocate_jetcoeffs!(::Val{RNp1BP_pN_A_J23E_J2S_eph_t dq[3] = Taylor1(identity(constant_term(q[6])), order) newtonianNb_Potential[N] = Taylor1(identity(constant_term(zero_q_1)), order) tmp1317 = Array{Taylor1{_S}}(undef, size(dq)) - tmp1317 .= Taylor1(zero(_S), order) + tmp1317 .= Taylor1(zero(constant_term(q[1])), order) tmp1319 = Array{Taylor1{_S}}(undef, size(ui)) - tmp1319 .= Taylor1(zero(_S), order) + tmp1319 .= Taylor1(zero(constant_term(q[1])), order) tmp1322 = Array{Taylor1{_S}}(undef, size(dq)) - tmp1322 .= Taylor1(zero(_S), order) + tmp1322 .= Taylor1(zero(constant_term(q[1])), order) tmp1324 = Array{Taylor1{_S}}(undef, size(vi)) - tmp1324 .= Taylor1(zero(_S), order) + tmp1324 .= Taylor1(zero(constant_term(q[1])), order) tmp1327 = Array{Taylor1{_S}}(undef, size(dq)) - tmp1327 .= Taylor1(zero(_S), order) + tmp1327 .= Taylor1(zero(constant_term(q[1])), order) tmp1329 = Array{Taylor1{_S}}(undef, size(wi)) - tmp1329 .= Taylor1(zero(_S), order) + tmp1329 .= Taylor1(zero(constant_term(q[1])), order) pn2x = Array{Taylor1{_S}}(undef, size(X)) - pn2x .= Taylor1(zero(_S), order) + pn2x .= Taylor1(zero(constant_term(q[1])), order) pn2y = Array{Taylor1{_S}}(undef, size(Y)) - pn2y .= Taylor1(zero(_S), order) + pn2y .= Taylor1(zero(constant_term(q[1])), order) pn2z = Array{Taylor1{_S}}(undef, size(Z)) - pn2z .= Taylor1(zero(_S), order) + pn2z .= Taylor1(zero(constant_term(q[1])), order) tmp1337 = Array{Taylor1{_S}}(undef, size(UU)) - tmp1337 .= Taylor1(zero(_S), order) + tmp1337 .= Taylor1(zero(constant_term(q[1])), order) tmp1340 = Array{Taylor1{_S}}(undef, size(X)) - tmp1340 .= Taylor1(zero(_S), order) + tmp1340 .= Taylor1(zero(constant_term(q[1])), order) tmp1342 = Array{Taylor1{_S}}(undef, size(Y)) - tmp1342 .= Taylor1(zero(_S), order) + tmp1342 .= Taylor1(zero(constant_term(q[1])), order) tmp1343 = Array{Taylor1{_S}}(undef, size(tmp1340)) - tmp1343 .= Taylor1(zero(_S), order) + tmp1343 .= Taylor1(zero(constant_term(q[1])), order) tmp1345 = Array{Taylor1{_S}}(undef, size(Z)) - tmp1345 .= Taylor1(zero(_S), order) + tmp1345 .= Taylor1(zero(constant_term(q[1])), order) tmp1353 = Array{Taylor1{_S}}(undef, size(pn2x)) - tmp1353 .= Taylor1(zero(_S), order) + tmp1353 .= Taylor1(zero(constant_term(q[1])), order) tmp1354 = Array{Taylor1{_S}}(undef, size(tmp1353)) - tmp1354 .= Taylor1(zero(_S), order) + tmp1354 .= Taylor1(zero(constant_term(q[1])), order) tmp1449 = Array{Taylor1{_S}}(undef, size(ui)) - tmp1449 .= Taylor1(zero(_S), order) + tmp1449 .= Taylor1(zero(constant_term(q[1])), order) tmp1451 = Array{Taylor1{_S}}(undef, size(vi)) - tmp1451 .= Taylor1(zero(_S), order) + tmp1451 .= Taylor1(zero(constant_term(q[1])), order) tmp1452 = Array{Taylor1{_S}}(undef, size(tmp1449)) - tmp1452 .= Taylor1(zero(_S), order) + tmp1452 .= Taylor1(zero(constant_term(q[1])), order) tmp1454 = Array{Taylor1{_S}}(undef, size(wi)) - tmp1454 .= Taylor1(zero(_S), order) + tmp1454 .= Taylor1(zero(constant_term(q[1])), order) tmp1365 = Array{Taylor1{_S}}(undef, size(X)) - tmp1365 .= Taylor1(zero(_S), order) + tmp1365 .= Taylor1(zero(constant_term(q[1])), order) tmp1367 = Array{Taylor1{_S}}(undef, size(Y)) - tmp1367 .= Taylor1(zero(_S), order) + tmp1367 .= Taylor1(zero(constant_term(q[1])), order) tmp1369 = Array{Taylor1{_S}}(undef, size(Z)) - tmp1369 .= Taylor1(zero(_S), order) + tmp1369 .= Taylor1(zero(constant_term(q[1])), order) tmp1371 = Array{Taylor1{_S}}(undef, size(t31)) - tmp1371 .= Taylor1(zero(_S), order) + tmp1371 .= Taylor1(zero(constant_term(q[1])), order) tmp1528 = Array{Taylor1{_S}}(undef, size(sin_ϕ)) - tmp1528 .= Taylor1(zero(_S), order) + tmp1528 .= Taylor1(zero(constant_term(q[1])), order) tmp1529 = Array{Taylor1{_S}}(undef, size(ϕ)) - tmp1529 .= Taylor1(zero(_S), order) + tmp1529 .= Taylor1(zero(constant_term(q[1])), order) tmp1381 = Array{Taylor1{_S}}(undef, size(sin2_ϕ)) - tmp1381 .= Taylor1(zero(_S), order) + tmp1381 .= Taylor1(zero(constant_term(q[1])), order) tmp1387 = Array{Taylor1{_S}}(undef, size(sin_ϕ)) - tmp1387 .= Taylor1(zero(_S), order) + tmp1387 .= Taylor1(zero(constant_term(q[1])), order) tmp1389 = Array{Taylor1{_S}}(undef, size(sin3_ϕ)) - tmp1389 .= Taylor1(zero(_S), order) + tmp1389 .= Taylor1(zero(constant_term(q[1])), order) tmp1393 = Array{Taylor1{_S}}(undef, size(sin2_ϕ)) - tmp1393 .= Taylor1(zero(_S), order) + tmp1393 .= Taylor1(zero(constant_term(q[1])), order) tmp1396 = Array{Taylor1{_S}}(undef, size(r_p2)) - tmp1396 .= Taylor1(zero(_S), order) + tmp1396 .= Taylor1(zero(constant_term(q[1])), order) tmp1397 = Array{Taylor1{_S}}(undef, size(Λ2)) - tmp1397 .= Taylor1(zero(_S), order) + tmp1397 .= Taylor1(zero(constant_term(q[1])), order) tmp1400 = Array{Taylor1{_S}}(undef, size(r_p1d2)) - tmp1400 .= Taylor1(zero(_S), order) + tmp1400 .= Taylor1(zero(constant_term(q[1])), order) tmp1401 = Array{Taylor1{_S}}(undef, size(Λ3)) - tmp1401 .= Taylor1(zero(_S), order) + tmp1401 .= Taylor1(zero(constant_term(q[1])), order) tmp1403 = Array{Taylor1{_S}}(undef, size(cos_ϕ)) - tmp1403 .= Taylor1(zero(_S), order) + tmp1403 .= Taylor1(zero(constant_term(q[1])), order) tmp1405 = Array{Taylor1{_S}}(undef, size(cos_ϕ)) - tmp1405 .= Taylor1(zero(_S), order) + tmp1405 .= Taylor1(zero(constant_term(q[1])), order) tmp1408 = Array{Taylor1{_S}}(undef, size(Λ2j_div_r4)) - tmp1408 .= Taylor1(zero(_S), order) + tmp1408 .= Taylor1(zero(constant_term(q[1])), order) tmp1412 = Array{Taylor1{_S}}(undef, size(Λ3j_div_r5)) - tmp1412 .= Taylor1(zero(_S), order) + tmp1412 .= Taylor1(zero(constant_term(q[1])), order) tmp1415 = Array{Taylor1{_S}}(undef, size(X)) - tmp1415 .= Taylor1(zero(_S), order) + tmp1415 .= Taylor1(zero(constant_term(q[1])), order) tmp1417 = Array{Taylor1{_S}}(undef, size(Y)) - tmp1417 .= Taylor1(zero(_S), order) + tmp1417 .= Taylor1(zero(constant_term(q[1])), order) tmp1419 = Array{Taylor1{_S}}(undef, size(Z)) - tmp1419 .= Taylor1(zero(_S), order) + tmp1419 .= Taylor1(zero(constant_term(q[1])), order) #= C:\Users\luisi\.julia\dev\NEOs.jl\cosa.jl:245 =# Threads.@threads for i = 1:Nm1 ui[i] = Taylor1(identity(constant_term(ss16asteph_t[3 * ((N - 1) + i) - 2])), order) vi[i] = Taylor1(identity(constant_term(ss16asteph_t[3 * ((N - 1) + i) - 1])), order) @@ -1974,13 +1974,13 @@ function TaylorIntegration._allocate_jetcoeffs!(::Val{RNp1BP_pN_A_J23E_J2S_eph_t tmp1462 = Taylor1(constant_term(q[6]) ^ float(constant_term(2)), order) v2[N] = Taylor1(constant_term(tmp1460) + constant_term(tmp1462), order) temp_004 = Array{Taylor1{_S}}(undef, size(newtonian1b_Potential)) - temp_004 .= Taylor1(zero(_S), order) + temp_004 .= Taylor1(zero(constant_term(q[1])), order) tmp1465 = Array{Taylor1{_S}}(undef, size(μ)) - tmp1465 .= Taylor1(zero(_S), order) + tmp1465 .= Taylor1(zero(constant_term(q[1])), order) tmp1467 = Array{Taylor1{_S}}(undef, size(μ)) - tmp1467 .= Taylor1(zero(_S), order) + tmp1467 .= Taylor1(zero(constant_term(q[1])), order) tmp1469 = Array{Taylor1{_S}}(undef, size(μ)) - tmp1469 .= Taylor1(zero(_S), order) + tmp1469 .= Taylor1(zero(constant_term(q[1])), order) for i = 1:Nm1 temp_004[i] = Taylor1(constant_term(newtonian1b_Potential[i]) + constant_term(newtonianNb_Potential[N]), order) newtonianNb_Potential[N] = Taylor1(identity(constant_term(temp_004[i])), order) @@ -1998,35 +1998,35 @@ function TaylorIntegration._allocate_jetcoeffs!(::Val{RNp1BP_pN_A_J23E_J2S_eph_t end _4ϕj[N] = Taylor1(constant_term(4) * constant_term(newtonianNb_Potential[N]), order) tmp1475 = Array{Taylor1{_S}}(undef, size(v2)) - tmp1475 .= Taylor1(zero(_S), order) + tmp1475 .= Taylor1(zero(constant_term(q[1])), order) tmp1477 = Array{Taylor1{_S}}(undef, size(vi_dot_vj)) - tmp1477 .= Taylor1(zero(_S), order) + tmp1477 .= Taylor1(zero(constant_term(q[1])), order) tmp1478 = Array{Taylor1{_S}}(undef, size(tmp1475)) - tmp1478 .= Taylor1(zero(_S), order) + tmp1478 .= Taylor1(zero(constant_term(q[1])), order) Xij_t_Ui = Array{Taylor1{_S}}(undef, size(X)) - Xij_t_Ui .= Taylor1(zero(_S), order) + Xij_t_Ui .= Taylor1(zero(constant_term(q[1])), order) Yij_t_Vi = Array{Taylor1{_S}}(undef, size(Y)) - Yij_t_Vi .= Taylor1(zero(_S), order) + Yij_t_Vi .= Taylor1(zero(constant_term(q[1])), order) Zij_t_Wi = Array{Taylor1{_S}}(undef, size(Z)) - Zij_t_Wi .= Taylor1(zero(_S), order) + Zij_t_Wi .= Taylor1(zero(constant_term(q[1])), order) tmp1484 = Array{Taylor1{_S}}(undef, size(Xij_t_Ui)) - tmp1484 .= Taylor1(zero(_S), order) + tmp1484 .= Taylor1(zero(constant_term(q[1])), order) Rij_dot_Vi = Array{Taylor1{_S}}(undef, size(tmp1484)) - Rij_dot_Vi .= Taylor1(zero(_S), order) + Rij_dot_Vi .= Taylor1(zero(constant_term(q[1])), order) tmp1487 = Array{Taylor1{_S}}(undef, size(Rij_dot_Vi)) - tmp1487 .= Taylor1(zero(_S), order) + tmp1487 .= Taylor1(zero(constant_term(q[1])), order) pn1t7 = Array{Taylor1{_S}}(undef, size(tmp1487)) - pn1t7 .= Taylor1(zero(_S), order) + pn1t7 .= Taylor1(zero(constant_term(q[1])), order) tmp1490 = Array{Taylor1{_S}}(undef, size(pn1t7)) - tmp1490 .= Taylor1(zero(_S), order) + tmp1490 .= Taylor1(zero(constant_term(q[1])), order) pn1t2_7 = Array{Taylor1{_S}}(undef, size(ϕs_and_vs)) - pn1t2_7 .= Taylor1(zero(_S), order) + pn1t2_7 .= Taylor1(zero(constant_term(q[1])), order) tmp1497 = Array{Taylor1{_S}}(undef, size(pNX_t_X)) - tmp1497 .= Taylor1(zero(_S), order) + tmp1497 .= Taylor1(zero(constant_term(q[1])), order) tmp1498 = Array{Taylor1{_S}}(undef, size(tmp1497)) - tmp1498 .= Taylor1(zero(_S), order) + tmp1498 .= Taylor1(zero(constant_term(q[1])), order) tmp1499 = Array{Taylor1{_S}}(undef, size(tmp1498)) - tmp1499 .= Taylor1(zero(_S), order) + tmp1499 .= Taylor1(zero(constant_term(q[1])), order) #= C:\Users\luisi\.julia\dev\NEOs.jl\cosa.jl:447 =# Threads.@threads for i = 1:10 ϕi_plus_4ϕj[i] = Taylor1(constant_term(newtonianNb_Potential_t[i]) + constant_term(_4ϕj[N]), order) tmp1475[i] = Taylor1(constant_term(2) * constant_term(v2[i]), order) @@ -2059,23 +2059,23 @@ function TaylorIntegration._allocate_jetcoeffs!(::Val{RNp1BP_pN_A_J23E_J2S_eph_t pNZ_t_pn3[i] = Taylor1(constant_term(acceph_t[3i]) * constant_term(pn3[i]), order) end tmp1507 = Array{Taylor1{_S}}(undef, size(U_t_pn2)) - tmp1507 .= Taylor1(zero(_S), order) + tmp1507 .= Taylor1(zero(constant_term(q[1])), order) termpnx = Array{Taylor1{_S}}(undef, size(X_t_pn1)) - termpnx .= Taylor1(zero(_S), order) + termpnx .= Taylor1(zero(constant_term(q[1])), order) sumpnx = Array{Taylor1{_S}}(undef, size(termpnx)) - sumpnx .= Taylor1(zero(_S), order) + sumpnx .= Taylor1(zero(constant_term(q[1])), order) tmp1510 = Array{Taylor1{_S}}(undef, size(V_t_pn2)) - tmp1510 .= Taylor1(zero(_S), order) + tmp1510 .= Taylor1(zero(constant_term(q[1])), order) termpny = Array{Taylor1{_S}}(undef, size(Y_t_pn1)) - termpny .= Taylor1(zero(_S), order) + termpny .= Taylor1(zero(constant_term(q[1])), order) sumpny = Array{Taylor1{_S}}(undef, size(termpny)) - sumpny .= Taylor1(zero(_S), order) + sumpny .= Taylor1(zero(constant_term(q[1])), order) tmp1513 = Array{Taylor1{_S}}(undef, size(W_t_pn2)) - tmp1513 .= Taylor1(zero(_S), order) + tmp1513 .= Taylor1(zero(constant_term(q[1])), order) termpnz = Array{Taylor1{_S}}(undef, size(Z_t_pn1)) - termpnz .= Taylor1(zero(_S), order) + termpnz .= Taylor1(zero(constant_term(q[1])), order) sumpnz = Array{Taylor1{_S}}(undef, size(termpnz)) - sumpnz .= Taylor1(zero(_S), order) + sumpnz .= Taylor1(zero(constant_term(q[1])), order) for i = 1:10 tmp1507[i] = Taylor1(constant_term(U_t_pn2[i]) + constant_term(pNX_t_pn3[i]), order) termpnx[i] = Taylor1(constant_term(X_t_pn1[i]) + constant_term(tmp1507[i]), order)