Recently the work I have been doing requires a higher degree of accuracy than which the haversine method I was using could provide. This prompted me to implement a Python version of the Vincenty’s inverse formula. Unlike the Haversine method (which I posted about previously) of directly calculating the great-circle distance between two points on a perfectly spherical Earth, Vincenty’s formulae is an iterative method which more realistically assumes Earth as an oblate spheroid.

Before proceeding further, several values need to be defined:

\( a = 6378137.0 \) : length of semi-major axis of the ellipsoid aka radius of the equator in meters (WGS-84) \( f = 1 / 298.257223563 \) : flattening of the ellipsoid (WGS-84) \( b = \left ( 1 - f \right ) a \) : length of semi-minor axis of the ellipsoid (radius at the poles, 6356752.314245 meters in WGS-84) \( \phi_1, \phi_2 \) : latitude of the points \( U_1 = \arctan \left [ \left ( 1 - f \right ) \tan \phi_1 \right ] \) : reduced latitude (latitude on the auxiliary sphere) \( U_2 = \arctan \left [ \left ( 1 - f \right ) \tan \phi_2 \right ] \) : reduced latitude (latitude on the auxiliary sphere) \( L = L_2 - L_1 \) : difference in longitude of two points \( \lambda_1, \lambda_2 \) : longitude of the points on the auxiliary sphere \( \alpha_1, \alpha_2 \) : forward azimuths at the points \( \alpha \) : azimuth at the equator \( s \) : ellipsoidal distance between the two points \( \sigma \) : arc length between points on the auxiliary sphere

Using the values defined above, and the eight equations below, the process of iterating towards a solution can begin. Starting with the first equation on top, sequentially progress towards the last equation at the bottom where the output is \( \lambda \). Once through the first iteration, simply use this newly calculated value of \( \lambda \) and dump it back into the first equation and repeat the process. The end game here is to minimize the value of \( \lambda \).

When the difference between the current value of \( \lambda \) and the value of \( \lambda \) from the previous iteration is less than the convergence tolerance, we can move on to the final stage (a convergence tolerance of 10e-12 translates to an accuracy of approximately 0.06mm).