Empirical SRP: ECOM¶

The cannonball is satkit's default solar-radiation-pressure model. For GNSS-class precise orbits the empirical CODE Orbit Model (ECOM) is available as an addition to it: a small set of constant and harmonic accelerations in a Sun-oriented frame whose coefficients are estimated in orbit determination rather than derived from a surface model.

Experimental

The ECOM interface (EcomParams / ecomparams, satproperties(ecom=...), SatProperties::srp_ecom) is new and marked experimental: it may be reshaped in a minor release (for example into a general empirical-acceleration hook) as it finds users. The physics and conventions below are stable; the API surface is not yet covered by the usual compatibility promise.

Model¶

ECOM is the empirical SRP parameterization used by CODE and most IGS analysis centres for GNSS precise orbit determination (Beutler et al. 1994; Springer et al. 1999; Arnold et al. 2015). It expresses the non-gravitational acceleration of a nominally yaw-steering satellite in a Sun-oriented DYB frame:

\[ \hat{e}_D = \frac{\vec{r}_\odot - \vec{r}}{|\vec{r}_\odot - \vec{r}|}, \qquad \hat{e}_Y = \frac{\hat{e}_D \times \hat{r}}{|\hat{e}_D \times \hat{r}|}, \qquad \hat{e}_B = \hat{e}_D \times \hat{e}_Y \]

(\(\hat{e}_D\) points from the satellite to the Sun; \(\hat{e}_Y\) is the solar-panel rotation axis.) The acceleration is

\[ \vec{a}_\text{ECOM} = \nu \left(\frac{\text{AU}}{d}\right)^2 \left[ D(\varphi)\,\hat{e}_D + Y(\varphi)\,\hat{e}_Y + B(\varphi)\,\hat{e}_B \right] \]

\[ \begin{aligned} D(\varphi) &= D_0 + D_c\cos\varphi + D_s\sin\varphi + D_{2c}\cos 2\varphi + D_{2s}\sin 2\varphi + D_{4c}\cos 4\varphi + D_{4s}\sin 4\varphi \\ Y(\varphi) &= Y_0 + Y_c\cos\varphi + Y_s\sin\varphi \\ B(\varphi) &= B_0 + B_c\cos\varphi + B_s\sin\varphi \end{aligned} \]

Coefficients¶

Coefficient	Axis	Term	In which model	Typical size (GPS, nm/s²)	Field name
\(D_0\)	\(\hat{e}_D\) (toward the Sun)	constant	all	\(-80\) to \(-110\) (≈ \(-P_\odot C_R A/m\); negative because \(\hat{e}_D\) points at the Sun)	`d0`
\(Y_0\)	\(\hat{e}_Y\) (solar-panel axis)	constant	all	\(\sim 1\) (attitude/thermal "Y-bias")	`y0`
\(B_0\)	\(\hat{e}_B\)	constant	all	\(\sim 1\)–\(5\), varies with \(\beta\)	`b0`
\(D_c, D_s\)	\(\hat{e}_D\)	\(\cos\varphi, \sin\varphi\)	ECOM1 only (\(\varphi = u\))	\(\lesssim 1\)	`dc`, `ds`
\(Y_c, Y_s\)	\(\hat{e}_Y\)	\(\cos\varphi, \sin\varphi\)	ECOM1 only	\(\lesssim 1\)	`yc`, `ys`
\(B_c, B_s\)	\(\hat{e}_B\)	\(\cos\varphi, \sin\varphi\)	reduced ECOM and ECOM1 (\(\varphi = u\)); ECOM2 as \(B_{1c}, B_{1s}\) (\(\varphi = \Delta u\))	\(\lesssim 2\)	`bc`, `bs`
\(D_{2c}, D_{2s}\)	\(\hat{e}_D\)	\(\cos 2\Delta u, \sin 2\Delta u\)	ECOM2	few, mostly in eclipse seasons	`d2c`, `d2s`
\(D_{4c}, D_{4s}\)	\(\hat{e}_D\)	\(\cos 4\Delta u, \sin 4\Delta u\)	ECOM2	few, mostly in eclipse seasons	`d4c`, `d4s`

The constructors set the fields for you: reduced(d0, y0, b0, bc, bs), ecom1(d0, y0, b0, dc, ds, yc, ys, bc, bs) and ecom2(d0, y0, b0, b1c, b1s, d2c, d2s, d4c, d4s) (ECOM2's \(B_{1c}, B_{1s}\) are stored in bc, bs; ecom2 sets sun_relative=True). Any coefficient left at zero costs nothing, so a 7-parameter ECOM2 is ecom2(..., d4c=0, d4s=0). All values are accelerations in m/s² at 1 AU (see Sun-distance scaling).

Conventions¶

where \(d\) is the satellite–Sun distance (Sun-distance scaling) and \(\nu\) is the same shadow function as the cannonball term, applied to all three axes — the CODE/Bernese convention ("the acceleration due to the solar radiation pressure is switched off when the satellite is in the Earth's shadow", Bernese GNSS Software v5.2 §2.2.2.3), so coefficients taken from CODE products keep their meaning. The argument \(\varphi\) is selected by sun_relative:

`sun_relative`	\(\varphi\)	Model family
`False`	argument of latitude \(u\) from the ascending node (x-axis projection for an equatorial orbit)	ECOM1 (Beutler 1994), reduced ECOM (Springer 1999)
`True`	\(\Delta u = u - u_\odot\) measured from orbit noon — zero at the point closest to the Sun's projection into the orbit plane, \(\pi\) at midnight; computed node-free and regular at all inclinations	ECOM2 (Arnold 2015)

Because \(\hat{e}_D\) points at the Sun, the physical \(D_0\) is negative — about \(-1\times10^{-7}\) m/s² for a GPS satellite (\(C_R A/m \approx 0.02\) m²/kg), and 10–30 nm/s² when ECOM is applied as a residual on top of an a-priori model. \(Y_0\) and the B terms are typically \(\sim10^{-9}\) m/s². The coefficients are estimated in orbit determination; satkit propagates with the values you supply and adds the ECOM term to the cannonball, so use craoverm=0 for a pure ECOM model:

Sun-distance scaling¶

The whole ECOM acceleration is scaled by \((\text{AU}/d)^2\), with \(d\) the satellite–Sun distance, exactly like the cannonball pressure. The solar flux, and with it every radiation-pressure coefficient, varies by ±3.4 % over the year with Earth's orbital eccentricity. Coefficients referred to 1 AU stay constant across seasons, and \(D_0 = -P_\odot\,C_R A/m\) (-consts.solar_pressure_1au * craoverm) reproduces the cannonball term exactly, at any epoch.

The effect shows up in long predictions. For GPS G20, the reduced ECOM was re-fitted on rolling 3-day windows from January to March 2024. Unscaled, \(D_0\) drifted by +1.5 nm/s² (1.4 %), and three-quarters of that was the \(1/d^2\) change in flux. With the scaling, the error of a 60-day prediction from a 3-day fit drops by 11 % (622 → 551 m), and the best single coefficient set for the whole two months fits twice as well (9.0 → 4.2 m RMS). The remaining drift follows the Sun elevation β over the orbit plane, which constant coefficients cannot capture.

Public descriptions of ECOM (Bernese GNSS Software v5.2 §2.2.2.3; Arnold et al. 2015) write the model without a distance factor, and implementations differ: Orekit's ECOM2 applies none, while Ginan scales the D/Y/B terms by \((\text{AU}/d)^2\) as satkit does. Over the 1–3-day arcs on which ECOM is normally estimated the factor changes by ≤ 0.1 %, so the choice only matters when coefficients move between programs or are used over long spans. To use coefficients estimated with an unscaled implementation, multiply them by \((d/\text{AU})^2\) at the arc epoch. satkit 0.23.1 and earlier applied ECOM coefficients unscaled.

Usage¶

import satkit as sk

# Reduced ECOM (D0, Y0, B0, Bc, Bs) in m/s^2 — e.g. from your own fit or a CODE product
ecom = sk.ecomparams.reduced(-1.06e-7, 1.0e-9, -3.2e-9, 1.2e-9, 0.3e-9)
props = sk.satproperties(craoverm=0.0, ecom=ecom)
res = sk.propagate(state, t0, t1, propsettings=settings, satproperties=props)

# ECOM2: D0, Y0, B0, B1c, B1s, D2c, D2s, D4c, D4s (argument Δu from orbit noon)
ecom2 = sk.ecomparams.ecom2(-1.06e-7, 1e-9, -3e-9, 0, 0, -2e-9, -3.6e-9, 1.9e-9, -0.7e-9)

Rust users can also implement SatProperties::srp_ecom(&self, tm, state) -> Option<EcomParams> to supply coefficients that change over the propagation (per-arc CODE values, an attitude-mode switch). Like the cannonball term, ECOM contributes no partials to the state transition matrix.

What to expect¶

Fitting an initial state plus the reduced 5-parameter ECOM to 3 days of IGS final GPS orbits (python/examples/ecom_gps_validation.py, 12×12 gravity, no a-priori box-wing) reproduces the orbit at the accuracy of the IGS product itself — 5 cm 3D RMS versus 3.8 m for the cannonball — and predicts the next 24 h to a median 6–7 cm 3D across the constellation (2-day fits, ten satellites), against ~5 cm for the IGS ultra-rapid predicted product and 8–10 cm reported by Duan & Hugentobler (2021) from 3-day arcs with a full analysis-centre force model. Beyond a few days the error grows along-track roughly as \(t^2\): ~0.9 m at 7 days, 10 m after ~12 days and ~95 m at 30 days from a 3-day fit (~15 days and ~60 m from a 7-day fit), because the true coefficients drift with the Sun elevation angle β over weeks — which is exactly why analysis centres re-estimate them every day. Constant ECOM coefficients are a short-arc (days) model, not a month-long one. The worked example notebook walks through the fit, the prediction and the constellation benchmark with plots.

Practical notes¶

SP3 epochs are in GPS time (scale=timescale.GPS), not UTC — reading them as UTC rotates the truth by 18 s of Earth rotation relative to the Sun/Moon geometry and quadruples the fit residual. And for fits whose arcs cross Earth's shadow, use integrator = gauss_jackson8. The adaptive Runge–Kutta steppers integrate through eclipses without trouble, but where their steps fall relative to the shadow edges shifts with every small change of the initial state, so the propagated orbit is not a smooth function of it: over 3 days of G08 (7 % in shadow) at 1e-11 tolerance the end position wanders by ~4 cm, against 0.003 mm on the fixed Gauss–Jackson grid or without eclipses. That noise corrupts finite-difference Jacobians (an rkv98 fit of G08 stalls at 5.6 m RMS), while Gauss–Jackson 8 fits the eclipsing satellite as well as any other (4.7 cm fit, 5.6 cm at 24 h). Plain propagation through eclipses is fine with either integrator.