- Author
- Charles F. F. Karney (charl.nosp@m.es@k.nosp@m.arney.nosp@m..com)
- Version
- 1.44
Abstract
This is a Fortran implementation of the geodesic algorithms from GeographicLib. This is a self-contained library which makes it easy to do geodesic computations for an ellipsoid of revolution in a Fortran program. It is written in Fortran 77 (avoiding features which are now deprecated) and should compile correctly with just about any Fortran compiler.
Downloading the source
The Fortran library is part of GeographicLib which available for download at
as either a compressed tar file (tar.gz) or a zip file. After unpacking the source, the Fortran library can be found in GeographicLib-1.44/legacy/Fortran. The library consists of the file geodesic.for.
Library documentation
The interface to the library is documented via doxygen in the source file. To access this, see geodesic.for.
Sample programs
Also included are 3 small test programs:
- geoddirect.for is a simple command line utility for solving the direct geodesic problem;
- geodinverse.for is a simple command line utility for solving the inverse geodesic problem;
- planimeter.for is a simple command line utility for computing the area of a geodesic polygon given its vertices.
Here, for example, is geodinverse.for
implicit none
include 'geodesic.inc'
double precision a, f, lat1, lon1, azi1, lat2, lon2, azi2, s12,
+ dummy1, dummy2, dummy3, dummy4, dummy5
integer omask
a = 6378137d0
f = 1/298.257223563d0
omask = 0
10 continue
read(*, *, end=90, err=90) lat1, lon1, lat2, lon2
if (abs(lat1) .gt. 90 .or. abs(lat2) .gt. 90) then
print 15
15 format(1x,'lat1 and lat2 must be in [-90,90]')
go to 10
end if
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
To compile, link, and run this, you would typically use
f95 -o geodinverse geodinverse.for geodesic.for
echo 30 0 29.5 179.5 | ./geodinverse
These sample programs can also be built with the supplied cmake file, CMakeLists.txt, as follows
mkdir BUILD
cd BUILD
cmake ..
make
echo 30 0 29.5 179.5 | ./geodinverse
Finally, the two programs
which are also built with cmake, provide drop-in replacements for replacements for the NGS tools FORWARD and INVERSE available from http://www.ngs.noaa.gov/PC_PROD/Inv_Fwd/. These cure two problems of the Vincenty algorithms used by NGS:
- the accuracy is "only" 0.1 mm;
- the inverse program sometimes goes into an infinite loop.
The corresponding source files
- ngsforward.for
- ngsinverse.for
- ngscommon.for
are derived from the NGS source files
- forward.for, version 2.0, dated 2002-08-21
- inverse.for, version 3.0, dated 2012-11-04
and are therefore in the public domain.
Using the library
- Optionally put in declaration section of your subroutines.
- make calls to the geodesic routines from your code. The interface to the library is documented in geodesic.for.
- Compile and link as described above.
External links
Change log
- Version 1.44 (released 2015-08-14)
- Improve accuracy of calculations by evaluating trigonometric functions more carefully and replacing the series for the reduced length with one with a smaller truncation error.
- The allowed ranges for longitudes and azimuths is now unlimited; it used to be [−540°, 540°).
- The sample programs, geoddirect and geodinverse, enforce the restriction of latitude to [−90°, 90°].
- The inverse calculation sets s12 to zero for coincident points at pole (instead of returning a tiny quantity).