Preface

While I was in the process of starting work with a new company, I spent my time improving my knowledge and my tools.

During this I translated my astrodynamics library from Python to C++ to improve its performance. Lot of algorithms there are based on those from the book of Mr. David Vallado “Fundamentals of Astrodynamics and Applications”.

Some of the algorithms in the book may look tricky at first glance and work not as expected. One of such an algorithms is an “Algorithm #13: ECEFToLatLon“, based on work of Borkowski (1989). In a nutshell, this algorithm allows to calculate geodetic coordinates from Cartesian position vector without iterating (just directly). And when I tried to implement it as it is provided in the book, it didn’t work.

Knowing from other examples that some algorithms in “Fundamentals” require deeper study, I dug into original sources of the algorithm. The author’s name looked quite Slavic to me, so I expected to find some specific instructions in original work, which may have been lost during translation.

And I found it when I read the original Polish text.

GEOD

Originally, algorithm was named “An Algorithm of Replacement of Carthesian Coordinates with Ellipsoidal”¹, and was developed by Mr. Kazimierz M. Borkowski from Nicolaus Copernicus University in Toruń, Poland. The algorithm provides exact analytical solution of coordinate conversion, and looks pretty elegant. In the core of algorithm lays the calculation of fourth-degree polynomial to evaluate results without iteration.

The source code of the algorithm is presented as a FORTRAN procedure on the site of University of Copernicys.

The procedure is named GEOD, and, if you look at the picture illustrating the logic of the algorithm, it seems clear why it is named like this (the name resembles the word “geode”, which fits nicely with the crystal-like geometry shown in the diagram).

Notes about implementation

To make the algorithm work properly, you must:

1. Provide to it elements of vector $\vec{r}$, projected onto the plane defined by the cross product of $\vec{r}$ and unit vector $\hat{K}$ (from $IJK$CS).
2. The equatorial and polar components of $\vec{r}$ must be in meters (not in km like in convention used in “Fundamentals”). This is important, because conversion polinominal is not equivalent, and wrong scale will provide extremely incorrect output.
3. $\vec{r}$ must not be much shorter than radius of Earth. In this case algorithm will provide errors because of logic of conversion (see picture).

Output will be geodetic latitude (in radians) and height from reference ellipsoid of Earth (in meters). Longitude should be calculated with different algorithm.

1. All rights to the referenced books, sites, pictures and algorithms belong to their respective authors (see ULTRADYNE.SPACE LEGAL NOTICE AND TERMS OF USE). ↩︎