Overview

There are a number of equations of motion implemented in HighFidelityEphemerisModel.jl.

• functions starting with eom_ integrates the translational state ([x,y,z,vx,vy,vz])
• functions starting with eom_stm_ integrates both the translational state and the flattened 6-by-6 STM.

x_stm_tf = sol.u[end]                       # concatenated state & STM (flattened)
x_tf     = x_stm_tf[1:6]                    # final state [x,y,z,vx,vy,vz]
STM_tf   = reshape(sol.u[end][7:42],6,6)'   # final 6-by-6 STM

Dynamics model

In HighFidelityEphemerisModel.jl, the dynamcis consists of the central gravitational term, together with the following perturbations:

• third-body perturbations
• spherical harmonics
• solar radiation pressure (todo)
• drag (todo)

\[\dot{\boldsymbol{x}}(t) = \begin{bmatrix} \dot{\boldsymbol{r}}(t) \\ \dot{\boldsymbol{v}}(t) \end{bmatrix} = \begin{bmatrix} \boldsymbol{v}(t) \\ -\dfrac{\mu}{\| \boldsymbol{r}(t) \|_2^3}\boldsymbol{r}(t) \end{bmatrix} + \sum_{i} \begin{bmatrix} \boldsymbol{0}_{3 \times 1} \\ \boldsymbol{a}_{\mathrm{3bd},i}(t) \end{bmatrix} + \begin{bmatrix} \boldsymbol{0}_{3 \times 1} \\ \boldsymbol{a}_{\mathrm{SH},n_{\max}}(t) \end{bmatrix} + \begin{bmatrix} \boldsymbol{0}_{3 \times 1} \\ \boldsymbol{a}_{\mathrm{SRP}}(t) \end{bmatrix}\]

Third-body perturbation

The third-body perturbation due to body $i$, $\boldsymbol{a}_{\mathrm{3bd},i}$, is given by

\[\boldsymbol{a}_{\mathrm{3bd},i}(t) = -\mu_i \left( \dfrac{\boldsymbol{r}(t) - \boldsymbol{r}_i(t)}{\| \boldsymbol{r}(t) - \boldsymbol{r}_i(t) \|_2^3} + \dfrac{\boldsymbol{r}_i(t)}{\| \boldsymbol{r}_i(t) \|_2^3} \right)\]

where $\boldsymbol{r}_i$ is the position vector of the perturbing body. In HEFM.jl, this term is implemented using Battin's $F(q)$ function:

\[\boldsymbol{a}_{\mathrm{3bd},i}(t) = -\dfrac{\mu_i}{\| \boldsymbol{r}(t) - \boldsymbol{r}_i(t) \|_2^3} (\boldsymbol{r}(t) + F(q_i)\boldsymbol{r}_i(t))\]

where $F(q_i)$ is given by

\[F(q_i) = q_i \left( \dfrac{3 + 3q_i + q_i^2}{1 + (\sqrt{1 + q_i})^3} \right) ,\quad q_i = \dfrac{\boldsymbol{r}(t)^T (\boldsymbol{r}(t) - 2\boldsymbol{r}_i(t))}{\boldsymbol{r}_i(t)^T \boldsymbol{r}_i(t)}\]

Spherical Harmonics

The spherical harmonics perturbation $\boldsymbol{a}_{\mathrm{SH},n_{\max}}$ is given by

\[\boldsymbol{a}_{\mathrm{SH},n_{\max}} = \sum_{n=2}^{n_{\max}} \sum_{m=0}^n \boldsymbol{a}_{\mathrm{SH},nm}\]

where $\boldsymbol{a}_{\mathrm{SH},nm}$ is given by

\[\boldsymbol{a}_{\mathrm{SH},nm} = \begin{bmatrix} \ddot{x}_{n m} \\ \ddot{y}_{n m} \\ \ddot{z}_{n m} \end{bmatrix}\]

where

\[\begin{aligned} & \ddot{x}_{n m} = \begin{cases} \frac{G M}{R_{\oplus}^2} \cdot\left\{-C_{n 0} V_{n+1,1}\right\} & m = 0 \\[1.0em] \frac{G M}{R_{\oplus}^2} \cdot \frac{1}{2} \cdot\left\{\left(-C_{n m} V_{n+1, m+1}-S_{n m} W_{n+1, m+1}\right) + \frac{(n-m+2)!}{(n-m)!} \cdot\left(+C_{n m} V_{n+1, m-1}+S_{n m} W_{n+1, m-1}\right)\right\} & m > 0 \end{cases} \\[2.5em] & \ddot{y}_{n m} = \begin{cases} \frac{G M}{R_{\oplus}^2} \cdot\left\{-C_{n 0} W_{n+1,1}\right\} & m = 0 \\[1.0em] \frac{G M}{R_{\oplus}^2} \cdot \frac{1}{2} \cdot\left\{\left(-C_{n m} \cdot W_{n+1, m+1}+S_{n m} \cdot V_{n+1, m+1}\right) + \frac{(n-m+2)!}{(n-m)!} \cdot\left(-C_{n m} W_{n+1, m-1}+S_{n m} V_{n+1, m-1}\right)\right\} & m > 0 \end{cases} \\[2.5em] & \ddot{z}_{n m} = \frac{G M}{R_{\oplus}^2} \cdot\left\{(n-m+1) \cdot\left(-C_{n m} V_{n+1, m}-S_{n m} W_{n+1, m}\right)\right\} \end{aligned}\]

and

\[V_{n m}=\left(\frac{R_{\oplus}}{r}\right)^{n+1} \cdot P_{n m}(\sin \phi) \cdot \cos m \lambda , \quad W_{n m}=\left(\frac{R_{\oplus}}{r}\right)^{n+1} \cdot P_{n m}(\sin \phi) \cdot \sin m \lambda\]

(c.f. Montenbruck & Gill Chapter 3.2)

Solar Radiation Pressure

Cannonball model

TODO

List of equations of motion in HighFidelityEphemerisModel.jl

The table below summarizes the equations of motion. Note:

• Nbody: central gravity term + third-body perturbations ($\boldsymbol{a}_{\mathrm{3bd},i}$)
• NbodySH: central gravity term + third-body perturbations + spherical harmonics perturbations up to nmax degree ($\boldsymbol{a}_{\mathrm{SH},n_{\max}}$)
• The STM is integrated with the Jacobian, which is computed either analytically (using symbolic derivative) or via ForwardDiff (functions containing _fd)

Initializing the parameter

We first need to define the parameter struct to be parsed as argument to the equations of motion.

Below is the most general example compatible with eom_NbodySH_Interp!/eom_stm_NbodySH_Interp_fd!:

using OrdinaryDiffEq
using HighFidelityEphemerisModel

# load SPICE kernels
spice_dir = ENV["SPICE"]
furnsh(joinpath(spice_dir, "lsk", "naif0012.tls"))
furnsh(joinpath(spice_dir, "spk", "de440.bsp"))
furnsh(joinpath(spice_dir, "pck", "gm_de440.tpc"))
furnsh(joinpath(spice_dir, "pck", "moon_pa_de440_200625.bpc"))
furnsh(joinpath(spice_dir, "fk", "moon_de440_250416.tf"))

naif_ids = ["301", "399", "10"]        # NAIF IDs of bodies to be included; first ID is of the central body
GMs = [bodvrd(ID, "GM", 1)[1] for ID in naif_ids]      # in km^3/s^2
naif_frame = "J2000"
abcorr = "NONE"
DU = 1e5                               # canonical distance unit, in km

nmax = 4                               # using up to 4-by-4 spherical harmonics
filepath_spherical_harmonics = "HighFidelityEphemerisModel.jl/data/luna/gggrx_1200l_sha_20x20.tab"

et0 = str2et("2026-01-05T00:00:00")    # reference epoch
etf = et0 + 30 * 86400.0
interpolate_ephem_span = [et0, etf]    # range of epoch to interpolate ephemeris
interpolation_time_step = 1000.0       # time-step to sample ephemeris for interpolation

parameters = HighFidelityEphemerisModel.HighFidelityEphemerisModelParameters(
    et0, DU, GMs, naif_ids, naif_frame, abcorr;
    interpolate_ephem_span=interpolate_ephem_span,
    interpolation_time_step = interpolation_time_step,
    filepath_spherical_harmonics = filepath_spherical_harmonics,
    nmax = nmax,
    frame_PCPF = "MOON_PA",
)

Note:

• NAIF body IDs are defined according to: https://naif.jpl.nasa.gov/pub/naif/toolkitdocs/C/req/naifids.html
• if using _SPICE equations of motion, you do not need to parse interpolate_ephem_span and interpolation_time_step
• if using Nbody dynamics instead of NbodySH, you do not need to parse filepath_spherical_harmonics, nmax, and frame_PCPF

Solving an Initial Value Problem

The integration is done with the OrdinaryDiffEq.jl library (or equivalently with DifferentialEquations.jl).

# initial state (in canonical scale)
x0 = [1.05, 0.0, 0.3, 0.5, 1.0, 0.0]

# time span (in canonical scale)
tspan = (0.0, 6 * 3600/parameters.TU)

# solve with SPICE
prob_spice = ODEProblem(HighFidelityEphemerisModel.eom_NbodySH_SPICE!, x0, tspan, parameters)
sol_spice = solve(prob_spice, Vern8(), reltol=1e-14, abstol=1e-14)

# solve with interpolation
prob_interp = ODEProblem(HighFidelityEphemerisModel.eom_NbodySH_Interp!, x0, tspan, parameters)
sol_interp = solve(prob_interp, Vern8(), reltol=1e-14, abstol=1e-14)