      subroutine ignewt(del,stp,grav,tilt)
c  returns the integrated green function for newtonian gravity and tilt;
c the function is the integral over the interval cntered at del with
c width stp (ie, the interval [del-stp/2,del+stp/2],
c del and stp both being in radians
c  the height correction is included in the green functions,
c the station height in meters being passed as ht in the common block
c stloc
c
c  $$$$$$calls only system routines
c
      save iflg,eps,eps1,eps2,eps3,g2,em
      common/stloc/ct,st,rlam,ht
      data gn/6.67e-11/,a/6.371e6/,iflg/0/
c  on first call do part of height correction, and find constant parts
      if(iflg.eq.0) then
         iflg = 1
         eps = ht/a
         eps1=1.+eps
         eps2=eps*eps
         eps3=2*eps1*eps1
         g2 = 2*gn/(1.+1.5*eps)
	 g = 9.7803327*(1+.005279*ct*ct) - 3.08e-6*ht
	 em = gn/g
      endif
      c1 = sin(.25*stp)
      sthl = .5*stp
c sign for gravity makes + acceleration up
c use different approximations if < or > 0.05 rad (2.9 deg)
      if(del.ge.0.05) then
        grav = -g2*cos(del/2.)*c1*(1. + eps/(cos(sthl)-cos(del)))
        tilt = em*(-2.*sin(del/2.)*c1 +
     1    alog(tan((del+sthl)/4.)/tan((del-sthl)/4.)))
      endif
      if(del.lt.0.05) then
         d1 = del + stp/2
         d2 = del - stp/2
         grav = -gn*(eps1*d1*d1-2*eps)/(eps3*sqrt(eps1*d1*d1+eps2))
         grav = grav 
     1         + gn*(eps1*d2*d2-2*eps)/(eps3*sqrt(eps1*d2*d2+eps2))
         tilt = em*(d2/sqrt(eps2+d2*d2) - d1/sqrt(eps2+d1*d1) + 
     1      alog( (d1+sqrt(eps2+d1*d1))/(d2+sqrt(eps2+d2*d2)) ) )
      endif
      return
      end
