I'm working on a project where I'm trying to describe the orbits of the planets in the solar system using the polar equation of an ellipse.
Below is an equation I got from here. $$ r = \frac{b^{2}}{a - c \cos \theta} $$
Where $a$ is the semi-major axis, $b$ is the semi-minor axis, $c$ is the distance between the center of the orbit and a focus point of the orbit $(c^2 = a^2 - b^2)$. The sun is at one of these foci, so $c$ is the distance between the center of the orbit/ellipse and the sun, the origin of the graph. This form of the equation also assumes no displacement of the ellipse in the $y$ direction from the origin.
This equation does relate the radius $r$ of a point on the ellipse as a function of the angle $\theta$ it makes with the origin.
However, I am trying to look for an equation that models the angle $\theta$ as a function of time $t$. For example, if $T$ was the period of one full orbit, then after $T$ seconds, the change in angle should be $2\pi$ radians. In other words, where if $f(t) = \theta$, $f(T) = 2\pi$ or $0$.
Considering that I have the values for $a$, $b$, and $c$, the orbital period $T$, and the linear speed $v$ of each planet, and since I cannot use the radius directly (because it is not constant across the ellipse), how could I derive the polar angle $\theta$ as a function of time $t$ without $r$?
For the time dynamic the Newton gravity law additionally results in the second Kepler law $$ r^2\dotθ=L=const.\tag1 $$ that allows to solve the angle-time problem via numerical integration.
One can also try to solve it symbolically, or at least reduce it to a simpler non-ODE form. Eliminating first $r$ gives $$ L=\frac{b^4\dotθ}{(a-c\cosθ)^2}.\tag2 $$ For integrals of fractions in trigonometric functions the recommendation is to use the half-angle tangent $u=\tan\fracθ2$ with $\dot u=\frac12(1+u^2)\dotθ$ and $\cos θ=\frac{1-u^2}{1+u^2}$. The equation simplifies to $$ \frac{L}{b^4}=\frac{2(1+u^2)\dot u}{((a-c)+(a+c)u^2)^2}\tag3 $$ This form now strongly suggests to set (using again a half-angle in the tangent for symmetry) \begin{align} \tan\fracε2&=\sqrt{\frac{a+c}{a-c}}u\tag4 \\ \implies (a-c)+(a+c)u^2&=(a-c)\sec^2\fracε2\tag{5a}\\ 1+u^2&=\frac{(a+c)\cos^2\fracε2+(a-c)\sin^2\fracε2}{a+c}\sec^2\fracε2=\frac{a+c\cosε}{a+c}\sec^2\fracε2\tag{5b} \\ 1-u^2&=\frac{(a+c)\cos^2\fracε2-(a-c)\sin^2\fracε2}{a+c}\sec^2\fracε2=\frac{c+a\cosε}{a+c}\sec^2\fracε2\tag{5c} \end{align} With the derivative formulas $$\begin{align} \frac{\dotε}2\left(1+\frac{a+c}{a-c}u^2\right)&=\sqrt{\frac{a+c}{a-c}}\dot u\tag{6a} \\\text{ or }~~ \dotε((a-c)+(a+c)u^2)&=2b\dot u\tag{6b} \end{align}$$ one can insert and eliminate $u$ for $ε$, simplifying the equation to $$\begin{align} \frac{L}{b^3}&=\frac{(1+u^2) \dotε}{(a-c)+(a+c)u^2} =\frac{a+c\cosε}{a^2-c^2}\dotε\tag{7a} \\~\\ \frac{L}{b}&=(a+c\cosε)\dotε \tag{7b} \end{align}$$ Now with an arbitrary change of the momentum constant to $L=ab\omega$, integration results in the Kepler equation $$\tag8 ωt=ε+\frac{c}{a}\sinε, ~~~ \tan\fracθ2=\sqrt{\frac{a-c}{a+c}}\tan \fracε2 $$ The integration constant was left out which is correct if time $t=0$ corresponds to apoapsis with $θ=0$, where $r=\frac{b^2}{a-c}=a+c$ and thus also $ε=0$. In this context $ε$ is called* the eccentric anomaly (and $θ$ the true anomaly, $M=ωt$ the mean anomaly, there might have been something lost in the historic process, one would think that anomaly refers to the deviation from circular motion.)
For the Cartesian coordinates one gets from the half-angle formulas \begin{alignat}{1} x=r\cosθ&=\frac{b^2(1-u^2)}{(a-c)+(a+c)u^2}&=c+a\cosε\tag{9a}\\ y=r\sinθ&=\frac{b^2·2u}{(a-c)+(a+c)u^2}&=b\sinε\tag{9b} \end{alignat}
* Note that this deviates from the common definitions in that the anomalies are the angles at the gravity center relative to the pericenter position, not as here to the apocenter position.