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
