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