What's wrong with this algorithm to determine Kepler parameters from orbital vectors?

Question

I am writing computer code which instantiates a Kepler orbit given cartesian position and velocity vectors. I am working from this reference document.

It says:

Calculate orbital momentum vector

$\mathbf{h} = \mathbf{r} \times \mathbf{\dot r}$

Determine the vector $\mathbf n$ pointing towards the ascending node

$\mathbf n = (0,0,1)^T \times \mathbf{h} = (-h_y, h_x, 0)^T$

Obtain the longitude of the ascending node

$\Omega = \mathrm{arccos}(\frac{n_x}{||\mathbf{n}||})$

Suppose we have $\mathbf{r} = (X,0,0)$ and $\mathbf{\dot{r}}=(0,Y,0)$, we will end up with $\mathbf{h}=(0,0,Z)$, but this causes $\mathbf n$ to collapse to zero, and then $\Omega$ and other parameters can't be calculated.

Is the document wrong, or am I reading it wrong, or do I just need to detect and handle certain special cases? If so, what are these special cases, exhaustively?

Answer 1

Nothing is wrong with that algorithm; everything is wrong with the longitude of the ascending node. It is undefined when the reference plane and the orbital plane are the same. See the wikipedia page for details. Even when they are close, and the algorithm gives you finite answers, you should expect bad numerical accuracy. While these sorts of singular parameters are reasonable for observational astronomy because the corresponding observations are also singular, for computational astronomy you might prefer a non-singular representation.

Answer 2

Exactly as Mike said, it is undefined in the case of non-inclined orbits. Other way of looking at it is that all the directions in the plane of reference / orbital plane are of equal importance and therefore you can choose any of them. In case you are struggling with the implementation here's something which you may find useful. The code works for all conic intersections.