Wir berechnen im Folgenden die Entfernung zweier Punkte auf der Erdoberfläche und die einzuschlagende Richtung, um vom ersten zum zweiten zu gelangen. Dabei bezeichnen wir den Ausgangspunkt mit seiner geografischen Länge L und Breite B und das Ziel durch die Längendifferenz l und Breitendifferenz b, d.h. der Zielpunkt hat die geografische Länge L+l und Breite B+b. Ganz offenbar können wir den Erdball um seine Achse drehen bis L = 0 ohne dass sich bei gegebenem l etwas an der einzuschlagenen Richtung oder Entfernung ändert; deswegen taucht L in den Formeln nicht mehr auf.
Wir legen nun ein cartesisches, also rechtwinkliges, Koordinatensystem durch die Erde mit dem Erdmittelpunkt als Ursprung, dem ersten Einheitsvektor von dort auf denjenigen Punkt des Äquators weisend, der auf dem Längengrad des Ausgangspunkts liegt, dem zweiten auf den um 90° östlich davon liegenden Punkt des Äquators, und dem dritten auf den Nordpol. In diesem System hat dann der Zielpunkt die Koordinaten <cos (B+b) cos l; cos (B+b) sin l; sin (B+b)>. Drehen wir nun die Kugel um die Achse <0; 1; 0> so weit, dass der Ausgangspunkt auf den Nordpol verschoben wird, dann ergibt sich die gesuchte Richtung und Entfernung als Richtung und Entfernung vom Nordpol. Das ist dann viel einfacher zu rechnen.
Die Spaltenvektoren der Drehmatrix sind die Koordinaten der Punkte, auf die die Enden der Einheitsvektoren (vom Ursprung aus gesehen) abgebildet werden. Die Koordinaten des Zielpunktes ergeben sich dann durch einfache Matrixmultiplikation.
Jetzt betrachten wir den Ausgangspunkt als konstant und interessieren uns für die Koordinaten des Zielpunktes in Abhängigkeit von l und b. Mit B sind auch dessen Sinus und Cosinus Konstante, denen wir Namen geben: S = sin B und C = cos B. Damit vereinfachen sich die Komponenten des der Zielvektors zu
x = S cos (B+b) cos l – C sin (B+b)
y = cos (B+b) sin l
z = C cos (B+b) cos l + S sin (B+b)
Die Formel für x hat eine numerische Tücke. Für l = 0 und b = 0 ergibt sich x = S cos B – C sin B = SC – CS = 0. Das ist zweifellos richtig, aber in der unmittelbaren Nähe von l = 0 und b = 0 erscheint eine Differenz zweier fast gleicher, von Null deutlich verschiedener Größen. Sind diese mit Rundungsfehlern behaftet, so bleiben nach der Subtraktion nur die Rundungsfehler übrig. Um diesen Effekt zu vermeiden, müssen wir ein wenig Trigonometrie treiben. Dazu führen wir ein paar Hilfsgrößen ein.
u = 1 – cos l
v = sin b
w = cos b
Mit diesen sieht die Anwendung der Additionstheoreme für Sinus und Cosinus recht übersichtlich aus.
sin (B+b) = sin B cos b + cos B sin b = S cos b + C sin b = Sw + Cv
cos (B+b) = cos B cos b – sin B sin b = C cos b – S sin b = Cw – Sv
und für x ergibt sich
x = S cos (B+b) cos l – C sin (B+b) =
= S (Cw – Sv) (1 – u) – C (Sw + Cv) =
= CSw – S²v – CSwu + S²vu – CSw – C²v =
= –v – CSwu + S²vu =
= (S²u – 1) v – CSwu
Das ist zwar immer noch eine Differenz, aber deren Operanden werden mit kleinem l ebenfalls klein, und mit ihnen eventuelle Rundungsfehler, falls diese proportional zur Größe selbst sind. Das war vorher anders. Die Formel für z enthält ebenfalls eine Differenz, aber z = 0 tritt nur bei großen Entfernungen auf, wo kleine Fehler keine Rolle mehr spielen. Allerdings muss man in der Umgebung von z = ±1 mit numerischer Instabilität rechnen; wir werden das berücksichtigen, indem wir dort die Entfernungen ohne Verwendung von z berechnen.