Modeling With Ephemeris

Update: August 2024

All of the ephemeris fetching and parsing functionality has been moved into a super-package, EphemerisSources.jl. Check it out! The thumbnail image for this post is the JuliaAstro logo.

Astrodynamical Models

Models are combinations of simplifying assumptions, and equations; an earlier post introduces the topic of dynamics. There are many c ommon approximations astrodynamicists use to write equations — models — that describe our world. A spacecraft orbiting a celestial body may be approximated as massless: this produces the Restricted Two Body Problem. Two gravitational bodies in a circular orbit about their common center of mass can results in the Circular Restricted Three Body Problem; for elliptical orbits about their center of mass, the Elliptical Restricted Three Body Problem. Many of these models are provided within the Julia ecosystem, including through AstrodynamicalModels.jl.

We can simulate each model in code by manually applying initial conditions, and parameter values. For example, to simulate the Earth-Moon system within the Circular Restricted Three Body Problem, one might look up the two bodies’ masses, and the average distance between them. This introduces error into your calculations! In the real world, the average distance between both bodies will differ from the true distance at any given time. This approximation is usually fine for rough calculations, but what about space mission planning which relies on pre-computing a trajectory at a specific time?

Ephemeris Models

The positions and velocities of solar system bodies are commonly referred to as ephemeris. NASA provides two sources of ephemeris data for free: Horizons, which provides CSV-formatted data for large celestial bodies, and SPICE Kernels, which provide polynomial models which closely track a body’s true position and velocity in the solar system.

Both of these ephemeris sources are available in Julia through several packages, including HorizonsAPI.jl, HorizonsEphemeris.jl, SPICE.jl, SPICEKernels.jl, and SPICEBodies.jl. When these packages are used in combination with one another, they provide a programatic interface to fetching a celestial body’s position and velocity at any instance in time. This capability allows us to generate astrodynamical models which are accurate to a specific epoch!

Ephemeris Data

To build models at specific instances in time – specific epochs — we’ll need to download ephemeris data. This post will use SPICE kernel data, which provides state interpolation out of the box. NASA provides Generic Kernels for downloading, but SPICEKernels.jl can download (and cache) these kernel files for you! The SPICE toolkit’s furnsh function loads each kernel into the local kernel pool: a space in memory where the SPICE toolkit looks for kernel data.

Code

using SPICE, SPICEKernels

furnsh(
    de440(),                    # position and velocity data for major solar system bodies
    latest_leapseconds_lsk(),   # astronomical timekeeping
    gm_de440(),                 # mass parameters for major solar system bodies
    pck00011(),                 # physical properties of major solar system bodies
)

Other SPICE toolkit functions are available to read ephemeris data from the kernel pool, i.e. spkezr, spkpos, and many more. Once again, a Julia package — SPICEBodies.jl — offers a more convenient interface. Kernel bodies can be constructed using the KernelBody type, which can be called like a function to produce the body’s state (position and velocity) at any given epoch. All values are provided in kilometers, and kilometers per second.

As an example, let’s use an epoch during NASA’s Artemis 1 mission: midnight at December 1st, 2022. To build a Circular Restricted Three Body (CR3BP) model at this epoch, we’ll need the distance between the Earth and its Moon, \(a\). As with any CR3BP model construction, we’ll also need the reduced mass \(\mu\), defined by Equation 1.

\[ \mu = \frac{\text{min}(\mu_1, \mu_2)}{\mu_1 + \mu_2} \tag{1}\]

Code

using LinearAlgebra
using Unitful, UnitfulAstro
using SPICEBodies, AstroTime

epoch = TTEpoch(2022, 12, 1, 0, 0, 0.0)
earth = KernelBody("earth")
moon = KernelBody("moon")

a = let R = SPICEBodies.pos(moon, epoch; wrt=earth)
    norm(R) * u"km"
end

μ₁ = gm(earth) * u"km^3/s^2"
μ₂ = gm(moon) * u"km^3/s^2"

@show a μ₁ μ₂;

a = 375501.1650094356 km
μ₁ = 398600.4355070226 km³ s⁻²
μ₂ = 4902.800118457549 km³ s⁻²

The mass parameters for Earth and its (our) Moon are constant. The model parameter which varies over time is the distance between the two bodies, \(a\). The code above finds \(a\) at our desired epoch. How much does this value differ from the average distance between the Earth and the Moon? As the code below shows, setting our model parameter \(a\) to a specific epoch gives us almost 10 km of accuracy!

Code

using Statistics: mean

aₘ = let
    timepoints = range(epoch - 50years, epoch + 50years, step=1days)

    mean(
        time -> norm(moon(time; wrt=earth)),
        timepoints
    ) * u"km"
end

Δa = abs(a - aₘ)

@show Δa;

Δa = 9503.033361179056 km

Model Building

We have all the necessary information to build two Circular Restricted Three Body Problem models: one with the average distance between Earth and its Moon, and one with the specific distance at our desired Artemis 1 epoch. We can simulate both models. By comparing the relative error between the two, we can determine how useful applying a specific epoch to our model is.

Code

using GeneralAstrodynamics

EarthMoonAverage = CR3BPParameters(μ₁, μ₂, aₘ);
EarthMoonEpoch = CR3BPParameters(μ₁, μ₂, a);

Model Comparisons

Three body dynamics are known to be highly sensitive, chaotic systems. With a length unit difference of almost 10 km, EarthMoonAverage and EarthMoonEpoch will produce different simulation results! Let’s test this hypothesis using halo orbits.

Halo orbits are periodic orbits about lagrange points in the Three Body Problem. These orbits require some effort to find; most trajectories within CR3BP dynamics will not produce closed, periodic orbits! Iterative solvers are required to converge on periodic solutions. Let’s find a periodic orbit in the EarthMoonAverage model. The iterative solver, provided by the halo function in GeneralAstrodynamics.jl, provides an initial condition and orbital period. The initial condition provided by halo is normalized by the modeled distance \(a\) between the two bodies, and the synodic period as calculated in Equation 2.

\[ T_S = 2\pi \sqrt{\frac{a^3}{\mu_1 + \mu_2}} \tag{2}\]

Please note that the code below will take a couple of minutes to run the first time you execute it! The code should finish within a few seconds on subsequent computations.

Code

using OrdinaryDiffEq

ic, T = let sys = EarthMoonAverage
    orbit, period = halo(sys, Az=5_000u"km", L=2)

    DU = lengthunit(sys)
    DT = timeunit(sys)

    state = CartesianState(
        orbit.state.r * DU, orbit.state.v * DU / DT
    )

    state, period
end

@show ic;
@show T;

ic = Cartesian State with eltype Float64

  x = 430956.90838446154 km
  y = 0.0 km
  z = 4587.346204072134 km
  ẋ = 0.0 km s⁻¹
  ẏ = 0.029169837693280888 km s⁻¹
  ż = -0.0 km s⁻¹

T = 3.413299065325121

The initial condition ic is periodic (with period T) in the CR3BP model with the parameters defined by EarthMoonAverage. As Figure 2 shows, these conditions are not periodic at our desired epoch in dimensioned units!

We need to take care of our units here. As previously stated, CR3BP dynamics are typically numerically integrated with states that are normalized by the parameters specific to the system. One benefit of this normalization is the normalized quantities are relevant across systems with the same mass parameters! The dimensioned quantities, however, are not. For illustrative purposes, we will consider the dimensioned quantities for the remainder of this post.

Code

orbit = (;
    average=Orbit(ic, EarthMoonAverage),
    epoch=Orbit(ic, EarthMoonEpoch)
);

Figure 1: Orbit Simulated in EarthMoonAverage

Figure 2: Orbit Simulated in EarthMoonEpoch

This result is important. A periodic orbit found by applying the average distance between the Earth and the Moon was not periodic at our desired mission epoch. This is why ephemeris data is important! Let’s use the ephemeris data to solve for an orbit which is periodic at our desired epoch.

Note that the code below is unnecessarily verbose. The same Orbit constructor is used throughout this post for clarity, but the variable orbit below should be identical to variable orbit_at_epoch.

Code

ic_at_epoch, period_at_epoch = let sys = EarthMoonEpoch
    orbit, period = halo(sys, Az=5_000u"km", L=2)

    DU = lengthunit(sys)
    DT = timeunit(sys)

    state = CartesianState(
        orbit.state.r * DU, orbit.state.v * DU / DT
    )

    state, period
end

orbit_at_epoch = Orbit(ic_at_epoch, EarthMoonEpoch)

Circular Restricted Three-body Orbit with eltype Float64

  Cartesian State

    x = 1.1193033959124965
    y = 0.0
    z = 0.012216225670010598
    ẋ = 0.0
    ẏ = 0.17918236349802066
    ż = -0.0

  CR3BP parameters

    μ = 0.012150584395829193

Figure 3: Periodic Orbit Simulated in EarthMoonEpoch

Conclusions

Ephemeris data is important, and necessary to construct models which accurately describe astrodynamical systems at specific epochs. As you may notice, this whole process can be automated through a package interface. I plan to integrate epoch-specific modeling, just as described here, within AstrodynamicalModels.jl!

I’m grateful for the free information we have available at all times. In short, if you want to know where the planets are, check out JPL Horizons and the Generic SPICE Kernels. If you want to apply such information to astrodynamical simulations, check out Julia!