geodesic_module Module

Implementation of geodesic routines in modern Fortran

See also

• This module contains modernized versions of the Fortran code from geographiclib-fortran. The geographiclib subroutines are documented at: sourceforge These are Fortran implementation of the geodesic algorithms described in: C. F. F. Karney, Algorithms for geodesics, J. Geodesy 87, 43--55 (2013). [Apr 23, 2022 version]
• Some of the code was also split off from the Fortran Astrodynamics Toolkit (see geodesy_module.f90). [Nov 6, 2022 version]

geographiclib Notes

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

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 subroutine direct.
• the inverse problem -- given lat1, lon1, lat2, lon2, determine s12, azi1, and azi2. This is solved by the subroutine inverse.

The ellipsoid is specified by its equatorial radius a (typically in meters) and flattening f. The routines are accurate to round off with real(wp) 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:

• SS12 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 geodesic 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.
• MM12 and MM21 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 MM12 dt at point 2. MM21 is defined similarly (with the geodesics being parallel to one another at point 2). On a flat surface, we have MM12 = MM21 = 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 deg.

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

• s13 = s12 + s23
• a13 = a12 + a23
• SS13 = SS12 + SS23
• m13 = m12 MM23 + m23 MM21
• MM13 = MM12 MM23 - (1 - MM12 MM21) m23 / m12
• MM31 = MM32 MM21 - (1 - MM23 MM32) 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], [MM12, MM21] --> [MM21, MM12], SS12 --> -SS12. (This occurs when the longitude difference is near +/- 180 deg for oblate ellipsoids.)
• lon2 = lon1 +/- 180 deg (with neither point at a pole). If azi1 = 0 deg or +/- 180 deg, the geodesic is unique. Otherwise there are two geodesics and the second one is obtained by setting [azi1, azi2] --> [-azi1, -azi2], SS12 --> -SS12. (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.

These routines are a simple transcription of the corresponding C++ classes in geographiclib. Because of the limitations of Fortran 77, the classes have been replaced by simple subroutines with no attempt to save "state" across subroutine calls. Most of the internal comments have been retained. However, in the process of transcription some documentation has been lost and the documentation for the C++ classes, geodesic_module::Geodesic, geodesic_module::GeodesicLine, and geodesic_module::PolygonAreaT, should be consulted. The C++ code remains the "reference implementation". Think twice about restructuring the internals of the Fortran code since this may make porting fixes from the C++ code more difficult.

License

geographiclib-fortran Copyright (c) Charles Karney (2012-2022) charles@karney.com and licensed under the MIT/X11 License. For more information, see https://geographiclib.sourceforge.io/