Fortran library for Geodesics  2.0
Introduction to the geodesic problem

The shortest path between two points on the ellipsoid at (lat1, lon1) and (lat2, lon2) is called the geodesic. Its length is s12 and the geodesic from point 1 to point 2 has forward azimuths azi1 and azi2 at the two end points.

Traditionally two geodesic problems are considered:

  • the direct problem – given lat1, lon1, s12, and azi1, determine lat2, lon2, and azi2. This is solved by the function geod_direct().
  • the inverse problem – given lat1, lon1, and lat2, lon2, determine s12, azi1, and azi2. This is solved by the function geod_inverse().

The ellipsoid is specified by its equatorial radius a (typically in meters) and flattening f. The routines are accurate to round off with double precision arithmetic provided that |f| < 1/50; for the WGS84 ellipsoid, the errors are less than 15 nanometers. (Reasonably accurate results are obtained for |f| < 1/5.) For a prolate ellipsoid, specify f < 0.

The routines also calculate several other quantities of interest

  • S12 is the area between the geodesic from point 1 to point 2 and the equator; i.e., it is the area, measured counter-clockwise, of the quadrilateral with corners (lat1,lon1), (0,lon1), (0,lon2), and (lat2,lon2).
  • m12, the reduced length of the geodesic is defined such that if the initial azimuth is perturbed by dazi1 (radians) then the second point is displaced by m12 dazi1 in the direction perpendicular to the geodesic. On a curved surface the reduced length obeys a symmetry relation, m12 + m21 = 0. On a flat surface, we have m12 = s12.
  • M12 and M21 are geodesic scales. If two geodesics are parallel at point 1 and separated by a small distance dt, then they are separated by a distance M12 dt at point 2. M21 is defined similarly (with the geodesics being parallel to one another at point 2). On a flat surface, we have M12 = M21 = 1.
  • a12 is the arc length on the auxiliary sphere. This is a construct for converting the problem to one in spherical trigonometry. a12 is measured in degrees. The spherical arc length from one equator crossing to the next is always 180°.

If points 1, 2, and 3 lie on a single geodesic, then the following addition rules hold:

  • s13 = s12 + s23
  • a13 = a12 + a23
  • S13 = S12 + S23
  • m13 = m12 M23 + m23 M21
  • M13 = M12 M23 − (1 − M12 M21) m23 / m12
  • M31 = M32 M21 − (1 − M23 M32) m12 / m23

The shortest distance returned by the solution of the inverse problem is (obviously) uniquely defined. However, in a few special cases there are multiple azimuths which yield the same shortest distance. Here is a catalog of those cases:

  • lat1 = −lat2 (with neither point at a pole). If azi1 = azi2, the geodesic is unique. Otherwise there are two geodesics and the second one is obtained by setting [azi1, azi2] → [azi2, azi1], [M12, M21] → [M21, M12], S12 → −S12. (This occurs when the longitude difference is near ±180° for oblate ellipsoids.)
  • lon2 = lon1 ± 180° (with neither point at a pole). If azi1 = 0° or ±180°, the geodesic is unique. Otherwise there are two geodesics and the second one is obtained by setting [azi1, azi2] → [−azi1, −azi2], S12 → −S12. (This occurs when lat2 is near −lat1 for prolate ellipsoids.)
  • Points 1 and 2 at opposite poles. There are infinitely many geodesics which can be generated by setting [azi1, azi2] → [azi1, azi2] + [d, −d], for arbitrary d. (For spheres, this prescription applies when points 1 and 2 are antipodal.)
  • s12 = 0 (coincident points). There are infinitely many geodesics which can be generated by setting [azi1, azi2] → [azi1, azi2] + [d, d], for arbitrary d.

The area of a geodesic polygon can be determined by summing −S12 for successive edges of the polygon (S12 is negated so that clockwise traversal of a polygon gives a positive area). However, if the polygon encircles a pole, the sum must be adjusted by ±A/2, where A is the area of the full ellipsoid, with the sign chosen to place the result in (−A/2, A/2].

The principal advantages of these algorithms over previous ones (e.g., Vincenty, 1975) are

  • accurate to round off for |f| < 1/50;
  • the solution of the inverse problem is always found;
  • differential and integral properties of geodesics are computed.