00001 *> @file planimeter.for 00002 *! @brief A test program for area() 00003 00004 *> A simple program to compute the area of a geodesic polygon. 00005 *! 00006 *! This program reads in up to 100000 lines with lat, lon for each vertex 00007 *! of a polygon. At the end of input, the program prints the number of 00008 *! vertices, the perimeter of the polygon and its area (for the WGS84 00009 *! ellipsoid). 00010 00011 program planimeter 00012 implicit none 00013 00014 include 'geodesic.inc' 00015 00016 integer maxpts 00017 parameter (maxpts = 100000) 00018 double precision a, f, lats(maxpts), lons(maxpts), S, P 00019 integer n 00020 00021 * WGS84 values 00022 a = 6378137d0 00023 f = 1/298.257223563d0 00024 00025 n = 0 00026 10 continue 00027 if (n .ge. maxpts) go to 20 00028 read(*, *, end=20, err=20) lats(n+1), lons(n+1) 00029 n = n+1 00030 go to 10 00031 20 continue 00032 call area(a, f, lats, lons, n, S, P) 00033 print 30, n, P, S 00034 30 format(i6, 1x, f20.8, 1x, f20.3) 00035 stop 00036 end