Fortran library for Geodesics
1.37
|
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, arcmod, 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... | |
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::PolygonArea, 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-2013) charl and licensed under the MIT/X11 License. For more information, see es@k arney .comhttp://geographiclib.sourceforge.net/
This library was distributed with GeographicLib 1.31.
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, | ||
logical | arcmod, | ||
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 arcmod is false, this is the distance between point 1 and point 2 (meters); otherwise it is the arc length between point 1 and point 2 (degrees); it can be negative. |
[in] | arcmod | logical flag determining the meaning of the s12a12. |
[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 arcmod is false, this is the arc length between point 1 and point 2 (degrees); otherwise it is the distance between point 1 and 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). |
omask is an integer in [0, 16) whose binary bits are interpreted as follows
lat1 should be in the range [−90°, 90°]; lon1 and azi1 should be in the range [−540°, 540°). The values of lon2 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 = 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°.)
Definition at line 170 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 between point 1 and 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 of between point 1 and 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°]; lon1 and lon2 should be in the range [−540°, 540°). 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.
Definition at line 507 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°]; lons should be in the range [−540°, 540°).
Only simple polygons (which are not self-intersecting) are allowed. 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 924 of file geodesic.for.