Deriving the Gauss Variational Equations

Author

Devendra Ghate

Published

October 12, 2026

1 Overview

The Gauss variational equations give the rates of change of the six osculating elements in terms of the perturbing acceleration’s radial, along-track and cross-track components. All six follow from one fact: a perturbation changes the velocity but never the position, so it acts on the elements like a continuous stream of infinitesimal impulses.

The problem is the perturbed two-body equation, with the perturbing specific force resolved in the RSW frame:

\[ \ddot{\mathbf r} = -\frac{\mu}{r^{3}}\,\mathbf r + \mathbf f, \qquad \mathbf f = f_R\,\hat{\mathbf R} + f_S\,\hat{\mathbf S} + f_W\,\hat{\mathbf W} \]

The result has the structure below. The 6 × 3 matrix B depends only on the current elements, and only the mean anomaly carries a Keplerian drift. The Summary section lists B in full.

\[ \frac{d}{dt} \begin{pmatrix} a \\ e \\ i \\ \Omega \\ \omega \\ M \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ n \end{pmatrix} + \mathsf{B}(a, e, i, \omega, M) \begin{pmatrix} f_R \\ f_S \\ f_W \end{pmatrix} \]

The Gauss form needs only the force, so it covers drag, radiation pressure and thrust as readily as gravity. Lagrange’s planetary equations need a disturbing potential. When one exists, writing the force as its gradient turns the Gauss form into Lagrange’s.

The derivation has three steps.

  1. Show that a perturbation changes any Keplerian integral only through that integral’s velocity gradient.
  2. For each element, differentiate the state function whose velocity gradient is simplest: angular momentum for h, p, i and Ω; energy for a; the pair e cos ν and e sin ν for e and ν; the difference u − ν for ω; Kepler’s equation for M.
  3. Cross-check the results against the eccentricity vector, a J2 orbit average and a finite-difference test.
flowchart TD
  R["<b>Fundamental rule</b><br/>rate of g = Kepler part + (∂g/∂v)·f"]
  EN["<b>Energy</b><br/>v²/2 − μ/r"]
  H["<b>Angular momentum</b><br/>h = r × v, torque r × f"]
  P["<b>Anomaly pair</b><br/>(e cos ν, e sin ν)"]
  A["<b>Semi-major axis a</b><br/>only the power v·f"]
  HP["<b>h and p</b><br/>only along-track"]
  IO["<b>i and Ω</b><br/>only cross-track"]
  EV["<b>e and ν</b><br/>only in-plane forces"]
  U["<b>u = ω + ν</b><br/>argument of latitude"]
  M["<b>Mean anomaly M</b><br/>chain rule on M(ν, e)"]
  W["<b>Argument of periapsis ω</b><br/>u − ν: Kepler parts cancel"]
  X["<b>Eccentricity vector (cross-check)</b><br/>confirms the e, ω, i and Ω rates"]
  R --> EN & H & P
  EN --> A
  H --> HP & IO
  P --> EV
  IO -- node shift --> U
  EV --> M
  U --> W
  EV --> W
  X -.-> W
Figure 1: Derivation map: from state functions to elements

Each arrow runs from a state function to the rates it yields. The dashed link marks the eccentricity-vector cross-check of the e, ω, i and Ω rates.

2 Osculating elements and the fundamental rule

A perturbing force changes a Keplerian integral only through that integral’s gradient with respect to velocity. This one rule generates every Gauss equation.

2.1 Variation of parameters

Without the perturbation, the motion is a conic fixed by six constants, the element vector α:

\[ \mathbf r = \boldsymbol\phi(t;\boldsymbol\alpha), \qquad \mathbf v = \boldsymbol\psi(t;\boldsymbol\alpha) \equiv \frac{\partial\boldsymbol\phi}{\partial t}, \qquad \frac{\partial\boldsymbol\psi}{\partial t} = -\frac{\mu}{\phi^{3}}\,\boldsymbol\phi \]

With the perturbation on, keep the same functional form but let the constants vary:

\[ \mathbf r(t) = \boldsymbol\phi\big(t;\boldsymbol\alpha(t)\big) \quad\Longrightarrow\quad \dot{\mathbf r} = \boldsymbol\psi + \frac{\partial\boldsymbol\phi}{\partial\boldsymbol\alpha}\,\dot{\boldsymbol\alpha} \]

Six unknown functions against three equations of motion leave three conditions free. The osculating choice demands that the velocity also keep its Keplerian form, which is Lagrange’s constraint:

\[ \frac{\partial\boldsymbol\phi}{\partial\boldsymbol\alpha}\,\dot{\boldsymbol\alpha} = \mathbf 0 \quad\Longleftrightarrow\quad \dot{\mathbf r} = \boldsymbol\psi\big(t;\boldsymbol\alpha(t)\big) \]

Differentiate the velocity once more and substitute the equation of motion. The Keplerian acceleration cancels on both sides:

\[ \ddot{\mathbf r} = \underbrace{\frac{\partial\boldsymbol\psi}{\partial t}}_{-\mu\mathbf r/r^{3}} + \frac{\partial\boldsymbol\psi}{\partial\boldsymbol\alpha}\,\dot{\boldsymbol\alpha} \quad\Longrightarrow\quad \frac{\partial\boldsymbol\psi}{\partial\boldsymbol\alpha}\,\dot{\boldsymbol\alpha} = \mathbf f \]

Stacking the two conditions gives a 6 × 6 linear system whose matrix is the Jacobian of the state with respect to the elements. Its inverse is the Jacobian of the elements with respect to the state, and the zero block kills the position gradient:

\[ \begin{bmatrix} \partial\boldsymbol\phi/\partial\boldsymbol\alpha \\ \partial\boldsymbol\psi/\partial\boldsymbol\alpha \end{bmatrix} \dot{\boldsymbol\alpha} = \begin{bmatrix} \mathbf 0 \\ \mathbf f \end{bmatrix} \quad\Longrightarrow\quad \boxed{\;\dot{\boldsymbol\alpha} = \frac{\partial\boldsymbol\alpha}{\partial\mathbf r}\,\mathbf 0 + \frac{\partial\boldsymbol\alpha}{\partial\mathbf v}\,\mathbf f = \frac{\partial\boldsymbol\alpha}{\partial\mathbf v}\,\mathbf f\;} \]

The three osculating conditions are a gauge choice. They make the elements describe the conic tangent to the true path: the orbit the body would fly if the perturbation vanished at that instant. Other gauges are legitimate but give non-osculating elements.

2.2 The fundamental rule

The same result holds for any smooth state function g of position and velocity. The chain rule along the true path splits its rate into a Keplerian part and a perturbation part:

\[ \dot g = \underbrace{\frac{\partial g}{\partial \mathbf r}\cdot\mathbf v - \frac{\mu}{r^{3}}\,\frac{\partial g}{\partial \mathbf v}\cdot\mathbf r}_{\textstyle \mathcal K g} \;+\; \underbrace{\frac{\partial g}{\partial \mathbf v}\cdot\mathbf f}_{\textstyle \mathcal D g} \]

For a Keplerian integral (a, e, i, Ω, ω, the angular momentum, the energy) the Keplerian part vanishes, so the rate is just 𝒟g. The fast angles keep a known Keplerian part:

\[ \mathcal K\nu = \mathcal K u = \frac{h}{r^{2}}, \qquad \mathcal K M = n \]

The operator 𝒟 is linear and obeys the product and chain rules. Its action on the basic quantities is

\[ \mathcal D\mathbf r = \mathbf 0, \qquad \mathcal D\mathbf v = \mathbf f, \qquad \mathcal D\,G(\mathbf r) = 0 \;\;\text{for any function of position alone.} \]

Physically, 𝒟g dt is the change in g caused by the velocity increment f dt applied at a fixed position. The Gauss equations are therefore the first-order impulsive-maneuver equations divided by dt, and a small impulse changes any element by

\[ \Delta\alpha \simeq \frac{\partial\alpha}{\partial\mathbf v}\cdot\Delta\mathbf v \]

Nothing here assumes the force is conservative or independent of velocity, so drag and thrust are covered. The rest of the derivation evaluates 𝒟 on well-chosen state functions.

3 Frames, notation and identities

Forces are resolved in the RSW frame, which rides with the body: R along the radius, W along the angular momentum, and S completing the triad in the direction of motion. The same frame is often called RTN or the Hill frame.

\[ \hat{\mathbf R} = \frac{\mathbf r}{r}, \qquad \hat{\mathbf W} = \frac{\mathbf h}{h} = \frac{\mathbf r\times\mathbf v}{\lVert \mathbf r\times\mathbf v\rVert}, \qquad \hat{\mathbf S} = \hat{\mathbf W}\times\hat{\mathbf R} \]

The triad is right-handed, so cross products cycle R → S → W:

\[ \hat{\mathbf R}\times\hat{\mathbf S} = \hat{\mathbf W}, \qquad \hat{\mathbf S}\times\hat{\mathbf W} = \hat{\mathbf R}, \qquad \hat{\mathbf W}\times\hat{\mathbf R} = \hat{\mathbf S} \]

3.1 Notation

The elements are a, e, i, Ω, ω and M, with true anomaly ν and eccentric anomaly E. The derived quantities are

\[ \begin{aligned} u &= \omega + \nu && \text{argument of latitude} \\ p &= a\,(1-e^{2}) = h^{2}/\mu && \text{semi-latus rectum} \\ n &= \sqrt{\mu/a^{3}} && \text{mean motion} \\ b &= a\sqrt{1-e^{2}}, \qquad h = n\,a\,b && \text{semi-minor axis} \end{aligned} \]

3.2 Unit vectors

With K the inertial polar axis, n̂ the ascending-node direction and m̂ its in-plane companion 90° ahead, the inertial components are

\[ \hat{\mathbf n} = (\cos\Omega,\ \sin\Omega,\ 0), \qquad \hat{\mathbf m} = \hat{\mathbf W}\times\hat{\mathbf n} = (-\cos i\sin\Omega,\ \cos i\cos\Omega,\ \sin i), \qquad \hat{\mathbf W} = (\sin i\sin\Omega,\ -\sin i\cos\Omega,\ \cos i) \]

The radial and along-track directions sit at angle u from the node, and the perifocal directions P (to periapsis) and Q sit at angle ω:

\[ \begin{aligned} \hat{\mathbf R} &= \cos u\,\hat{\mathbf n} + \sin u\,\hat{\mathbf m}, & \hat{\mathbf S} &= -\sin u\,\hat{\mathbf n} + \cos u\,\hat{\mathbf m}, \\ \hat{\mathbf P} &= \cos\omega\,\hat{\mathbf n} + \sin\omega\,\hat{\mathbf m} = \cos\nu\,\hat{\mathbf R} - \sin\nu\,\hat{\mathbf S}, & \hat{\mathbf Q} &= -\sin\omega\,\hat{\mathbf n} + \cos\omega\,\hat{\mathbf m} = \sin\nu\,\hat{\mathbf R} + \cos\nu\,\hat{\mathbf S} \end{aligned} \]

Projecting onto the polar axis gives the three components used for i, Ω and u:

\[ \hat{\mathbf K}\cdot\hat{\mathbf R} = \sin i\sin u, \qquad \hat{\mathbf K}\cdot\hat{\mathbf S} = \sin i\cos u, \qquad \hat{\mathbf K}\cdot\hat{\mathbf W} = \cos i, \qquad \hat{\mathbf K} = \cos i\,\hat{\mathbf W} + \sin i\,\hat{\mathbf m} \]

3.3 Orbit geometry and velocity

On the osculating conic, the radius and the velocity components are

\[ r = \frac{p}{1+e\cos\nu}, \qquad \mathbf v = v_R\,\hat{\mathbf R} + v_S\,\hat{\mathbf S}, \qquad v_R = \frac{\mu}{h}\,e\sin\nu, \qquad v_S = \frac{h}{r} = \frac{\mu}{h}\,(1+e\cos\nu) \]

The energy integral and the eccentric-anomaly relations complete the toolkit:

\[ v^{2} = \mu\left(\frac{2}{r} - \frac{1}{a}\right), \qquad r = a\,(1-e\cos E), \qquad \cos E = \frac{e+\cos\nu}{1+e\cos\nu}, \qquad \sin E = \frac{\sqrt{1-e^{2}}\,\sin\nu}{1+e\cos\nu} \]

3.4 Three identities

Three identities do most of the algebra. (I1) is the orbit equation rearranged, (I2) restates the definition of p, and (I3) follows from the expression for cos E.

\[ \begin{aligned} &p - r = r\,e\cos\nu && \text{(I1)} \\[4pt] &\frac{h}{\mu} = \frac{p}{h} && \text{(I2)} \\[4pt] &\frac{r}{p}\,(e+\cos\nu) = \cos E && \text{(I3)} \end{aligned} \]

4 Angular momentum: h, p and the orbit normal

Only the along-track force changes the size of the angular momentum, and only the cross-track force tilts it.

Angular momentum is a Keplerian integral, so its rate is 𝒟 of it. Position is untouched by 𝒟, which leaves the torque of the perturbing force:

\[ \dot{\mathbf h} = \mathcal D(\mathbf r\times\mathbf v) = \mathbf r\times\mathcal D\mathbf v = \mathbf r\times\mathbf f = r\,\hat{\mathbf R}\times\big(f_R\,\hat{\mathbf R} + f_S\,\hat{\mathbf S} + f_W\,\hat{\mathbf W}\big) = r f_S\,\hat{\mathbf W} - r f_W\,\hat{\mathbf S} \]

Projecting on the normal gives the magnitude; the remainder gives the rotation of the orbit normal:

\[ \boxed{\;\dot h = \hat{\mathbf W}\cdot\dot{\mathbf h} = r f_S\;} \qquad \dot{\hat{\mathbf W}} = \frac{\dot{\mathbf h} - \dot h\,\hat{\mathbf W}}{h} = -\frac{r f_W}{h}\,\hat{\mathbf S} \]

The semi-latus rectum follows from its definition:

\[ p = \frac{h^{2}}{\mu} \quad\Longrightarrow\quad \boxed{\;\dot p = \frac{2h}{\mu}\,\dot h = 2r\sqrt{\frac{p}{\mu}}\;f_S\;} \]

A radial force exerts no torque, so it cannot change h or p. It reshapes the orbit only by trading a against e at fixed p. The normal’s rotation feeds the inclination and node equations later.

5 Energy and the semi-major axis

The semi-major axis changes only through the power the perturbation delivers, so only the velocity-aligned part of the force matters.

The specific energy is a Keplerian integral whose potential term depends on position alone. Applying 𝒟 gives the work–energy theorem:

\[ \mathcal E = \tfrac12\,\mathbf v\cdot\mathbf v - \frac{\mu}{r} \quad\Longrightarrow\quad \dot{\mathcal E} = \mathcal D\mathcal E = \mathbf v\cdot\mathcal D\mathbf v = \mathbf v\cdot\mathbf f \]

In RSW components, using the velocity components of the osculating conic,

\[ \dot{\mathcal E} = v_R f_R + v_S f_S = \frac{\mu}{h}\left[e\sin\nu\,f_R + (1+e\cos\nu)\,f_S\right] = \frac{\mu}{h}\left[e\sin\nu\,f_R + \frac{p}{r}\,f_S\right] \]

The energy fixes the semi-major axis, so

\[ a = -\frac{\mu}{2\mathcal E} \quad\Longrightarrow\quad \dot a = \frac{\mu}{2\mathcal E^{2}}\,\dot{\mathcal E} = \frac{2a^{2}}{\mu}\,\mathbf v\cdot\mathbf f \]

\[ \boxed{\;\dot a = \frac{2a^{2}}{h}\left(e\sin\nu\,f_R + \frac{p}{r}\,f_S\right) = \frac{2}{n\sqrt{1-e^{2}}}\left(e\sin\nu\,f_R + \frac{p}{r}\,f_S\right)\;} \]

The second form uses h = nab. On a circular orbit the radial term vanishes and only the along-track force changes a:

\[ e = 0:\qquad \dot a = \frac{2}{n}\,f_S \]

The mean motion follows from a, and its change is what makes the mean anomaly drift later:

\[ n = \sqrt{\mu/a^{3}} \quad\Longrightarrow\quad \dot n = -\frac{3n}{2a}\,\dot a \]

6 Eccentricity and true anomaly

Treating e cos ν and e sin ν as the Cartesian coordinates of one point makes the eccentricity and true-anomaly equations fall out together.

6.1 The two state functions

Both coordinates are functions of the state. The first is the orbit equation solved for e cos ν; the second comes from the radial velocity:

\[ e\cos\nu = \frac{h^{2}}{\mu r} - 1, \qquad e\sin\nu = \frac{h\,v_R}{\mu}, \qquad v_R = \hat{\mathbf R}\cdot\mathbf v \]

The radius and the radial unit vector depend on position alone, and h was handled in the angular-momentum section:

\[ \mathcal D r = 0, \qquad \mathcal D v_R = \hat{\mathbf R}\cdot\mathbf f = f_R, \qquad \mathcal D h = r f_S \]

\[ \mathcal D(e\cos\nu) = \frac{2h}{\mu r}\,\mathcal D h = \frac{2h}{\mu}\,f_S, \qquad \mathcal D(e\sin\nu) = \frac{v_R\,\mathcal D h + h\,\mathcal D v_R}{\mu} = \frac{h f_R + r v_R f_S}{\mu} \]

6.2 Polar decomposition

Rotating the variation of the pair by −ν separates its radial part (the change in e) from its angular part (the change in ν):

\[ \mathcal D e = \cos\nu\,\mathcal D(e\cos\nu) + \sin\nu\,\mathcal D(e\sin\nu), \qquad e\,\mathcal D\nu = \cos\nu\,\mathcal D(e\sin\nu) - \sin\nu\,\mathcal D(e\cos\nu) \]

6.3 Eccentricity

Substitute the two variations, the radial velocity from the osculating conic, and (I2):

\[ \dot e = \frac{2h\cos\nu}{\mu}\,f_S + \frac{\sin\nu\,(h f_R + r v_R f_S)}{\mu} = \frac{1}{h}\Big[p\sin\nu\,f_R + \big(2p\cos\nu + re\sin^{2}\nu\big)f_S\Big] \]

Identity (I1) puts the along-track coefficient into its standard form:

\[ 2p\cos\nu + re\sin^{2}\nu = 2p\cos\nu + re - \underbrace{re\cos\nu}_{p-r}\,\cos\nu = (p+r)\cos\nu + re \]

\[ \boxed{\;\dot e = \frac{1}{h}\Big\{p\sin\nu\,f_R + \big[(p+r)\cos\nu + re\big]f_S\Big\} = \sqrt{\frac{p}{\mu}}\,\Big[\sin\nu\,f_R + (\cos\nu + \cos E)\,f_S\Big]\;} \]

The second form uses (I3). Unlike the angle rates below, the eccentricity rate stays finite as e → 0.

6.4 True anomaly

The angular part uses the same substitutions, and (I1) collapses the along-track coefficient:

\[ e\,\mathcal D\nu = \frac{\cos\nu\,(h f_R + r v_R f_S)}{\mu} - \frac{2h\sin\nu}{\mu}\,f_S = \frac{1}{h}\Big[p\cos\nu\,f_R - \big(\underbrace{2p - re\cos\nu}_{p+r}\big)\sin\nu\,f_S\Big] \]

Adding the Keplerian part gives the full rate:

\[ \boxed{\;\dot\nu = \frac{h}{r^{2}} + \frac{1}{he}\Big[p\cos\nu\,f_R - (p+r)\sin\nu\,f_S\Big]\;} \]

The perturbative part is singular at e = 0, where periapsis, and with it ν, is undefined.

7 Orientation of the orbital plane: i and Ω

Only the cross-track force moves the orbital plane. It does so by rotating the whole conic about the current radius vector.

7.1 Inclination

The inclination is the angle between the orbit normal and the polar axis. Differentiate its cosine and use the normal’s rotation from the angular-momentum section:

\[ \cos i = \hat{\mathbf K}\cdot\hat{\mathbf W} \quad\Longrightarrow\quad -\sin i\;\dot i = \hat{\mathbf K}\cdot\dot{\hat{\mathbf W}} = -\frac{r f_W}{h}\,\hat{\mathbf K}\cdot\hat{\mathbf S} = -\frac{r f_W}{h}\,\sin i\cos u \]

\[ \boxed{\;\dot i = \frac{r\cos u}{h}\,f_W\;} \]

7.2 Node

The node is the azimuth of the normal’s equatorial projection. With the components of the normal from the notation section,

\[ h_x = h\sin i\sin\Omega, \quad h_y = -h\sin i\cos\Omega \quad\Longrightarrow\quad \dot\Omega = \frac{h_x\dot h_y - h_y\dot h_x}{h_x^{2} + h_y^{2}} = \frac{\hat{\mathbf K}\cdot(\mathbf h\times\dot{\mathbf h})}{h^{2}\sin^{2} i} \]

The cross product needs only the torque found earlier:

\[ \mathbf h\times\dot{\mathbf h} = h\,\hat{\mathbf W}\times\big(r f_S\,\hat{\mathbf W} - r f_W\,\hat{\mathbf S}\big) = h\,r f_W\,\hat{\mathbf R}, \qquad \hat{\mathbf K}\cdot\hat{\mathbf R} = \sin i\sin u \]

\[ \boxed{\;\dot\Omega = \frac{r\sin u}{h\sin i}\,f_W\;} \]

7.3 Geometric view: rigid rotation about the radius vector

A cross-track impulse leaves the position fixed and tilts the velocity. A small rotation of the entire state about the radial direction does exactly the same thing:

\[ \begin{aligned} \text{impulse:}\quad &(\mathbf r,\ \mathbf v) \;\mapsto\; \big(\mathbf r,\ \mathbf v + f_W\,dt\,\hat{\mathbf W}\big) \\ \text{rotation by } \delta\phi \text{ about } \hat{\mathbf R}:\quad &(\mathbf r,\ \mathbf v) \;\mapsto\; \big(\mathbf r,\ \mathbf v + \delta\phi\,\hat{\mathbf R}\times\mathbf v\big) = \big(\mathbf r,\ \mathbf v + \delta\phi\,v_S\,\hat{\mathbf W}\big) \end{aligned} \]

The two agree to first order for the angle below. The Kepler problem is rotation-invariant, so the new conic is the old one turned rigidly: a, e and ν are untouched, and the orientation spins with angular velocity

\[ \delta\phi = \frac{f_W\,dt}{v_S} = \frac{r f_W}{h}\,dt \qquad\Longrightarrow\qquad \boldsymbol\omega_{\perp} = \frac{r f_W}{h}\,\hat{\mathbf R} \]

Resolve this along the axes of the 3-1-3 Euler sequence: Ω about K, i about the node line, ω about the normal. Substituting the unit-vector relations from the notation section,

\[ \boldsymbol\omega_{\perp} = \frac{r f_W}{h}\left[\cos u\,\hat{\mathbf n} + \sin u\,\frac{\hat{\mathbf K} - \cos i\,\hat{\mathbf W}}{\sin i}\right] \equiv \dot\Omega\,\hat{\mathbf K} + \dot i\,\hat{\mathbf n} + \dot\omega_{\perp}\,\hat{\mathbf W} \]

Matching coefficients reproduces both boxed results and adds a third. The cross-track force also moves ω, because the node from which ω is measured slides along the orbit:

\[ \dot\Omega = \frac{r\sin u}{h\sin i}\,f_W, \qquad \dot i = \frac{r\cos u}{h}\,f_W, \qquad \dot\omega_{\perp} = -\frac{r\sin u\cos i}{h\sin i}\,f_W = -\cos i\;\dot\Omega \]

Per unit cross-track impulse, inclination responds most at the nodes and the node responds most 90° from them. Both effects grow with r, which is why plane changes are done far out when the geometry allows.

8 Argument of periapsis: ω

The argument of periapsis moves for two reasons. In-plane forces turn the apse line within the orbital plane, and the cross-track force shifts the node from which ω is measured.

8.1 From the difference u − ν

The argument of periapsis is the argument of latitude minus the true anomaly. Both angles carry the same Keplerian part h/r², so it cancels and only velocity gradients remain:

\[ \omega = u - \nu \quad\Longrightarrow\quad \dot\omega = (\mathcal K u + \mathcal D u) - (\mathcal K\nu + \mathcal D\nu) = \mathcal D u - \mathcal D\nu \]

The true-anomaly gradient comes from the previous section. For u, the polar-axis components of Ŝ and R̂ form a point at polar angle u, the device already used for e and ν:

\[ \big(\hat{\mathbf K}\cdot\hat{\mathbf S},\ \hat{\mathbf K}\cdot\hat{\mathbf R}\big) = \sin i\,(\cos u,\ \sin u) \quad\Longrightarrow\quad \sin^{2} i\;\mathcal D u = (\hat{\mathbf K}\cdot\hat{\mathbf S})\,\mathcal D(\hat{\mathbf K}\cdot\hat{\mathbf R}) - (\hat{\mathbf K}\cdot\hat{\mathbf R})\,\mathcal D(\hat{\mathbf K}\cdot\hat{\mathbf S}) \]

The radial direction depends on position alone, so its term vanishes. The along-track direction inherits the rotation of the normal:

\[ \mathcal D\hat{\mathbf R} = \mathbf 0, \qquad \mathcal D\hat{\mathbf S} = \mathcal D\hat{\mathbf W}\times\hat{\mathbf R} = -\frac{r f_W}{h}\,\hat{\mathbf S}\times\hat{\mathbf R} = \frac{r f_W}{h}\,\hat{\mathbf W} \quad\Longrightarrow\quad \mathcal D(\hat{\mathbf K}\cdot\hat{\mathbf S}) = \frac{r f_W}{h}\,\cos i \]

\[ \mathcal D u = -\frac{r\sin u\cos i}{h\sin i}\,f_W = -\cos i\;\dot\Omega \]

This is the ω̇⊥ of the rigid-rotation picture, as it must be. A rotation about R̂ leaves ν alone, so all of its change in u is a change in ω. Subtracting 𝒟ν gives the full rate:

\[ \boxed{\;\dot\omega = \frac{1}{he}\Big[-p\cos\nu\,f_R + (p+r)\sin\nu\,f_S\Big] - \frac{r\sin u\cos i}{h\sin i}\,f_W\;} \]

With p/h = √(p/μ) from (I2), this becomes the form quoted in most textbooks:

\[ \dot\omega = \frac{1}{e}\sqrt{\frac{p}{\mu}}\,\Big[-\cos\nu\,f_R + \Big(1+\frac{r}{p}\Big)\sin\nu\,f_S\Big] - \frac{r\sin u\cos i}{h\sin i}\,f_W \]

The in-plane terms alone give the rotation rate of the apse line within the orbital plane. This combination returns in the mean-anomaly equation:

\[ \dot\omega + \cos i\;\dot\Omega = \frac{1}{he}\Big[-p\cos\nu\,f_R + (p+r)\sin\nu\,f_S\Big] \]

8.2 Cross-check with the eccentricity vector

The eccentricity vector 𝐞 points to periapsis with length e, so its rate holds ė and the apse-line rotation together. It is a Keplerian integral, and 𝒟 acts only on its velocities:

\[ \mu\,\mathbf e = \mathbf v\times\mathbf h - \mu\,\hat{\mathbf R} = (\mathbf v\cdot\mathbf v)\,\mathbf r - (\mathbf r\cdot\mathbf v)\,\mathbf v - \mu\,\hat{\mathbf R} \quad\Longrightarrow\quad \mu\,\dot{\mathbf e} = 2(\mathbf v\cdot\mathbf f)\,\mathbf r - (\mathbf r\cdot\mathbf f)\,\mathbf v - (\mathbf r\cdot\mathbf v)\,\mathbf f \]

In RSW components, using r v_S = h, the radial force drops out of the radial component:

\[ \mu\,\dot{\mathbf e} = 2h f_S\,\hat{\mathbf R} - \big(h f_R + r v_R f_S\big)\,\hat{\mathbf S} - r v_R f_W\,\hat{\mathbf W} \]

On the element side, 𝐞 = e P̂ and the perifocal triad turns with the Euler-angle rates. Projecting on the triad separates the pieces:

\[ \dot{\mathbf e} = \dot e\,\hat{\mathbf P} + e\,\boldsymbol\omega_{\mathrm{PQW}}\times\hat{\mathbf P}, \qquad \boldsymbol\omega_{\mathrm{PQW}} = \dot\Omega\,\hat{\mathbf K} + \dot i\,\hat{\mathbf n} + \dot\omega\,\hat{\mathbf W} \]

\[ \hat{\mathbf P}\cdot\dot{\mathbf e} = \dot e, \qquad \hat{\mathbf Q}\cdot\dot{\mathbf e} = e\,\big(\dot\omega + \cos i\;\dot\Omega\big), \qquad \hat{\mathbf W}\cdot\dot{\mathbf e} = -e\,\big(\sin i\cos\omega\;\dot\Omega - \sin\omega\;\dot i\big) \]

Now project the component form on P̂ = cos ν R̂ − sin ν Ŝ and Q̂ = sin ν R̂ + cos ν Ŝ. The P̂ projection reproduces the unsimplified ė of the eccentricity section term by term. The Q̂ projection needs (I1) once:

\[ 2h\sin\nu - r v_R\cos\nu = \frac{\mu\sin\nu}{h}\,\big(2p - re\cos\nu\big) = \frac{\mu\sin\nu}{h}\,(p+r) \]

\[ e\,\big(\dot\omega + \cos i\;\dot\Omega\big) = \frac{1}{\mu}\Big[-h\cos\nu\,f_R + \frac{\mu\sin\nu}{h}\,(p+r)\,f_S\Big] = \frac{1}{h}\Big[-p\cos\nu\,f_R + (p+r)\sin\nu\,f_S\Big] \]

That is the apse-line rate found above. The Ŵ projection tests the plane equations: inserting i̇ and Ω̇ collapses the bracket to sin(u − ω) = sin ν,

\[ -e\,\big(\sin i\cos\omega\;\dot\Omega - \sin\omega\;\dot i\big) = -\frac{e\,r f_W}{h}\,\big(\sin u\cos\omega - \cos u\sin\omega\big) = -\frac{re\sin\nu}{h}\,f_W = -\frac{r v_R}{\mu}\,f_W \]

which matches the Ŵ component above. Three projections, three agreements: the eccentricity vector confirms ė, the in-plane part of ω̇, and the i and Ω equations together.

9 Mean anomaly: M

The mean anomaly keeps its Keplerian drift n and responds only to in-plane forces, because it depends on the state only through ν and e.

9.1 M as a function of ν and e

Kepler’s equation and the eccentric anomaly make M a function of ν and e alone. The chain rule splits its rate into the Keplerian drift and two perturbation terms:

\[ M = E - e\sin E, \qquad E = E(\nu, e) \quad\Longrightarrow\quad \dot M = \mathcal K M + \frac{\partial M}{\partial\nu}\bigg|_{e}\mathcal D\nu + \frac{\partial M}{\partial e}\bigg|_{\nu}\mathcal D e, \qquad \mathcal K M = n \]

The partial derivatives of E follow from the expressions for cos E and sin E. Then dM = (1 − e cos E) dE = (r/a) dE converts them:

\[ \frac{\partial E}{\partial\nu}\bigg|_{e} = \frac{\sqrt{1-e^{2}}}{1+e\cos\nu}, \qquad \frac{\partial E}{\partial e}\bigg|_{\nu} = -\frac{\sin\nu}{\sqrt{1-e^{2}}\,(1+e\cos\nu)} \]

\[ \frac{\partial M}{\partial\nu}\bigg|_{e} = \frac{r^{2}}{ab}, \qquad \frac{\partial M}{\partial e}\bigg|_{\nu} = \frac{r}{a}\,\frac{\partial E}{\partial e}\bigg|_{\nu} - \sin E = -\sqrt{1-e^{2}}\;\frac{r\,(p+r)}{p^{2}}\,\sin\nu \]

The first partial gives a quick check. It maps the Keplerian part of ν onto exactly the mean motion, (r²/ab)(h/r²) = h/ab = n.

9.2 Assembling the rate

Substitute 𝒟ν and 𝒟e from the eccentricity section and collect the coefficients of f_R and f_S. Both collapse after replacing e r cos ν by p − r, as (I1) allows, and a by p/(1 − e²):

\[ \begin{aligned} f_R:&\quad \frac{r^{2}\cos\nu}{a} - \frac{er\,(p+r)\sin^{2}\nu}{p} = p\cos\nu - 2re \\[4pt] f_S:&\quad \frac{r^{2}}{a} + \frac{er\,(p+r)\cos\nu + e^{2}r^{2}}{p} = p \end{aligned} \]

\[ \boxed{\;\dot M = n + \frac{b}{ahe}\Big[(p\cos\nu - 2re)\,f_R - (p+r)\sin\nu\,f_S\Big]\;} \]

Splitting off the apse-line rate from the previous section gives a more revealing form:

\[ \dot M = n - \frac{2r}{\sqrt{\mu a}}\,f_R - \sqrt{1-e^{2}}\,\big(\dot\omega + \cos i\;\dot\Omega\big) \]

Apart from one radial term, M falls back as the apse line advances, scaled by √(1 − e²). The 1/e singularity is inherited entirely from the apse line, so the sum M + ω stays regular as e → 0.

9.3 The mean-anomaly-at-epoch subtlety

The boxed rate is for M itself, which is what a numerical integrator should carry. Many analytical theories use the mean anomaly at epoch instead, and there the varying mean motion leaves an explicit time in the equation:

\[ M = M_0 + n\,(t - t_0) \quad\Longrightarrow\quad \dot M_0 = (\dot M - n) - \dot n\,(t - t_0) = (\dot M - n) + \frac{3n}{2a}\,(t - t_0)\,\dot a \]

Even a purely periodic ȧ then drives terms in M₀ that grow in proportion to t − t₀. These mixed secular terms are an artifact of the element choice, not of the physics.

The standard remedy integrates the mean motion separately and keeps only the perturbative part in the new element σ:

\[ M = \int_{t_0}^{t} n\,dt' + \sigma \quad\Longrightarrow\quad \dot\sigma = \dot M - n = \frac{b}{ahe}\Big[(p\cos\nu - 2re)\,f_R - (p+r)\sin\nu\,f_S\Big] \]

The integral of n is carried alongside a, and σ obeys a Gauss equation with no explicit time. The time of periapsis passage τ, defined by M = n(t − τ), has the same ṅ(t − τ) problem.

10 Summary: the complete equations in RSW form

Six equations give the rates of the osculating elements from the radial, along-track and cross-track force components. Only M carries a Keplerian term.

\[ \begin{aligned} \dot a &= \frac{2a^{2}}{h}\Big(e\sin\nu\,f_R + \frac{p}{r}\,f_S\Big) \\[6pt] \dot e &= \frac{1}{h}\Big\{p\sin\nu\,f_R + \big[(p+r)\cos\nu + re\big]\,f_S\Big\} \\[6pt] \dot i &= \frac{r\cos u}{h}\,f_W \\[6pt] \dot\Omega &= \frac{r\sin u}{h\sin i}\,f_W \\[6pt] \dot\omega &= \frac{1}{he}\Big[-p\cos\nu\,f_R + (p+r)\sin\nu\,f_S\Big] - \frac{r\sin u\cos i}{h\sin i}\,f_W \\[6pt] \dot M &= n + \frac{b}{ahe}\Big[(p\cos\nu - 2re)\,f_R - (p+r)\sin\nu\,f_S\Big] \end{aligned} \]

Here r = p/(1 + e cos ν), u = ω + ν, b = a√(1 − e²), h = √(μp) and n = √(μ/a³), all evaluated on the osculating conic. The auxiliary quantities derived along the way are:

\[ \begin{aligned} \dot h &= r f_S, \qquad \dot p = 2r\sqrt{\frac{p}{\mu}}\,f_S \\[6pt] \dot\nu &= \frac{h}{r^{2}} + \frac{1}{he}\Big[p\cos\nu\,f_R - (p+r)\sin\nu\,f_S\Big] \\[6pt] \dot u &= \frac{h}{r^{2}} - \frac{r\sin u\cos i}{h\sin i}\,f_W \end{aligned} \]

The argument of latitude u stays defined on a circular orbit, so it replaces ν or M as the fast angle when e is near zero.

10.1 Matrix form

The same equations as one linear map from the force to the element rates:

\[ \frac{d}{dt}\begin{pmatrix} a\\ e\\ i\\ \Omega\\ \omega\\ M \end{pmatrix} = \begin{pmatrix} 0\\ 0\\ 0\\ 0\\ 0\\ n \end{pmatrix} + \underbrace{\begin{pmatrix} \dfrac{2a^{2}e\sin\nu}{h} & \dfrac{2a^{2}p}{hr} & 0 \\[10pt] \dfrac{p\sin\nu}{h} & \dfrac{(p+r)\cos\nu + re}{h} & 0 \\[10pt] 0 & 0 & \dfrac{r\cos u}{h} \\[10pt] 0 & 0 & \dfrac{r\sin u}{h\sin i} \\[10pt] -\dfrac{p\cos\nu}{he} & \dfrac{(p+r)\sin\nu}{he} & -\dfrac{r\sin u\cos i}{h\sin i} \\[10pt] \dfrac{b\,(p\cos\nu - 2re)}{ahe} & -\dfrac{b\,(p+r)\sin\nu}{ahe} & 0 \end{pmatrix}}_{\textstyle \mathsf B(a,\,e,\,i,\,\omega,\,M)} \begin{pmatrix} f_R\\ f_S\\ f_W \end{pmatrix} \]

The matrix depends on a, e, i, ω and M through r, ν and u, but never on Ω. That is the Kepler problem’s symmetry about the polar axis; the force itself may still depend on Ω.

By the fundamental rule, B is the velocity gradient of the elements expressed in RSW axes. A small impulse Δv, resolved in RSW, therefore changes the elements by Δα ≈ B Δv.

11 Tangential–normal (TNW) form

Resolving the force along and across the velocity puts the entire change of the semi-major axis into the tangential component. This is the natural frame for thrust and drag.

11.1 The frame

T̂ points along the velocity and Ŵ along the orbit normal. N̂ = Ŵ × T̂ lies in the plane and points to the inside of the path, toward the center of curvature. With the flight-path angle γ defined by v_R = v sin γ and v_S = v cos γ,

\[ \hat{\mathbf T} = \frac{\mathbf v}{v} = \sin\gamma\,\hat{\mathbf R} + \cos\gamma\,\hat{\mathbf S}, \qquad \hat{\mathbf N} = \hat{\mathbf W}\times\hat{\mathbf T} = -\cos\gamma\,\hat{\mathbf R} + \sin\gamma\,\hat{\mathbf S} \]

\[ \sin\gamma = \frac{\mu e\sin\nu}{hv}, \qquad \cos\gamma = \frac{h}{rv}, \qquad v^{2} = \mu\Big(\frac{2}{r} - \frac{1}{a}\Big) \]

The force components transform by the same rotation:

\[ f_R = \sin\gamma\,f_T - \cos\gamma\,f_N, \qquad f_S = \cos\gamma\,f_T + \sin\gamma\,f_N \]

11.2 Converting the equations

The semi-major axis needs no algebra. The energy result ȧ = (2a²/μ) v·f already depends only on v·f = v f_T.

For e, substitute and use μp/h² = 1 and p/r = 1 + e cos ν. The f_N coefficient also needs (I1) in the form e r (p + r) cos ν = p² − r²:

\[ \begin{aligned} f_T:&\quad \frac{p\sin\nu\sin\gamma + \big[(p+r)\cos\nu + re\big]\cos\gamma}{h} = \frac{e\sin^{2}\nu + (2 + e\cos\nu)\cos\nu + e}{v} = \frac{2\,(e+\cos\nu)}{v} \\[8pt] f_N:&\quad \frac{-p\sin\nu\cos\gamma + \big[(p+r)\cos\nu + re\big]\sin\gamma}{h} = \frac{\sin\nu}{prv}\Big[-p^{2} + er\,(p+r)\cos\nu + e^{2}r^{2}\Big] = -\frac{r}{av}\,\sin\nu \end{aligned} \]

The ω and M coefficients reduce the same way:

\[ \begin{aligned} \omega,\ f_T:&\quad \frac{-p\cos\nu\sin\gamma + (p+r)\sin\nu\cos\gamma}{he} = \frac{2\sin\nu}{ev} \\[8pt] \omega,\ f_N:&\quad \frac{p\cos\nu\cos\gamma + (p+r)\sin\nu\sin\gamma}{he} = \frac{1}{ev}\Big(2e + \frac{r}{a}\cos\nu\Big) \\[8pt] M,\ f_T:&\quad \frac{b}{ahe}\Big[(p\cos\nu - 2re)\sin\gamma - (p+r)\sin\nu\cos\gamma\Big] = -\frac{2b}{eav}\Big(1 + \frac{e^{2}r}{p}\Big)\sin\nu \\[8pt] M,\ f_N:&\quad \frac{b}{ahe}\Big[-(p\cos\nu - 2re)\cos\gamma - (p+r)\sin\nu\sin\gamma\Big] = -\frac{b}{eav}\,\frac{r}{a}\cos\nu \end{aligned} \]

11.3 The TNW equations

The cross-track terms are untouched, since Ŵ is common to both frames:

\[ \begin{aligned} \dot a &= \frac{2a^{2}v}{\mu}\,f_T \\[6pt] \dot e &= \frac{1}{v}\Big[2\,(e+\cos\nu)\,f_T - \frac{r}{a}\,\sin\nu\,f_N\Big] \\[6pt] \dot i &= \frac{r\cos u}{h}\,f_W, \qquad \dot\Omega = \frac{r\sin u}{h\sin i}\,f_W \\[6pt] \dot\omega &= \frac{1}{ev}\Big[2\sin\nu\,f_T + \Big(2e + \frac{r}{a}\cos\nu\Big)f_N\Big] - \frac{r\sin u\cos i}{h\sin i}\,f_W \\[6pt] \dot M &= n - \frac{b}{eav}\Big[2\Big(1 + \frac{e^{2}r}{p}\Big)\sin\nu\,f_T + \frac{r}{a}\cos\nu\,f_N\Big] \end{aligned} \]

Some texts take the in-plane normal outward, T̂ × Ŵ. Every f_N coefficient then changes sign.

11.4 What the TNW form shows

Only f_T changes a, and at a rate proportional to the speed. A given tangential acceleration therefore raises a fastest at periapsis.

At an apsis, sin ν = 0. A tangential force there changes a and e but does not turn the apse line. A normal force there leaves a and e alone and only turns the apse line.

12 Structure, physical interpretation and singularities

The equations split cleanly by force direction, and each singularity sits exactly where an element’s geometric definition fails.

12.1 Which force moves which element

Element Radial f_R Along-track f_S Cross-track f_W
a yes, ∝ e sin ν yes, ∝ p/r no
e yes, ∝ sin ν yes, ∝ (p + r) cos ν + re no
h, p no yes, ∝ r no
i no no yes, ∝ r cos u
Ω no no yes, ∝ r sin u / sin i
ω yes, ∝ −cos ν / e yes, ∝ (p + r) sin ν / e yes, −cos i times the Ω coefficient
M yes, ∝ p cos ν − 2re yes, ∝ −(p + r) sin ν no

Three facts organize the table. The semi-major axis responds only to the power v·f, and the size of the angular momentum only to the torque r f_S.

A radial force therefore trades a against e at fixed p = a(1 − e²). Only f_W moves the plane, and it reaches ω solely because the node, the origin of ω, moves.

A radial impulse at periapsis shows the in-plane coupling most cleanly:

\[ \nu = 0:\qquad \Delta a = 0, \qquad \Delta e = 0, \qquad \Delta\omega = -\frac{p}{he}\,\Delta v_R, \qquad \Delta\nu = +\frac{p}{he}\,\Delta v_R \]

The conic keeps its size and shape and turns backward. The body, with its position unchanged, now sits just past its new periapsis.

12.3 Singularities

Condition Equations that blow up Cause Regular replacement
e → 0 ω, ν, M, through 1/e periapsis undefined u or mean longitude as the fast angle; e cos ϖ, e sin ϖ
i → 0 Ω, ω, through 1/sin i node undefined ϖ = Ω + ω; tan(i/2) cos Ω, tan(i/2) sin Ω
e → 1 a, M the conic becomes a parabola p in place of a; true longitude

The fixes for i → 0 fail at i = 180°, where a retrograde variant is used. In their p, r, ν forms, the equations for a, e, i, Ω, ω and ν also hold for hyperbolic orbits, with a < 0.

12.4 Regular combinations

Sums of angles that do not need the missing reference direction stay finite. The longitude of periapsis ϖ = Ω + ω is regular at i = 0, because 1 − cos i = sin i tan(i/2):

\[ \dot\varpi = \dot\Omega + \dot\omega = \frac{1}{he}\Big[-p\cos\nu\,f_R + (p+r)\sin\nu\,f_S\Big] + \frac{r\sin u}{h}\,\tan\frac{i}{2}\;f_W \]

The mean longitude λ = M + ω + Ω is regular at e = 0 as well. In the compact Ṁ, the apse-line rate is multiplied by 1 − √(1 − e²) ≈ e²/2, which cancels its 1/e:

\[ \dot\lambda = n - \frac{2r}{\sqrt{\mu a}}\,f_R + \Big(1-\sqrt{1-e^{2}}\Big)\big(\dot\omega + \cos i\;\dot\Omega\big) + \frac{r\sin u}{h}\,\tan\frac{i}{2}\;f_W \]

A fully regular set follows by applying 𝒟 to state functions with no singular geometry. The modified equinoctial elements of Walker, Ireland and Owens (1985) are the common choice. Here h and k are element names, not the angular momentum:

\[ p, \qquad f = e\cos(\omega+\Omega), \qquad g = e\sin(\omega+\Omega), \qquad h = \tan\tfrac{i}{2}\cos\Omega, \qquad k = \tan\tfrac{i}{2}\sin\Omega, \qquad L = \Omega + \omega + \nu \]

The results above convert directly. With w = 1 + f cos L + g sin L = p/r, the p equation is already in equinoctial form, and the true longitude follows from L̇ = Ω̇ + u̇:

\[ \dot p = \frac{2p}{w}\sqrt{\frac{p}{\mu}}\;f_S, \qquad \dot L = \frac{\sqrt{\mu p}}{r^{2}} + (1-\cos i)\,\dot\Omega = \sqrt{\mu p}\,\Big(\frac{w}{p}\Big)^{2} + \frac{1}{w}\sqrt{\frac{p}{\mu}}\,\big(h\sin L - k\cos L\big)\,f_W \]

The full equinoctial set is singular only at i = 180°.

13 Worked check: J2 nodal regression by orbit averaging

Averaging the node equation over one orbit recovers the classical J2 regression rate, −(3/2) n J2 (R_E/p)² cos i.

13.1 The J2 force in RSW components

The J2 disturbing function depends on position only through r and the height above the equator, z = r sin i sin u. Here R_E is the equatorial radius:

\[ \mathcal R = -\frac{\mu J_2 R_E^{2}}{2}\left(\frac{3z^{2}}{r^{5}} - \frac{1}{r^{3}}\right), \qquad z = \hat{\mathbf K}\cdot\mathbf r = r\sin i\sin u \]

Its gradient therefore has only a radial part and a polar-axis part. The K̂ projections from the notation section resolve it into RSW:

\[ \mathbf f = \nabla\mathcal R = \frac{\partial\mathcal R}{\partial r}\bigg|_{z}\hat{\mathbf R} + \frac{\partial\mathcal R}{\partial z}\bigg|_{r}\hat{\mathbf K}, \qquad \frac{\partial\mathcal R}{\partial z}\bigg|_{r} = -\frac{3\mu J_2 R_E^{2}\,z}{r^{5}} \]

\[ \begin{aligned} f_R &= \frac{3\mu J_2 R_E^{2}}{2r^{4}}\,\big(3\sin^{2} i\sin^{2} u - 1\big) \\[6pt] f_S &= -\frac{3\mu J_2 R_E^{2}}{r^{4}}\,\sin^{2} i\,\sin u\cos u \\[6pt] f_W &= -\frac{3\mu J_2 R_E^{2}}{r^{4}}\,\sin i\cos i\,\sin u \end{aligned} \]

13.2 Instantaneous and averaged rates

Inserting f_W into the node and inclination equations gives the instantaneous rates:

\[ \dot\Omega = -\frac{3\mu J_2 R_E^{2}\cos i}{h\,r^{3}}\,\sin^{2} u, \qquad \dot i = -\frac{3\mu J_2 R_E^{2}\sin i\cos i}{h\,r^{3}}\,\sin u\cos u \]

To first order, average over one osculating orbit with the elements frozen. Time converts to true anomaly through dt = r² dν / h, and h² = μp:

\[ \langle\dot\Omega\rangle = \frac{n}{2\pi}\int_{0}^{2\pi}\dot\Omega\,\frac{r^{2}}{h}\,d\nu = -\frac{3nJ_2R_E^{2}\cos i}{2\pi\,p^{2}}\int_{0}^{2\pi}\big(1+e\cos\nu\big)\sin^{2}(\omega+\nu)\,d\nu \]

The integral is π, because e cos ν sin²(ω + ν) has no constant part. The inclination integrand, (1 + e cos ν) sin u cos u, has none at all:

\[ \boxed{\;\langle\dot\Omega\rangle = -\frac{3}{2}\,n J_2\left(\frac{R_E}{p}\right)^{2}\cos i\;} \qquad \langle\dot i\rangle = 0 \]

13.3 Numbers

Prograde orbits regress and retrograde orbits advance. Take a circular orbit at a = 7000 km and i = 98°, with μ = 398600.4418 km³/s², R_E = 6378.137 km and J2 = 1.08262668 × 10⁻³:

\[ \langle\dot\Omega\rangle = +1.001325^{\circ}/\text{day} \]

That is close to the Sun-synchronous rate, 360° per 365.2422 days or 0.985647°/day. This orbit reaches it at i = 97.874°.

The same averaging applied to the ω equation gives the familiar apsidal rate, which vanishes at the critical inclination of 63.4°:

\[ \langle\dot\omega\rangle = \frac{3}{4}\,n J_2\left(\frac{R_E}{p}\right)^{2}\big(5\cos^{2} i - 1\big) \]

A numerical average of the instantaneous equations reproduces both formulas to six decimals, as the next section shows.

14 Numerical verification of every equation

Every equation in this document passes three independent numerical tests: finite-difference gradients, an end-to-end integration against Cartesian dynamics, and the J2 averages. The checks used Python with NumPy and SciPy, and the constants of the J2 section.

14.1 Method: gradients by finite differences

The fundamental rule turns verification into a finite-difference exercise. The perturbative part of each rate is the element change per unit velocity kick at fixed position:

\[ \frac{\alpha(\mathbf r,\ \mathbf v + \varepsilon\hat{\mathbf d}) - \alpha(\mathbf r,\ \mathbf v - \varepsilon\hat{\mathbf d})}{2\varepsilon} = \frac{\partial\alpha}{\partial\mathbf v}\cdot\hat{\mathbf d} + O(\varepsilon^{2}) \]

  1. Draw a random orbit: a from 7000 to 42000 km, e from 0.02 to 0.9, i from 0.1 rad to π − 0.1 rad, and uniform Ω, ω, ν.
  2. Kick the velocity by ±ε along each RSW or TNW axis, with ε = 10⁻⁶ v, and convert both states to elements with a standard state-to-element routine.
  3. Form the central difference, wrapping angle differences into (−π, π], and compare it with the matching column of B.
  4. Repeat for 500 orbits and keep the worst relative error of each element’s row of B.

14.2 Results

The six element rates agree with the finite differences to 10⁻⁷ or better in both frames:

Rate RSW form TNW form
ȧ 1.3 × 10⁻⁹ 1.4 × 10⁻⁹
ė 2.2 × 10⁻⁹ 2.3 × 10⁻⁹
i̇ 4.4 × 10⁻⁸ 4.4 × 10⁻⁸
Ω̇ 6.7 × 10⁻⁸ 6.7 × 10⁻⁸
ω̇ 4.7 × 10⁻⁹ 4.7 × 10⁻⁹
Ṁ 4.6 × 10⁻⁹ 4.6 × 10⁻⁹

The auxiliary rates, the alternate forms and the J2 components pass as well:

Quantity Compared with Worst relative error
ν̇ finite difference 4.7 × 10⁻⁹
u̇ finite difference 9.7 × 10⁻⁸
ḣ and ṗ finite difference 3.1 × 10⁻¹⁰
eccentricity-vector rate finite difference 1.8 × 10⁻⁷
ė, cos E form ė, p–r–ν form 5.6 × 10⁻¹⁶
Ṁ, compact form Ṁ, boxed form 5.9 × 10⁻¹⁵
J2 components f_R, f_S, f_W Cartesian J2 acceleration projected on RSW, 200 states 1.9 × 10⁻¹⁵

All gradient errors sit at the level expected of a central difference with this step. The alternate forms agree to rounding, as identities should.

14.3 End-to-end integration

The gradient test cannot see the Keplerian parts, so a second test integrates the complete equations. It propagates one orbit for five periods two ways and compares the final positions:

  1. Start from a = 9000 km, e = 0.3, i = 50°, Ω = 30°, ω = 40°, ν = 10°.
  2. Apply J2 plus a time-varying thrust with RSW components 2 × 10⁻⁶ cos(t / 900 s), 3 × 10⁻⁶ and 2.5 × 10⁻⁶ sin(t / 1300 s) km/s².
  3. Integrate the Cartesian equations, and separately the Gauss equations, with DOP853 at rtol = atol = 10⁻¹².
  4. For the Gauss run, carry either M, solving Kepler’s equation at each step, or ν as the sixth element.

The final positions differ by 0.129 mm with M and 0.494 mm with ν, at a radius of 7875 km. The relative agreement is 2 × 10⁻¹¹ and 6 × 10⁻¹¹, at the integrator’s tolerance.

14.4 J2 averages

Averaging the instantaneous Gauss rates numerically over one orbit reproduces the closed forms of the J2 section to all six decimals shown:

a (km) e i ω ⟨Ω̇⟩ (°/day) ⟨ω̇⟩ (°/day)
7000 0 98° — +1.001325 undefined at e = 0
9000 0.3 50° 40° −2.317389 +1.921363
12000 0.4 63° 25° −0.701805 +0.023603
26000 0.7 30° 250° −0.242575 +0.385140

The near-zero apsidal rate at 63° is the critical inclination at work. The averaged inclination rate stayed below 10⁻¹⁵ °/day in every case.