Hi guys, I pondered the question of calculating Kepler orbits recently and just for fun I implemented the algorithms implemented in "Astronomy On The Personal Computer" by Montenbruck & Pfleger. This is the code that I got by translating the code in the book: https://anonscm.debian.org/cgit/users/tomasz/orbits.git/. It is around 90% identical to what is in the book.
The tricky cases are:
1) (ecc ~ 1.0) and (e (ecc. anomaly) ~ 0 mod 2*pi) (i.e., the original case);
they cause problems because the Newton's method divides by
(1.0 - ecc * cos(e)) (elliptic case) or (ecc*cosh(e) - 1.0) (hyperbolic
case) which are ~ 0
2) ecc = 1.0; in this case the orbit is parabolic and the formula
does not work at all; in the code, however, there is a special case
for parabolic orbits which works great (it also applies a bit to ecc ~
1.0)
Stellarium takes care of the second case explicitely, but if ecc ~ 1.0,
it simply calls into elliptic or hyperbolic code, which although uses
a bit different algorithms, still suffers from the 1st tricky case.
There are a few ways to solve the problem:
1) approximate the 1st case with the second one when e ~ 0 mod 2*pi;
this is what the original code tries to do, sort of; the difficulty is to
find the cutoff point (the book chose quite arbitrary 5 degrees)
2) clamp ecc ~ 1.0 to 1.0 for *any* e and use parabolic formula (the
motivation
for that is that if ecc ~ 1.0 than the object has a huge period, so
there is a low chance that using the precise value of ecc will be
accurate;
moreover, it seems to me that eccentrities of such objects are very
approximate anyway)
3) use a better method for finding roots, a one that does not divide by zero
(bisection?) and apply it to the tricky case; there is an implicit
"trickiness" in the 1st case since the root is very close to the extremum
(try apply Newton's method to x^2 = 0!)
Any ideas or comments?
Cheers,
Tomasz
signature.asc
Description: Digital signature
