Fortran library for Geodesics  2.0
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

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]athe equatorial radius (meters).
[in]fthe flattening of the ellipsoid. Setting f = 0 gives a sphere. Negative f gives a prolate ellipsoid.
[in]lat1latitude of point 1 (degrees).
[in]lon1longitude of point 1 (degrees).
[in]azi1azimuth at point 1 (degrees).
[in]s12a12if 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]flagsa bitor'ed combination of the arcmode and unroll flags.
[out]lat2latitude of point 2 (degrees).
[out]lon2longitude of point 2 (degrees).
[out]azi2(forward) azimuth at point 2 (degrees).
[in]omaska bitor'ed combination of mask values specifying which of the following parameters should be set.
[out]a12s12if 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]m12reduced length of geodesic (meters).
[out]MM12geodesic scale of point 2 relative to point 1 (dimensionless).
[out]MM21geodesic scale of point 1 relative to point 2 (dimensionless).
[out]SS12area 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 lon2lon1 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
program geoddirect
A simple program to solve the direct geodesic problem.
Definition: geoddirect.for:9

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]athe equatorial radius (meters).
[in]fthe flattening of the ellipsoid. Setting f = 0 gives a sphere. Negative f gives a prolate ellipsoid.
[in]lat1latitude of point 1 (degrees).
[in]lon1longitude of point 1 (degrees).
[in]lat2latitude of point 2 (degrees).
[in]lon2longitude of point 2 (degrees).
[out]s12distance from point 1 to point 2 (meters).
[out]azi1azimuth at point 1 (degrees).
[out]azi2(forward) azimuth at point 2 (degrees).
[in]omaska bitor'ed combination of mask values specifying which of the following parameters should be set.
[out]a12arc length from point 1 to point 2 (degrees).
[out]m12reduced length of geodesic (meters).
[out]MM12geodesic scale of point 2 relative to point 1 (dimensionless).
[out]MM21geodesic scale of point 1 relative to point 2 (dimensionless).
[out]SS12area 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
program geodinverse
A simple program to solve the inverse geodesic problem.
Definition: geodinverse.for:9

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]athe equatorial radius (meters).
[in]fthe flattening of the ellipsoid. Setting f = 0 gives a sphere. Negative f gives a prolate ellipsoid.
[in]latsan array of the latitudes of the vertices (degrees).
[in]lonsan array of the longitudes of the vertices (degrees).
[in]nthe number of vertices.
[out]AAthe (signed) area of the polygon (meters2).
[out]PPthe 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]majorthe major version number.
[out]minorthe minor version number.
[out]patchthe 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.