Differential Equation of orbit of a planet

Polar coordinates

The movement vector $\vec{OP}=\vec r$ in polar coordinates depends of $\theta$ which in turn depends of $t$, so $r=r(\theta(t))$.

But we will see that if $\theta$ is a convenient angle to describe the trajectory curve it is not straightforward to describe how $P$ evolves in time (or equivalently, to have a formula for $\theta(t)$). Instead we will introduce a new "angle" $E$ named eccentric anomaly and we will see that there are formulas to relate $E\leftrightarrow\theta$ and also to relate $E\leftrightarrow t$.

Thus giving a way via $E$ intermediate to express both $r$ and $\theta$ as direct functions of time $t$.

Central force movement

gravity law : $m\vec{\ddot r}=\vec F=-\frac{\mathcal GmM}{r^3}\vec r\quad$ where $m$ is planet mass and $M$ is sun mass.

moment law : $\vec\sigma_O=\vec{OP}\wedge m\vec v=\vec r\wedge m\vec{\dot r}$

When derivating we get $\frac{\partial\vec\sigma_O}{\partial t}=\vec{\dot r}\wedge m\vec{\dot r}+\vec r\wedge m\vec{\ddot r}=0+\vec r\wedge\vec F=\vec 0\quad$ since $F$ is a radial force.

Thus $\vec\sigma_O$ is a constant $\vec C$, we name it the area constant.

$\bbox[5px,border:1px solid] {\vec C=\vec r\wedge\vec{\dot r}}$

In the frame $(\vec i,\vec j,\vec k)$ defined by $\vec i=\frac{\vec r}{r},\ \vec j=\frac{\partial\vec i}{\partial\theta},\ \vec k=\vec i\wedge\vec j$.

Note: $\vec i=(\cos\theta,\sin\theta,0)$ so $(\vec i,\vec j,\vec k)$ is a well defined orthonormal frame.

Let's multiply by $\vec v=\vec{\dot r}$ and calculate the kinetic energy : $m\,\vec{\ddot r}.\vec{\dot r}=-\mathcal GmM\,\frac{\vec r.\vec{\dot r}}{r^3}\ $

integrating this gives $\frac12mv^2=+\mathcal GmM\frac 1r+cst$

We call the arbitrary integration constant $h=\frac{cst}{m}$ and also $\bbox[5px,border:1px solid]{\mu=\mathcal GM}$

$\vec{\dot r}=\frac{\partial(r\vec i)}{\partial t}=\dot r\,\vec i+r\dot\theta\,\vec j\quad$ from this we get $\vec C=r^2\dot\theta\,\vec k\ $ and $\ v^2={\dot r}^2+r^2{\dot\theta}^2$

$\bbox[5px,border:1px solid] {\begin{array}{cc} C=r^2\dot\theta\\ \frac12({\dot r}^2+r^2{\dot\theta}^2)-\frac{\mu}{r}=h \end{array}}$

Conic trajectory

This is where the manipulation is a bit artificial, we express $\dot r=\dot\theta\frac{dr}{d\theta}$ in terms of $\dot\theta=\frac{C}{r^2}$

giving $\frac12(\frac{C^2}{r^4}(\frac{dr}{d\theta})^2+\frac{C^2}{r^2})-\frac{\mu}{r}=h,\quad$and we do the change of variable $\frac1r=u+\frac{\mu}{C^2}$.

I do not detail the calculacy but this gives :

$\bbox[5px,border:1px solid] {(\frac{du}{d\theta})^2+u^2=\frac{\mu^2e^2}{C^4}}\quad$ where $\ \bbox[5px,border:1px solid]{e=\sqrt{1+\frac{2hC^2}{\mu^2}}}\ $ is called the eccentricity.

This integrates into $u=\frac{\mu e}{C^2}\cos(\theta-\theta_0)$ and we will introduce yet some new notations.

$\bbox[5px,border:1px solid] {\begin{array}{c} \nu=\theta-\theta_0 \\ p=\frac{C^2}{\mu} \\ \end{array}}\quad$ giving $\quad\bbox[5px,border:3px solid blue] {r=\frac{p}{1+e\cos \nu}}$

We can remark that :

• $e=0\quad$ this is a circle

• $h<0\Rightarrow 0<e<1\quad$ this is an ellipse

• $h=0\Rightarrow e=1\quad$ this is a parabola

• $h>0\Rightarrow e>1\quad$ this is an hyperbola

The case which interest us currently is the ellipse. Let's call $2a$ the length of the major axis.

We see that $r$ is minimum for $\nu=0$ and $r_{min}=a(1-e)=\frac p{1+e}$

$\bbox[5px,border:1px solid] {p=a(1-e^2)}\quad$ remark that we have also $a=\frac{\mu}{-2h}$, we will use it later.

Note: for Earth $e\simeq0.0167$ and $a\sim1\,UA\simeq149\,598\,000\ km$

Third law of Kepler

Let's introduce our new angles.

• $\nu$ is called the true anomaly

• $E$ is called the eccentric anomaly

• $M=n(t-t0)=2\pi(\frac{t-t0}{T})$ is called the mean anomaly, $n$ the frequency and $T$ the period.

With the new angle $E$ we have

$\begin{cases} x=r\cos\nu=a(\cos E-e) \\ y=r\sin\nu=b\sin E=a\sqrt{1-e^2}\sin E \\ \end{cases}$

The equation of the ellipse becomes $\bbox[5px,border:3px solid blue] {r=a(1-e\cos E)}$

We deduce also without much trouble that $\bbox[5px,border:3px solid green] {\tan(\frac{\nu}{2})=\sqrt\frac{1+e}{1-e}\tan(\frac E2)}$

The tedious part is establishing the Kepler's equation. I will not detail everything here, but basically you define $dt=\frac{r}{\sqrt{-2h}}d\tilde E$ an elementary area element.

Then by replacing $r,\dot r,\ddot r$ by the appropriate terms $r,r'=\frac{dr}{d\tilde E},r''=\frac{d^2r}{d\tilde E^2}$ in the equation $r\ddot r+{\dot r}^2=2h+\frac{\mu}{r}$ obtained by eliminating $\dot \theta$ from the movement equation.

You eventually arrive to $\bbox[5px,border:1px solid] {r''+r=a}$ which combined to $\sqrt{-2h}(t-t_0)=\int\limits_0^{\tilde E} r(\tilde E)d\tilde E$ proves to be compatible with the $E=\tilde E$ and $r(E)$.

The latter integral relation gives you $\sqrt{-2h}(t-t_0)=a(E-e\sin E)$

We define $\bbox[5px,border:3px solid green] {M=2\pi(\frac{t-t_0}{T})=E-e\sin E}$ and gain $Ma=\frac{2\pi a}{T}=\sqrt{-2h}=\sqrt\frac{\mu}{a}$

which is the famous $\bbox[5px,border:1px solid] {\frac{T^2}{a^3}=\frac{4\pi^2}{\mu}}\ $.

As promised we have now the two blue movement trajectories $r(\theta),r(E)$ expressed in two different frames, and the angles $\theta,E$ are linked to the time $t$ thanks to the green equations.