Solution (source code)

= Solution

Treat $a,e$ as instantaneous <osculating orbital elements> of the fixed-central-parameter <Kepler orbit>. The perturbing <work> changes the <specific orbital energy> according to
$$
\dot\varepsilon=\boldsymbol v\cdot(\bar R\widehat{\boldsymbol r}+\bar T\widehat{\boldsymbol\theta})
 =\bar R\dot r+\bar T\frac hr,\qquad
\dot a=\frac{2a^2}{\mu}\dot\varepsilon.
$$
The <specific angular momentum> and conic equation give
$$
\dot r=\frac{\mu}{h}e\sin f,\qquad
\frac hr=\frac{\mu}{h}(1+e\cos f),\qquad h=\sqrt{\mu a(1-e^2)}.
$$
Substituting proves the <Gauss planetary equation for the semi-major axis>:
$$
\boxed{\dot a=2\sqrt{\frac{a^3}{\mu(1-e^2)}}
[\bar R e\sin f+\bar T(1+e\cos f)]}.
$$
Here $f$ denotes <true anomaly>, not the fractional <luminosity> used in Question 1. The equation follows from <work> alone; a normal perturbing <acceleration> has no instantaneous contribution because the <velocity> lies in the <orbital plane>. For slowly averaged evolution, the perturbation must be weak over one <orbital period>.