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