Solid tides

[ChristopherRabotin]

High level description

Solid tides modify the gravity field to account for the crust deformation due to the Moon and the Sun (in the case of Earth). This should be supported in Nyx, likely by building an AccelModel.

Requirements

• Support SolidTides as eom and dual_eom. Note that the solid tides only depend on the central body position, and not the spacecraft, so its partial is constant.

Test plans

Unsure.

Design

Let me read the uploaded files first before giving any advice.Good, I have a thorough picture of the code. Here's a detailed breakdown of what needs to be done and how.

---

Architecture of the change

The cleanest approach is to implement solid tides as a separate AccelModel — a SolidTides struct — rather than modifying Harmonics directly. This matches Nyx's existing compositional force model pattern, keeps Harmonics clean, and lets users opt in explicitly. The tidal corrections are computed as $\Delta\bar{C}{nm}$ / $\Delta\bar{S}{nm}$, then fed through the same spherical harmonic acceleration machinery already in sph_harmonics.rs.

---

What the new struct needs

pub struct SolidTides {
    compute_frame: Frame,    // ITRF / Earth body-fixed frame
    eq_radius_km: f64,       // R_E, from the gravity model
    mu_km3_s2: f64,          // GM_Earth
    k2: f64,                 // Elastic Love number, ~0.30
    k3: f64,                 // Degree-3 Love number, ~0.093
    // Frames for querying Moon and Sun positions from the almanac
    moon_frame: Frame,
    sun_frame: Frame,
}

The Love numbers are constants you can hard-code from the IERS 2010 Conventions (Table 6.3): $k_2 = 0.30190$, $k_3 = 0.093$.

---

The core computation — eom

The IERS Conventions §6.2 procedure, translated to what this code needs:

Step 1 — Get Moon and Sun positions in ECEF

let moon_state = almanac.transform_to(moon_orbit, self.compute_frame, None)?;
let sun_state  = almanac.transform_to(sun_orbit,  self.compute_frame, None)?;

You need their ECEF position vectors $\mathbf{r}_{\text{body}}$, geocentric latitude $\phi$, and longitude $\lambda$.

Step 2 — Compute $\Delta\bar{C}{nm}$, $\Delta\bar{S}{nm}$

For each tide-generating body and for degree $n \in {2, 3}$:

$$\Delta\bar{C}_{nm} - i,\Delta\bar{S}_{nm} = \frac{k_n}{2n+1} \cdot \frac{GM_{\text{body}}}{GM_\oplus} \left(\frac{R_E}{r_{\text{body}}}\right)^{n+1} \bar{P}_{nm}(\sin\phi_{\text{body}}) \cdot e^{-im\lambda_{\text{body}}}$$

Expanding $e^{-im\lambda} = \cos(m\lambda) - i\sin(m\lambda)$:

let delta_c_nm =  k_n / (2*n+1) * gm_ratio * ratio_pow * p_nm * cos_ml;
let delta_s_nm = -k_n / (2*n+1) * gm_ratio * ratio_pow * p_nm * sin_ml;

where gm_ratio = GM_body / GM_earth, ratio_pow = (R_E / r_body)^(n+1), and p_nm is the normalized associated Legendre polynomial evaluated at sin(lat_body).

You sum contributions from Moon and Sun into a small array of delta coefficients — only degrees 2 and 3, so at most ~16 (n,m) pairs.

Step 3 — Feed into the harmonic acceleration kernel

Rather than re-implementing the full Legendre recursion, build a tiny HarmonicsMem-like structure populated only with the $\Delta\bar{C}{nm}$, $\Delta\bar{S}{nm}$ values, then call the existing eom inner loop logic. Or more pragmatically for a first pass: compute the tidal potential gradient directly using the same rho_np1, a_nm, r_m/i_m machinery from sph_harmonics.rs, but substituting the delta coefficients in the inner loop where cs_nm(n, m) is called.

The cleanest path given Nyx's structure is to make HarmonicsMem mutable or provide an additive override path — or simply have SolidTides::eom duplicate the small inner kernel (it only loops over n=2,3 so it's trivial in size).

---

The Legendre polynomial you need

The normalized associated Legendre polynomials $\bar{P}_{nm}(\sin\phi)$ for just n=2,3 can be computed analytically or via a short recursion. For n=2:

P_20(x) = sqrt(5)/2   * (3x² - 1)
P_21(x) = sqrt(5/3)   * 3x*sqrt(1-x²)       [× 1/sqrt(2) for m=0 normalisation handled separately]
P_22(x) = sqrt(5/12)  * 3(1-x²)

where $x = \sin\phi_{\text{body}}$. These are the fully normalized forms consistent with the rest of the code's convention.

---

Wiring into AccelModel

impl AccelModel for SolidTides {
    fn eom(&self, osc: &Orbit, almanac: Arc<Almanac>) 
        -> Result<Vector3<f64>, DynamicsError> 
    {
        // 1. Query Moon and Sun positions in compute_frame (ECEF)
        // 2. Compute delta_c_nm, delta_s_nm for n=2,3, bodies=Moon+Sun
        // 3. Run the harmonic acceleration inner loop with only those corrections
        // 4. Rotate result back to integration frame via almanac.rotate(...)
    }
    
    fn dual_eom(&self, osc: &Orbit, almanac: Arc<Almanac>) 
        -> Result<(Vector3<f64>, Matrix3<f64>), DynamicsError> 
    {
        // Same but with OHyperdual — the delta coefficients are real-valued constants
        // at a given epoch, so they promote trivially: OHyperdual::from(delta_c_nm)
    }
}

---

Things to be careful about

Tide system consistency. The EGM2008 file loaded in gravity.rs is EGM2008_to2190_TideFree — the tide-free convention. The IERS solid tide correction formula assumes a zero-tide or tide-free reference depending on which permanent tide term you include. For tide-free EGM2008, you should use the tide-free Love number formulation and not add back the permanent tide. The IERS Conventions §6.2.1 cover this distinction explicitly — the simplest safe choice is to use the step-1/step-2 procedure from Table 6.3 directly, which is designed for use with tide-free models.

GM values for Moon and Sun. These need to come from your ANISE almanac / ephemeris constants, not be hard-coded. ANISE already carries DE440 constants so almanac.frame_from_uid(MOON_J2000)?.mu_km3_s2() should give you what you need.

The dual_eom STM. The tidal delta coefficients are computed from the body positions, not from the satellite position. So for the STM calculation they are constants at a given epoch — you don't need to hyperdual-differentiate through the Legendre evaluation of the Moon/Sun positions. Only the satellite position derivatives (through r_, s_, t_, u_, a_nm) need the hyperdual treatment, identical to what's already in sph_harmonics.rs.

Frequency-dependent corrections to $\bar{C}_{20}$. The IERS table provides ~12 additive corrections driven by Doodson/Delaunay tidal arguments. These require computing the fundamental arguments of nutation (linear functions of TDB time). This is a second-pass refinement — get the basic $k_2$/$k_3$ step working first, then layer these on.

---

Summary of files to touch

• New file solid_tides.rs — the SolidTides struct and its AccelModel impl. Roughly 150–200 lines.
• gravity.rs — no changes needed; HarmonicsMem is fine as-is.
• sph_harmonics.rs — no changes needed if you duplicate the small inner loop; minor refactor if you want to share it.
• mod.rs / dynamics exports — expose SolidTides alongside Harmonics.