⚠️ 🔧 Add orbit epoch to the propagators

The propagators structures now have a field called `epoch` that stores the Julian Day of the epoch in which the orbit was defined. This is important to the Deep Space algorithm of SGP4 that will be implemented. When calling `propagate!` or `step!`, the `Orbit` structure that is returned contains in the `t` parameter the epoch of those orbit elements in Julian Day. Notice that this is a **BREAKING CHANGE** if you were using `propagate!` with an orbit in which the epoch was not zero or if you were using the `orb.t` parameter.
JuliaSpace · May 30, 2018 · dea8cbf · dea8cbf
1 parent bceaacd
commit dea8cbf
Show file tree

Hide file tree

Showing 8 changed files with 70 additions and 59 deletions.
diff --git a/src/orbit/propagators/j2.jl b/src/orbit/propagators/j2.jl
@@ -95,14 +95,14 @@ j2_gc_wgs72 = J2_GravCte(
 ################################################################################
 
 """
-    function j2_init(j2_gc::J2_GravCte{T}, t_0::Number, n_0::Number, e_0::Number, i_0::Number, Ω_0::Number, ω_0::Number, M_0::Number, dn_o2::Number, ddn_o6::Number) where T
+    function j2_init(j2_gc::J2_GravCte{T}, epoch::Number, n_0::Number, e_0::Number, i_0::Number, Ω_0::Number, ω_0::Number, M_0::Number, dn_o2::Number, ddn_o6::Number) where T
 
 Initialize the data structure of J2 orbit propagator algorithm.
 
 # Args
 
 * `j2_gc`: J2 orbit propagator gravitational constants (see `J2_GravCte`).
-* `t_0`: Epoch of the orbital elements [s].
+* `epoch`: Epoch of the orbital elements [Julian Day].
 * `n_0`: Mean motion at epoch [rad/s].
 * `e_0`: Initial eccentricity.
 * `i_0`: Initial inclination [rad].
@@ -118,7 +118,7 @@ The structure `J2_Structure` with the initialized parameters.
 
 """
 function j2_init(j2_gc::J2_GravCte{T},
-                 t_0::Number,
+                 epoch::Number,
                  n_0::Number,
                  e_0::Number,
                  i_0::Number,
@@ -180,8 +180,8 @@ function j2_init(j2_gc::J2_GravCte{T},
     # Create the output structure with the data.
     J2_Structure{T}(
 
-        t_0, a_0, n_0, e_0, i_0, Ω_0, ω_0, M_0, dn_o2, ddn_o6, a_0, e_0, i_0,
-        Ω_0, ω_0, M_0, n_0, f_0, C1, C2, C3, C4, j2_gc
+        epoch, a_0, n_0, e_0, i_0, Ω_0, ω_0, M_0, 0, dn_o2, ddn_o6, a_0, e_0,
+        i_0, Ω_0, ω_0, M_0, n_0, f_0, C1, C2, C3, C4, j2_gc
 
     )
 end
@@ -216,7 +216,7 @@ function j2!(j2d::J2_Structure{T}, t::Number) where T
     @unpack_J2_GravCte   j2_gc
 
     # Time elapsed since epoch.
-    Δt = t - t_0
+    Δt = t
 
     # Propagate the orbital elements.
     a_k = a_0 - C1*Δt
@@ -234,6 +234,7 @@ function j2!(j2d::J2_Structure{T}, t::Number) where T
     (r_i_k, v_i_k) = kepler_to_rv(a_k*R0, e_k, i_k, Ω_k, ω_k, f_k)
 
     # Update the J2 orbit propagator structure.
+    j2d.Δt  = Δt
     j2d.a_k = a_k
     j2d.e_k = e_k
     j2d.i_k = i_k

diff --git a/src/orbit/propagators/j2_api.jl b/src/orbit/propagators/j2_api.jl
@@ -31,15 +31,15 @@ export init_orbit_propagator, step!, propagate!
 ################################################################################
 
 """
-    function init_orbit_propagator(::Type{Val{:J2}}, t_0::Number, n_0::Number, e_0::Number, i_0::Number, Ω_0::Number, ω_0::Number, M_0::Number, dn_o2::Number = 0, ddn_o6::Number = 0, j2_gc::J2_GravCte{T} = j2_gc_wgs84) where T
+    function init_orbit_propagator(::Type{Val{:J2}}, epoch::Number, n_0::Number, e_0::Number, i_0::Number, Ω_0::Number, ω_0::Number, M_0::Number, dn_o2::Number = 0, ddn_o6::Number = 0, j2_gc::J2_GravCte{T} = j2_gc_wgs84) where T
 
 Initialize the J2 orbit propagator using the initial orbit specified by the
-elements `t_0, `n_0, `e_0`, `i_0`, `Ω_0`, `ω_0`, and `M_0`, and the
+elements `epoch, `n_0, `e_0`, `i_0`, `Ω_0`, `ω_0`, and `M_0`, and the
 gravitational parameters `j2_gc` (see `J2_GravCte`).
 
 # Args
 
-* `t_0`: Initial orbit epoch [s].
+* `epoch`: Initial orbit epoch [Julian Day].
 * `n_0`: Initial angular velocity [rad/s].
 * `e_0`: Initial eccentricity.
 * `i_0`: Initial inclination [rad].
@@ -60,7 +60,7 @@ of the orbit propagator.
 
 """
 function init_orbit_propagator(::Type{Val{:J2}},
-                               t_0::Number,
+                               epoch::Number,
                                n_0::Number,
                                e_0::Number,
                                i_0::Number,
@@ -71,11 +71,11 @@ function init_orbit_propagator(::Type{Val{:J2}},
                                ddn_o6::Number = 0,
                                j2_gc::J2_GravCte{T} = j2_gc_wgs84) where T
     # Create the new Two Body propagator structure.
-    j2d = j2_init(j2_gc, t_0, n_0, e_0, i_0, Ω_0, ω_0, M_0, dn_o2, ddn_o6)
+    j2d = j2_init(j2_gc, epoch, n_0, e_0, i_0, Ω_0, ω_0, M_0, dn_o2, ddn_o6)
 
     # Create the `Orbit` structure.
     orb_0 =
-        Orbit{T,T,T,T,T,T,T}(t_0, j2d.a_0*j2_gc.R0, e_0, i_0, Ω_0, ω_0, j2d.f_k)
+        Orbit{T,T,T,T,T,T,T}(epoch, j2d.a_0*j2_gc.R0, e_0, i_0, Ω_0, ω_0, j2d.f_k)
 
     # Create and return the orbit propagator structure.
     OrbitPropagatorJ2(orb_0, j2d)
@@ -144,7 +144,7 @@ function init_orbit_propagator(::Type{Val{:J2}},
                                tle::TLE,
                                j2_gc::J2_GravCte = j2_gc_wgs84)
     init_orbit_propagator(Val{:J2},
-                          tle.epoch_day*24*60*60,
+                          tle.epoch,
                           tle.n*2*pi/(24*60*60),
                           tle.e,
                           tle.i*pi/180,
@@ -188,10 +188,10 @@ function step!(orbp::OrbitPropagatorJ2, Δt::Number)
     j2d = orbp.j2d
 
     # Propagate the orbit.
-    (r_i, v_i) = j2!(j2d, orb.t + Δt)
+    (r_i, v_i) = j2!(j2d, j2d.Δt + Δt)
 
     # Update the elements in the `orb` structure.
-    orb.t += Δt
+    orb.t += Δt/86400
     orb.a  = j2d.a_k*j2d.j2_gc.R0
     orb.e  = j2d.e_k
     orb.i  = j2d.i_k
@@ -245,10 +245,10 @@ function propagate!(orbp::OrbitPropagatorJ2, t::Vector)
 
     for k in t
         # Propagate the orbit.
-        (r_i_k, v_i_k) = j2!(j2d, j2d.t_0 + k)
+        (r_i_k, v_i_k) = j2!(j2d, k)
 
         # Update the elements in the `orb` structure.
-        orb.t = j2d.t_0 + k
+        orb.t = j2d.epoch + k/86400
         orb.a = j2d.a_k*j2d.j2_gc.R0
         orb.e = j2d.e_k
         orb.i = j2d.i_k

diff --git a/src/orbit/propagators/sgp4.jl b/src/orbit/propagators/sgp4.jl
@@ -85,14 +85,14 @@ sgp4_gc_wgs72 = SGP4_GravCte(
 ################################################################################
 
 """
-    function sgp4_init(spg4_gc::SGP4_GravCte{T}, t_0::Number, n_0::Number, e_0::Number, i_0::Number, Ω_0::Number, ω_0::Number, M_0::Number, bstar::Number) where T
+    function sgp4_init(spg4_gc::SGP4_GravCte{T}, epoch::Number, n_0::Number, e_0::Number, i_0::Number, Ω_0::Number, ω_0::Number, M_0::Number, bstar::Number) where T
 
 Initialize the data structure of SGP4 algorithm.
 
 # Args
 
 * `spg4_gc`: SPG4 gravitational constants (see `SGP4_GravCte`).
-* `t_0`: Epoch of the orbital elements [min].
+* `epoch`: Epoch of the orbital elements [Julian Day].
 * `n_0`: SGP type "mean" mean motion at epoch [rad/min].
 * `e_0`: "Mean" eccentricity at epoch.
 * `i_0`: "Mean" inclination at epoch [rad].
@@ -107,7 +107,7 @@ The structure `SGP4_Structure` with the initialized parameters.
 
 """
 function sgp4_init(sgp4_gc::SGP4_GravCte{T},
-                   t_0::Number,
+                   epoch::Number,
                    n_0::Number,
                    e_0::Number,
                    i_0::Number,
@@ -261,12 +261,12 @@ function sgp4_init(sgp4_gc::SGP4_GravCte{T},
     n_k = nll_0
 
     # Create the output structure with the data.
-    SGP4_Structure(
+    SGP4_Structure{T}(
 
-        t_0, n_0, e_0, i_0, Ω_0, ω_0, M_0, bstar, a_k, e_k, i_k, Ω_k, ω_k, M_k,
-        n_k, all_0, nll_0, AE, QOMS2T, β_0, ξ, η, sin_i_0, θ, θ2, θ3, θ4, A_30,
-        k_2, k_4, C1, C3, C4, C5, D2, D3, D4, dotM, dotω, dotΩ1, dotΩ, isimp,
-        sgp4_gc
+        epoch, n_0, e_0, i_0, Ω_0, ω_0, M_0, bstar, 0, a_k, e_k, i_k, Ω_k, ω_k,
+        M_k, n_k, all_0, nll_0, AE, QOMS2T, β_0, ξ, η, sin_i_0, θ, θ2, θ3, θ4,
+        A_30, k_2, k_4, C1, C3, C4, C5, D2, D3, D4, dotM, dotω, dotΩ1, dotΩ,
+        isimp, sgp4_gc
 
        )
 end
@@ -294,7 +294,7 @@ function sgp4!(sgp4d::SGP4_Structure{T}, t::Number) where T
     @unpack_SGP4_GravCte   sgp4_gc
 
     # Time elapsed since epoch.
-    Δt = t-t_0
+    Δt = t
 
     # Secular effects of atmospheric drag and gravitation.
     # ====================================================
@@ -483,6 +483,7 @@ function sgp4!(sgp4d::SGP4_Structure{T}, t::Number) where T
     v_TEME = (dot_r_k*Uv + r_dot_f_k*Vv)*R0/60
 
     # Update the variables.
+    sgp4d.Δt  = Δt
     sgp4d.a_k = a
     sgp4d.e_k = e
     sgp4d.i_k = i_k

diff --git a/src/orbit/propagators/sgp4_api.jl b/src/orbit/propagators/sgp4_api.jl
@@ -47,7 +47,7 @@ macro update_orb!(orbp, t)
         local sgp4d = $(esc(orbp)).sgp4d
         local R0    = $(esc(orbp)).sgp4_gc.R0
 
-        orb.t = $(esc(t))
+        orb.t = sgp4d.epoch + $(esc(t))/86400
         orb.a = sgp4d.a_k*R0*1000
         orb.e = sgp4d.e_k
         orb.i = sgp4d.i_k
@@ -62,17 +62,17 @@ end
 ################################################################################
 
 """
-    function init_orbit_propagator(::Type{Val{:sgp4}}, t_0::Number, n_0::Number, e_0::Number, i_0::Number, Ω_0::Number, ω_0::Number, M_0::Number, bstar::Number, sgp4_gc::SGP4_GravCte{T} = sgp4_gc_wgs84) where T
+    function init_orbit_propagator(::Type{Val{:sgp4}}, epoch::Number, n_0::Number, e_0::Number, i_0::Number, Ω_0::Number, ω_0::Number, M_0::Number, bstar::Number, sgp4_gc::SGP4_GravCte{T} = sgp4_gc_wgs84) where T
 
 Initialize the SGP4 orbit propagator using the initial orbit specified by the
-elements `t_0, `n_0, `e_0`, `i_0`, `Ω_0`, `ω_0`, and `M_0`, the B* parameter
+elements `epoch, `n_0, `e_0`, `i_0`, `Ω_0`, `ω_0`, and `M_0`, the B* parameter
 `bstar`, and the gravitational constants in the structure `sgp4_gc`.
 
 Notice that the orbit elements **must be** represented in TEME frame.
 
 # Args
 
-* `t_0`: Initial orbit epoch [s].
+* `epoch`: Initial orbit epoch [Julian Day].
 * `n_0`: Initial angular velocity [rad/s].
 * `e_0`: Initial eccentricity.
 * `i_0`: Initial inclination [rad]
@@ -89,7 +89,7 @@ information of the orbit propagator.
 
 """
 function init_orbit_propagator(::Type{Val{:sgp4}},
-                               t_0::Number,
+                               epoch::Number,
                                n_0::Number,
                                e_0::Number,
                                i_0::Number,
@@ -100,7 +100,7 @@ function init_orbit_propagator(::Type{Val{:sgp4}},
                                sgp4_gc::SGP4_GravCte{T} = sgp4_gc_wgs84) where T
     # Create the new SGP4 structure.
     sgp4d = sgp4_init(sgp4_gc,
-                      t_0/60.0,
+                      epoch,
                       n_0*60.0,
                       e_0,
                       i_0,
@@ -110,7 +110,7 @@ function init_orbit_propagator(::Type{Val{:sgp4}},
                       bstar)
 
     # Create the `Orbit` structure.
-    orb_0 = Orbit{T,T,T,T,T,T,T}(t_0,
+    orb_0 = Orbit{T,T,T,T,T,T,T}(epoch,
                                  sgp4d.a_k*sgp4_gc_wgs84.R0,
                                  e_0,
                                  i_0,
@@ -181,8 +181,9 @@ information of the orbit propagator.
 function init_orbit_propagator(::Type{Val{:sgp4}},
                                tle::TLE,
                                sgp4_gc::SGP4_GravCte = sgp4_gc_wgs84)
+
     init_orbit_propagator(Val{:sgp4},
-                          tle.epoch_day*24*60*60,
+                          tle.epoch,
                           tle.n*2*pi/(24*60*60),
                           tle.e,
                           tle.i*pi/180,
@@ -218,14 +219,14 @@ function step!(orbp::OrbitPropagatorSGP4{T}, Δt::Number) where T
     sgp4d = orbp.sgp4d
 
     # Propagate the orbit.
-    (r_teme, v_teme) = sgp4!(sgp4d, (orb.t + Δt)/60)
+    (r_teme, v_teme) = sgp4!(sgp4d, sgp4d.Δt + Δt/60)
 
     # Convert km to m.
     r_teme *= 1000
     v_teme *= 1000
 
     # Update the elements in the `orb` structure.
-    @update_orb!(orbp, orb.t + Δt)
+    @update_orb!(orbp, sgp4d.Δt*60)
 
     # Return the information about the step.
     (copy(orbp.orb), r_teme, v_teme)
@@ -265,14 +266,14 @@ function propagate!(orbp::OrbitPropagatorSGP4{T}, t::Vector) where T
 
     for k in t
         # Propagate the orbit.
-        (r_teme_k, v_teme_k) = sgp4!(sgp4d, sgp4d.t_0 + k/60)
+        (r_teme_k, v_teme_k) = sgp4!(sgp4d, k/60)
 
         # Convert km to m.
         r_teme_k *= 1000
         v_teme_k *= 1000
 
         # Update the elements in the `orb` structure.
-        @update_orb!(orbp, sgp4d.t_0*60 + k)
+        @update_orb!(orbp, k)
 
         push!(result_orb, copy(orb))
         push!(result_r,   r_teme_k)

diff --git a/src/orbit/propagators/twobody.jl b/src/orbit/propagators/twobody.jl
@@ -44,13 +44,14 @@ export twobody_init, twobody!
 
 # Copy for TwoBody_Structure.
 function Base.copy(m::TwoBody_Structure)
-    TwoBody_Structure(m.t_0,
+    TwoBody_Structure(m.epoch,
                       m.n_0,
                       m.e_0,
                       m.i_0,
                       m.Ω_0,
                       m.ω_0,
                       m.M_0,
+                      m.Δt,
                       m.a,
                       m.M_k,
                       m.f_k,
@@ -62,13 +63,13 @@ end
 ################################################################################
 
 """
-    function twobody_init(t_0::Number, n_0::Number, e_0::Number, i_0::Number, Ω_0::Number, ω_0::Number, M_0::Number, μ::T) where T
+    function twobody_init(epoch::Number, n_0::Number, e_0::Number, i_0::Number, Ω_0::Number, ω_0::Number, M_0::Number, μ::T) where T
 
 Initialize the data structure of Two Body orbit propagator algorithm.
 
 # Args
 
-* `t_0`: Epoch of the orbital elements [s].
+* `epoch`: Epoch of the orbital elements [s].
 * `n_0`: Mean motion at epoch [rad/s].
 * `e_0`: Initial eccentricity.
 * `i_0`: Initial inclination [rad].
@@ -82,7 +83,7 @@ Initialize the data structure of Two Body orbit propagator algorithm.
 The structure `TwoBody_Structure` with the initialized parameters.
 
 """
-function twobody_init(t_0::Number,
+function twobody_init(epoch::Number,
                       n_0::Number,
                       e_0::Number,
                       i_0::Number,
@@ -102,7 +103,7 @@ function twobody_init(t_0::Number,
     f = M_to_f(e_0, M_0)
 
     # Create and return the Two Body orbit propagator structure.
-    TwoBody_Structure{T}(t_0, n_0, e_0, i_0, Ω_0, ω_0, M_0, a, M_0, f, μ)
+    TwoBody_Structure{T}(epoch, n_0, e_0, i_0, Ω_0, ω_0, M_0, 0, a, M_0, f, μ)
 end
 
 """
@@ -130,7 +131,8 @@ from a TLE, then the inertial frame will be TEME.
 """
 function twobody!(tbd::TwoBody_Structure, t::Number)
     # Time elapsed since epoch.
-    Δt = t - tbd.t_0
+    Δt     = t
+    tbd.Δt = Δt
 
     # Update the mean anomaly.
     tbd.M_k = tbd.M_0 + tbd.n_0*Δt