subroutine gpoly(x0,z0,xcorn,zcorn,ncorn,rho,g)
c
c  Subroutine GPOLY computes the vertical attraction of a two-
c  dimensional body with polygonal cross section.  Axes are 
c  right-handed system with y axis parallel to long direction 
c  of body and z axis vertical down.
c
c  Input parameters:
c    Observation point is (x0,z0).  Arrays xcorn and zcorn (each
c    of length ncorn) contain the coordinates of the polygon 
c    corners, arranged in clockwise order when viewed with x axis 
c    to right.  Density of body is rho.  All distance parameters 
c    in units of km; rho in units of kg/(m**3).
c
c  Output parameters:
c    Vertical attraction of gravity, g, in mGal.
c
      real km2m
      dimension xcorn(ncorn),zcorn(ncorn)
      data gamma/6.670e-11/,si2mg/1.e5/,km2m/1.e3/
      sum=0.
      do 1 n=1,ncorn
         if(n.eq.ncorn)then
            n2=1
            else
               n2=n+1
               end if
         x1=xcorn(n)-x0
         z1=zcorn(n)-z0
         x2=xcorn(n2)-x0
         z2=zcorn(n2)-z0
         r1sq=x1**2+z1**2
         r2sq=x2**2+z2**2
         if(r1sq.eq.0.)pause 'GPOLY:  Field point on corner'
         if(r2sq.eq.0.)pause 'GPOLY:  Field point on corner'         
         denom=z2-z1
         if(denom.eq.0.)denom=1.e-6
         alpha=(x2-x1)/denom
         beta=(x1*z2-x2*z1)/denom
         factor=beta/(1.+alpha**2)
         term1=0.5*(alog(r2sq)-alog(r1sq))
         term2=atan2(z2,x2)-atan2(z1,x1)
         sum=sum+factor*(term1-alpha*term2)
    1    continue
      g=2.*rho*gamma*sum*si2mg*km2m
      return
      end