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.
- where x,y,z are the coordinates of a prism corner.
- the equation 2 is the numerical method, where $S(xi) = -1 for x< 0, 1 for x > 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] < 0: sgnX = -1 elif x[i-1] == 0: sgnX = 0 else: sgnX = 1 if y[j-1] < 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