✨ Support OrbitStateVector in mean element fitting

src/api/Propagators.jl

@@ -40,29 +40,57 @@ abstract type OrbitPropagator{Tepoch<:Number, T<:Number} end
 
 """
     fit_mean_elements(::Val{:propagator}, vjd::AbstractVector{Tjd}, vr_i::AbstractVector{Tv}, vv_i::AbstractVector{Tv}; kwargs...) where {Tjd<:Number, Tv<:AbstractVector} -> <Mean elements>
+    fit_mean_elements(::Val{:propagator}, vsv::OrbitStateVector{Tepoch, T}; kwargs...) where {Tepoch<:Number, T<:Number} -> <Mean elements>
 
-Fit a set of mean elements for the `propagator` using the osculating elements represented by
-a set of position vectors `vr_i` [m] and a set of velocity vectors `vv_i` [m / s]
-represented in an inertial reference frame at instants in the array `vjd` [Julian Day]. The
-keywords `kwargs` depends on the propagator type.
+Fit a set of mean elements for the `propagator` using the osculating state vector
+represented in an intertial reference frame. The state vector can be represented using a set
+of position vectors `vr_i` [m] and a set of velocity vectors `vv_i` [m / s] obtained at the
+instants in the array `vjd` [Julian Day], or an array of `OrbitStateVector` `vsv` [SI],
+containing the same information. The keywords `kwargs` depends on the propagator type.
 
 This function returns the set of mean elements used to initialize the `propagator`.
 """
 function fit_mean_elements end
 
+function fit_mean_elements(
+    prop::Val,
+    vsv::Vector{OrbitStateVector{Tepoch, T}};
+    kwargs...
+) where {Tepoch<:Number, T<:Number}
+    vjd  = map(x -> x.t, vsv)
+    vr_i = map(x -> x.r, vsv)
+    vv_i = map(x -> x.v, vsv)
+
+    return fit_mean_elements(prop, vjd, vr_i, vv_i; kwargs...)
+end
+
 """
     fit_mean_elements!(orbp::OrbitPropagator, vjd::AbstractVector{Tjd}, vr_i::AbstractVector{Tv}, vv_i::AbstractVector{Tv}; kwargs...) where {Tjd<:Number, Tv<:AbstractVector} -> <Mean elements>
+    fit_mean_elements!(orbp::OrbitPropagator, vsv::Vector{OrbitStateVector{Tepoch, T}}; kwargs...) where {Tepoch<:Number, T<:Number} -> <Mean elements>
 
-Fit a set of mean elements for the propagator `orbp` using the osculating elements
-represented by a set of position vectors `vr_i` [m] and a set of velocity vectors `vv_i` [m
-/ s] represented in an inertial reference frame at instants in the array `vjd` [Julian Day].
-The keywords `kwargs` depends on the propagator type.
+Fit a set of mean elements for the propagator `orbp` using the osculating state vector
+represented in an intertial reference frame. The state vector can be represented using a set
+of position vectors `vr_i` [m] and a set of velocity vectors `vv_i` [m / s] obtained at the
+instants in the array `vjd` [Julian Day], or an array of `OrbitStateVector` `vsv` [SI],
+containing the same information. The keywords `kwargs` depends on the propagator type.
 
 This function returns the set of mean elements used to initialize the `propagator` and also
 initializes `orbp` with the fitted mean elements.
 """
 function fit_mean_elements! end
 
+function fit_mean_elements!(
+    orbp::OrbitPropagator,
+    vsv::Vector{OrbitStateVector{Tepoch, T}};
+    kwargs...
+) where {Tepoch<:Number, T<:Number}
+    vjd  = map(x -> x.t, vsv)
+    vr_i = map(x -> x.r, vsv)
+    vv_i = map(x -> x.v, vsv)
+
+    return fit_mean_elements!(orbp, vjd, vr_i, vv_i; kwargs...)
+end
+
 """
     init(::Val{:propagator}, args...; kwargs...) -> OrbitPropagator
 

test/api.jl

@@ -955,6 +955,74 @@ end
     end
 end
 
+@testset "Fitting Mean Elements Using OrbitStateVector" verbose = true begin
+    for (T, j2c) in ((Float64, j2c_egm2008), (Float32, j2c_egm2008_f32))
+        @testset "$T" begin
+            jd₀ = date_to_jd(2023, 1, 1, 0, 0, 0)
+
+            orb_input = KeplerianElements(
+                jd₀,
+                T(8000e3),
+                T(0.015),
+                T(28.5) |> deg2rad,
+                T(100)  |> deg2rad,
+                T(200)  |> deg2rad,
+                T(45)   |> deg2rad
+            )
+
+            orbp_input = Propagators.init(Val(:J2osc), orb_input; j2c = j2c)
+
+            # Create the osculating state vectors that we will use to fit the mean elements
+            # and compare with the reference Keplerian elements.
+            vsv = Propagators.propagate!(orbp_input, 0:10:12_000, OrbitStateVector)
+
+            # == Fit the Orbit =============================================================
+
+            # Create a new, uninitialized propagator.
+            j2d = J2Propagator{Float64, T}()
+            j2d.j2c = j2c
+            j2oscd = J2OsculatingPropagator{Float64, T}()
+            j2oscd.j2d = j2d
+            orbp = OrbitPropagatorJ2Osculating(j2oscd)
+
+            orb, _ = redirect_stdout(devnull) do
+                Propagators.fit_mean_elements!(
+                    orbp,
+                    vsv;
+                    max_iterations = 10,
+                    mean_elements_epoch = jd₀
+                )
+            end
+
+            @test orb.t ≈ orb_input.t
+            @test orb.a ≈ orb_input.a
+            @test orb.e ≈ orb_input.e
+            @test orb.i ≈ orb_input.i
+            @test orb.Ω ≈ orb_input.Ω
+            @test orb.ω ≈ orb_input.ω
+            @test orb.f ≈ orb_input.f atol = 1e-5
+
+            # Fit the mean elements without explicitly initializing the propagator.
+            orb, _ = redirect_stdout(devnull) do
+                Propagators.fit_mean_elements(
+                    Val(:J2osc),
+                    vsv;
+                    max_iterations = 10,
+                    mean_elements_epoch = jd₀
+                )
+            end
+
+            @test orb.t ≈ orb_input.t
+            @test orb.a ≈ orb_input.a
+            @test orb.e ≈ orb_input.e
+            @test orb.i ≈ orb_input.i
+            @test orb.Ω ≈ orb_input.Ω
+            @test orb.ω ≈ orb_input.ω
+            @test orb.f ≈ orb_input.f atol = 1e-5
+        end
+    end
+end
+
 @testset "Show" verbose = true begin
     T = Float64