Terrain correction for gravity data using Nagy, 1973 I am working with an algorithm for terrain correction of gravity data. I am using Nagy, 1973. However, I could not get the results I need.

So, this is the problem. I have this function



where x,y,z are the coordinates of a prism corner.
the equation 2 is the numerical method, where $S(xi) =  -1 for x&lt; 0,  1 for x &gt; 0$ and $0$ for $x=0$. same applies $y$
but $F(x,y,z) = 0$ when $z = 0$ 
The equation 3 gives the $z$ component of the gravitational attraction 


It looks pretty simple to implement but I can't get a consitent result. 

then, I know that I have to multiply by $G$ the gravitational constant and the density of the earth. 

I have the following coordinates


$(x1,x2) = (-500,500)$
$(y1,y2) = (500,1500)$
$(z1,z2) = (500,1500)$


The earth density is $1000\ kg/m^3$. I should get a result of $2.360\ miliGals$. But I didn't. 

I attach the code in python, Anyone with an answer?

def nagy1973(x, y, z):
from math import log, sqrt, atan

# computes the total gravitational atraction of a rightrectangular Prism
# input x =  tuple with x1 and x2 coordinates
#       y =  tuple with y1 and y2 coordinates
#       z =  tuple with z1 and z2 coordinates
# output GravCorrect = a terrain correction for a right rectangular prims.
def Fz(a,b,c):
    # this is the numerical method for the gravitational
    # attraction over the Z component
    # input a = x
    #       b = y
    #       b = z
    # output F = GravCorrect    terrain correction
    a = abs(a)
    b = abs(b)
    c = abs(c)
    r = sqrt(a**2 + b**2 + c**2)
    r0 = sqrt(a**2 + b**2)

    GravCorrect = a * log((b+r0)/(b+r)) + b * log((a + r0)/(a + r)) + c* atan(a*b/(c*r))
    return GravCorrect

GravCorrAcum = 0 # this gives total z component of a prism
for k in range(1, 3):

    for j in range(1, len(y)+1):

        for i in range(1, len(x)+1):

            if x[i-1] &lt; 0:
                sgnX = -1
            elif x[i-1] == 0:
                sgnX = 0
            else:
                sgnX = 1

            if y[j-1] &lt; 0:
                sgnY = -1
            elif y[j-1] == 0:
                sgnY = 0
            else:
                sgnY = 1

            sgnIJK = (-1)**(i+j+k)

            RealSign = sgnIJK * sgnX * sgnY

            #print(sgnIJK, sgnX, sgnY, RealSign)

            GravCorr = Fz(x[i-1], y[j-1], z[k-1])

            FinalFz =  RealSign * GravCorr

            GravCorrAcum += FinalFz*6.67408e-11*1000*100000 #multiply by the gravitational constant 
                                                            #density 1000 kg/m3,100000 conversion constant to miligals 
                                                            # output miligals
           # print  FinalFz, GravCorrAcum

return GravCorrAcum


Call the function:

x= (-500, 500)
y = (500, 1500)
z = (500, 1500)
FZ1 = nagy1973(x, y, z)
print FZ1
Instead of the Nagy prism formula, I suggest you to use the formula quoted in the following paper:


  B. Banerjee, S.P. Das Gupta (1977): 
  "Gravitational attraction of a rectangular parallelepiped" 
  Geophysics, vol. 42, n. 5, pp. 1053-1055 
  doi: 10.1190/1.1440766


I wrote, in Fortran, a TC program based on that formula. If you read the paper you better understand why the algorithm proposed by Banerjee-das Gupta is easier to program, but in any case, I would like to point out some features that are described at the pp. 1055:

1)"...our expression contains only a tan inverse term facilitating computer work in all respects".

2)"The argument of the tan inverse term of our expression is simpler than the argument of sin inverse term in Nagy’s formula"...

3)"...in the Nagy's formula when the signs of x1 and x2, for example, are opposite (i.e., the y-axis is crossed), the prism is divided in two, and the integral is separately evaluated from 0 to x1, and 0 to x2, and finally they are added. In this process, some logical control is required in programming his formula. But from our expression (3), it easily follows that separation of the r.h.s. term into two is not required; similarly, when the x-axis is crossed or both x- and y-axes are crossed. This is because the signs of x’s and y’s are automatically taken cam of in our subroutine".

The formula to be programmed is quoted at pp. 1055 n° 3), but in any case I could provide you the code (in Fortran) that I wrote 30 years ago writing to palmieri.gravimetry@libero.it

Kind regards and best wishes to all, Francesco