div by 0 in FindOrbitalElements
in that function, when hz >> hy+hx
then sqrt(hx*hx + hy*hy + hz*hz) \approx hz
, so inc = acos(hz/h) \approx acos(1) \approx 0
As a result, a few lines down, in z/sin(inc)
we get a div by 0.
So in
232 u = atan2 ( z/sin(inc) , x*cos(capom) + y*sin(capom));
which is, in this case atan2(\infty, 1)
, do we want to detect the case and set u
to \pi/2
?