geodesic.for File Reference

Implementation of geodesic routines in Fortran. More...

Go to the source code of this file.

Functions/Subroutines
subroutine	direct (a, f, lat1, lon1, azi1, s12a12, flags, lat2, lon2, azi2, omask, a12s12, m12, MM12, MM21, SS12)
	Solve the direct geodesic problem. More...
subroutine	invers (a, f, lat1, lon1, lat2, lon2, s12, azi1, azi2, omask, a12, m12, MM12, MM21, SS12)
	Solve the inverse geodesic problem. More...
subroutine	area (a, f, lats, lons, n, AA, PP)
	Determine the area of a geodesic polygon. More...
subroutine	geover (major, minor, patch)
	Return the version numbers for this package. More...
subroutine	geoini
	Initialize some constants. More...

Detailed Description

Implementation of geodesic routines in Fortran.

This is a Fortran implementation of the geodesic algorithms described in

• C. F. F. Karney, Algorithms for geodesics, J. Geodesy 87, 43–55 (2013); DOI: 10.1007/s00190-012-0578-z; addenda.

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 invers().

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

• 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°.

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° 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], 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, GeographicLib::Geodesic, GeographicLib::GeodesicLine, and GeographicLib::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.

Copyright (c) Charles Karney (2012-2022) charl.nosp@m.es@k.nosp@m.arney.nosp@m..com and licensed under the MIT/X11 License. For more information, see https://geographiclib.sourceforge.io/

Definition in file geodesic.for.

Function/Subroutine Documentation

direct()

subroutine direct	(	double precision	a,
		double precision	f,
		double precision	lat1,
		double precision	lon1,
		double precision	azi1,
		double precision	s12a12,
		integer	flags,
		double precision	lat2,
		double precision	lon2,
		double precision	azi2,
		integer	omask,
		double precision	a12s12,
		double precision	m12,
		double precision	MM12,
		double precision	MM21,
		double precision	SS12
	)

Solve the direct geodesic problem.

Parameters

[in]	a	the equatorial radius (meters).
[in]	f	the flattening of the ellipsoid. Setting f = 0 gives a sphere. Negative f gives a prolate ellipsoid.
[in]	lat1	latitude of point 1 (degrees).
[in]	lon1	longitude of point 1 (degrees).
[in]	azi1	azimuth at point 1 (degrees).
[in]	s12a12	if arcmode is not set, this is the distance from point 1 to point 2 (meters); otherwise it is the arc length from point 1 to point 2 (degrees); it can be negative.
[in]	flags	a bitor'ed combination of the arcmode and unroll flags.
[out]	lat2	latitude of point 2 (degrees).
[out]	lon2	longitude of point 2 (degrees).
[out]	azi2	(forward) azimuth at point 2 (degrees).
[in]	omask	a bitor'ed combination of mask values specifying which of the following parameters should be set.
[out]	a12s12	if arcmode is not set, this is the arc length from point 1 to point 2 (degrees); otherwise it is the distance from point 1 to point 2 (meters).
[out]	m12	reduced length of geodesic (meters).
[out]	MM12	geodesic scale of point 2 relative to point 1 (dimensionless).
[out]	MM21	geodesic scale of point 1 relative to point 2 (dimensionless).
[out]	SS12	area under the geodesic (meters2).

flags is an integer in [0, 4) whose binary bits are interpreted as follows

• 1 the arcmode flag
• 2 the unroll flag

If arcmode is not set, s12a12 is s12 and a12s12 is a12; otherwise, s12a12 is a12 and a12s12 is s12. It unroll is not set, the value lon2 returned is in the range [−180°, 180°]; if unroll is set, the longitude variable is "unrolled" so that lon2 − lon1 indicates how many times and in what sense the geodesic encircles the ellipsoid.

omask is an integer in [0, 16) whose binary bits are interpreted as follows

• 1 return a12
• 2 return m12
• 4 return MM12 and MM21
• 8 return SS12

lat1 should be in the range [−90°, 90°]. The value azi2 returned is in the range [−180°, 180°].

If either point is at a pole, the azimuth is defined by keeping the longitude fixed, writing lat = lat = ±(90° − ε), and taking the limit ε → 0+. An arc length greater that 180° signifies a geodesic which is not a shortest path. (For a prolate ellipsoid, an additional condition is necessary for a shortest path: the longitudinal extent must not exceed of 180°.)

Example of use:

*> @file geoddirect.for
*! @brief A test program for direct()
 
*> A simple program to solve the direct geodesic problem.
*!
*! This program reads in lines with lat1, lon1, azi1, s12 and prints out
*! lines with lat2, lon2, azi2 (for the WGS84 ellipsoid).
 
      program geoddirect
      implicit none
 
      include 'geodesic.inc'
 
      double precision a, f, lat1, lon1, azi1, lat2, lon2, azi2, s12,
     +    dummy1, dummy2, dummy3, dummy4, dummy5
      integer flags, omask
 
* WGS84 values
      a = 6378137d0
      f = 1/298.257223563d0
 
      flags = 0
      omask = 0
 
 10   continue
      read(*, *, end=90, err=90) lat1, lon1, azi1, s12
      call direct(a, f, lat1, lon1, azi1, s12, flags,
     +    lat2, lon2, azi2, omask,
     +    dummy1, dummy2, dummy3, dummy4, dummy5)
      print 20, lat2, lon2, azi2
 20   format(1x, f20.15, 1x, f20.15, 1x, f20.15)
      go to 10
 90   continue
 
      stop
      end

Definition at line 182 of file geodesic.for.

invers()

subroutine invers	(	double precision	a,
		double precision	f,
		double precision	lat1,
		double precision	lon1,
		double precision	lat2,
		double precision	lon2,
		double precision	s12,
		double precision	azi1,
		double precision	azi2,
		integer	omask,
		double precision	a12,
		double precision	m12,
		double precision	MM12,
		double precision	MM21,
		double precision	SS12
	)

Solve the inverse geodesic problem.

Parameters

[in]	a	the equatorial radius (meters).
[in]	f	the flattening of the ellipsoid. Setting f = 0 gives a sphere. Negative f gives a prolate ellipsoid.
[in]	lat1	latitude of point 1 (degrees).
[in]	lon1	longitude of point 1 (degrees).
[in]	lat2	latitude of point 2 (degrees).
[in]	lon2	longitude of point 2 (degrees).
[out]	s12	distance from point 1 to point 2 (meters).
[out]	azi1	azimuth at point 1 (degrees).
[out]	azi2	(forward) azimuth at point 2 (degrees).
[in]	omask	a bitor'ed combination of mask values specifying which of the following parameters should be set.
[out]	a12	arc length from point 1 to point 2 (degrees).
[out]	m12	reduced length of geodesic (meters).
[out]	MM12	geodesic scale of point 2 relative to point 1 (dimensionless).
[out]	MM21	geodesic scale of point 1 relative to point 2 (dimensionless).
[out]	SS12	area under the geodesic (meters2).

omask is an integer in [0, 16) whose binary bits are interpreted as follows

• 1 return a12
• 2 return m12
• 4 return MM12 and MM21
• 8 return SS12

lat1 and lat2 should be in the range [−90°, 90°]. The values of azi1 and azi2 returned are in the range [−180°, 180°].

If either point is at a pole, the azimuth is defined by keeping the longitude fixed, writing lat = ±(90° − ε), and taking the limit ε → 0+.

The solution to the inverse problem is found using Newton's method. If this fails to converge (this is very unlikely in geodetic applications but does occur for very eccentric ellipsoids), then the bisection method is used to refine the solution.

Example of use:

*> @file geodinverse.for
*! @brief A test program for invers()
 
*> A simple program to solve the inverse geodesic problem.
*!
*! This program reads in lines with lat1, lon1, lon2, lat2 and prints
*! out lines with azi1, azi2, s12 (for the WGS84 ellipsoid).
 
      program geodinverse
      implicit none
 
      include 'geodesic.inc'
 
      double precision a, f, lat1, lon1, azi1, lat2, lon2, azi2, s12,
     +    dummy1, dummy2, dummy3, dummy4, dummy5
      integer omask
 
* WGS84 values
      a = 6378137d0
      f = 1/298.257223563d0
 
      omask = 0
 
 10   continue
      read(*, *, end=90, err=90) lat1, lon1, lat2, lon2
      call invers(a, f, lat1, lon1, lat2, lon2,
     +    s12, azi1, azi2, omask,
     +    dummy1, dummy2, dummy3, dummy4, dummy5)
      print 20, azi1, azi2, s12
 20   format(1x, f20.15, 1x, f20.15, 1x, f19.10)
      go to 10
 90   continue
 
      stop
      end

Definition at line 524 of file geodesic.for.

area()

subroutine area	(	double precision	a,
		double precision	f,
		double precision, dimension(n)	lats,
		double precision, dimension(n)	lons,
		integer	n,
		double precision	AA,
		double precision	PP
	)

Determine the area of a geodesic polygon.

Parameters

[in]	a	the equatorial radius (meters).
[in]	f	the flattening of the ellipsoid. Setting f = 0 gives a sphere. Negative f gives a prolate ellipsoid.
[in]	lats	an array of the latitudes of the vertices (degrees).
[in]	lons	an array of the longitudes of the vertices (degrees).
[in]	n	the number of vertices.
[out]	AA	the (signed) area of the polygon (meters2).
[out]	PP	the perimeter of the polygon.

lats should be in the range [−90°, 90°].

Arbitrarily complex polygons are allowed. In the case of self-intersecting polygons the area is accumulated "algebraically", e.g., the areas of the 2 loops in a figure-8 polygon will partially cancel. There's no need to "close" the polygon by repeating the first vertex. The area returned is signed with counter-clockwise traversal being treated as positive.

Definition at line 962 of file geodesic.for.

geover()

subroutine geover	(	integer	major,
		integer	minor,
		integer	patch
	)

Return the version numbers for this package.

Parameters

[out]	major	the major version number.
[out]	minor	the minor version number.
[out]	patch	the patch number.

This subroutine was added with version 1.44.

Definition at line 1029 of file geodesic.for.

geoini()

subroutine geoini

Initialize some constants.

This routine is called by the main geodesic routines to initialize some constants, so usually there's no need for it to be called explicitly. However, if accessing some of the other routines, it may need to be called first.

Definition at line 1047 of file geodesic.for.