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< 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
