I am running an animation of a satellite in an elliptic orbit (defined by a parametric equation for $x$ and $y$ as a function of $t$) and want to make sure the spacecraft is traveling at the right speeds at different points in its orbit. That is, it should go slower at is apoapsis and much faster at its periapsis. I can easily calculate the tangential speed of the satellite using this equation:
$v=\sqrt{GM(\cfrac{2}{r}-\cfrac{1}{a})}$
How do I convert this to the angular speed of the satellite at this point?
I've done extensive research (hours and hours) but haven't found anything of value. The closest thing was this expression of Kepler's Third Law:
$\cfrac{dA}{dt}=\cfrac{1}{2}r^2\omega$
Since $\cfrac{dA}{dt}$ is a rate (area swept out per second) I rewrote this equation as
$\cfrac{A}{P}=\cfrac{1}{2}r^2\omega$
where $A$ is the area of the elliptic orbit (given by $A=\pi ab$ where $a$ and $b$ are the semi-major and semi-minor axes of the ellipse, respectively), and $P$ is the period of the elliptic orbit (given by $P=2 \pi \sqrt{\cfrac{a^3}{GM}}$). Solving this for $\omega$ yields:
$\omega=\cfrac{2A}{Pr^2}$
For each time step in my simulation I use the satellite's current position in this equation to compute $\omega$ and then use the result to update the current $\theta$. This updated $\theta$ is then plugged into the parametric equation mentioned above to get the satellite's $x$ and $y$ position.
I can't find my mistake anywhere and would really appreciate it if someone could point it out to me.
Best Answer
The formula $$ \dot{\theta} = \omega = \frac{2A}{Pr^2} $$ is correct; it can also be derived from the specific angular momentum $h$: $$ h = r^2\omega = \sqrt{GMa(1-e^2)} = b\sqrt{\frac{GM}{a}}, $$ with $e = \sqrt{(a^2-b^2)/a^2}$ the orbital eccentricity. However, this doesn't solve the Kepler problem, because both $\omega$ and $r$ depend on $t$ in a complicated way, which isn't specified by the above formula. In other words, the above formula gives you $\omega(r)$, but not $\omega(t)$ and $r(t)$.
Also, note that $\theta$ is the true anomaly, which is the angle between the direction of periapsis and the current position of the body, as seen from the main focal point (where the attracting body is). And $r$ is the distance between the current position and the focal point. If you want to use cartesian coordinates $(x,y)$, it is better to parametrize them using the eccentric anomaly $E$: $$ x = a\cos E,\qquad y = b\sin E. $$
So how to find $E(t)$? For this, we need to introduce another parameter, called the mean anomaly $M$. The mean anomaly increases linearly with time: $$ M(t) = \frac{2\pi}{P}t = \sqrt{\frac{GM}{a^3}}t. $$ From $M(t)$, we can calculate the eccentric anomaly $E(t)$ using Kepler's equation $$ M = E - e\sin E, $$ which you have to solve numerically. Once $E(t)$ is known, $(x(t),y(t))$ follow immediately.
For completeness, the true anomaly can be calculated from the eccentric anomaly: $$ \tan\frac{\theta}{2} = \sqrt{\frac{1+e}{1-e}}\tan\frac{E}{2}, $$ and the distance $r$ to the attracting body at the focal point is $$ r = \frac{a(1-e^2)}{1+e\cos\theta}. $$