Celestial Programming : Two Star Celestial Navigation Fix

A two star fix, using only the measured altitude of two stars, usually produces two possible location for an observer. Usually one of the two points can be disregarded based on the observer's knowledge of their approximate location, with one point being near by, and the other being thousands of miles away. It's also possible to perform a three star fix, which is just doing the two star fix twice, but with a third star to disambiguate which of the two points is the correct one.

\(\lambda_i, \delta_i, \zeta_i\) are the geographic position longitude, declination, and zenith distance of the two observations. $$ \begin{align*} a_i &= \cos \lambda_i \cos \delta_i \\ b_i &= \sin \lambda_i \cos \delta_i \\ c_i &= \sin \delta_i \\ p_i &= \cos \zeta_i \\ \\ A &= a_1 b_2 - a_2 b_1 \\ B &= b_2 c_1 - b_1 c_2 \\ C &= b_2 p_1 - b_1 p_2 \\ D &= a_1 c_2 - a_2 c_1 \\ E &= b_1 c_2 - b_2 c_1 \\ F &= c_2 p_1 - c_1 p_2 \\ \\ 0 & = x^2\left(1+\frac{A^2}{B^2} + \frac{D^2}{E^2}\right) + x\left(\frac{-2AC}{B^2} - \frac{2FD}{E^2}\right) + \frac{C^2}{B^2} + \frac{F^2}{E^2} - 1 \\ \end{align*} $$ The above equation is a quadradic equation, and is solved using the quadradic formula, which yields two solutions: \(x_1, x_2\). $$ \begin{align*} a &=1+\frac{A^2}{B^2} + \frac{D^2}{E^2} \\ \\ b &=\frac{-2AC}{B^2} - \frac{2FD}{E^2} \\ \\ c &=\frac{C^2}{B^2} + \frac{F^2}{E^2} - 1 \\ \\ x_1 &=\frac{-b+\sqrt{b^2 - 4ac}}{2a} \\ \\ x_2 &=\frac{-b-\sqrt{b^2 - 4ac}}{2a} \\ \end{align*} $$ Now use the values of \(x_i\) to compute the y and z values. $$ \begin{align*} y_i &=\frac{F-Dx_i}{E} \\ \\ z_i &=\frac{C-Ax_i}{B} \\ \end{align*} $$ \(x_i, y_i, z_i\) are the rectangular coordinates of the matching points. Convert them to polar coordinates to get latitude and longitude.