00001 *> @file geodinverse.for
00002 *! @brief A test program for invers()
00003
00004 *> A simple program to solve the inverse geodesic problem.
00005 *!
00006 *! This program reads in lines with lat1, lon1, lon2, lat2 and prints
00007 *! out lines with azi1, azi2, s12 (for the WGS84 ellipsoid).
00008
00009 program geodinverse
00010 implicit none
00011
00012 include 'geodesic.inc'
00013
00014 double precision a, f, lat1, lon1, azi1, lat2, lon2, azi2, s12,
00015 + dummy
00016 integer omask
00017
00018 * WGS84 values
00019 a = 6378137d0
00020 f = 1/298.257223563d0
00021
00022 omask = 0
00023
00024 10 continue
00025 read(*, *, end=90, err=90) lat1, lon1, lat2, lon2
00026 call invers(a, f, lat1, lon1, lat2, lon2,
00027 + s12, azi1, azi2, omask, dummy, dummy, dummy, dummy, dummy)
00028 print 20, azi1, azi2, s12
00029 20 format(f20.15, 1x, f20.15, 1x, f19.10)
00030 go to 10
00031 90 continue
00032 stop
00033 end
