Fortran library for Geodesics
2.0
|
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... | |
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
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 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
If points 1, 2, and 3 lie on a single geodesic, then the following addition rules hold:
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:
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 and licensed under the MIT/X11 License. For more information, see es@k arney .comhttps://geographiclib.sourceforge.io/
Definition in file geodesic.for.
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.
[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
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
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:
Definition at line 182 of file geodesic.for.
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.
[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
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:
Definition at line 524 of file geodesic.for.
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.
[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.
subroutine geover | ( | integer | major, |
integer | minor, | ||
integer | patch | ||
) |
Return the version numbers for this package.
[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.
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.