I need a function to accurately calculate the distance, initial azimuth, and final azimuth between two points defined by latitude and longitude on a reference ellipsoid.

The correct solution is given by the inverse geodetic problem of Vincenty's formulae. I have found implementations in java and C#:

https://github.com/mgavaghan/geodesya Fortran version from NOAA:

https://www.ngs.noaa.gov/TOOLS/Inv_Fwd/Inv_Fwd.htmland in C++, java, js, matlab, python, (legacy) Fortran, and (legacy) C:

https://geographiclib.sourceforge.io/The first is likely the easiest to translate, the second is spaghetti code (apologies to the coder, but looks like a lot of work to translate), the third seems the most bombproof in terms of unusual cases, but also looks like a lot of work to translate.

Is anyone aware of solid Delphi/Pascal code for this problem? I don't want to dive in if it is already done.