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[ 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².

Conventions¶

where \(\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:

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.8 m at 7 days, 10 m after ~12 days and ~100 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 ECOM Solar Radiation Pressure tutorial 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 arcs that cross Earth's shadow, use integrator = gauss_jackson8: the adaptive Runge–Kutta steppers can abort at a shadow boundary with too many consecutive step rejections, whereas the fixed-step multistep integrator is immune and fits an eclipsing satellite just as well (G08, 8% umbra: 4.7 cm fit, 5.6 cm at 24 h).