# 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 1. where x,y,z are the coordinates of a prism corner.
2. 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$
3. 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

• Despite it being related to earth science, this question is about coding. It may be better suited to another site. – Fred Jan 28 '16 at 0:13
• You state the density of the Earth as 1 g/cm3. That's the density of water. The density of Earth is 5.515 g/cm3 (planetfacts.org/density-and-mass-of-earth) – Fred Jan 28 '16 at 0:14
• Yes, it is true about the density of the Earth but for this example they use 1 g/cm3. I believe it was just for exemplification. It is not really about coding. the coding part is not that hard. It is some little detail in the theory I am missing. – jaime garbanzo Jan 28 '16 at 2:14
• I humbly disagree with Fred. This is exactly the kind of question that I would like to use this website to get help with. A substantial chunk of geophysics is implementing the theory by writing code, and this is where a lot of confusion can arise. Good luck with this! – Antonio Mar 1 '16 at 6:56

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

• This 'answer' does not actually answer the question that was asked. However, it might give valuable information for people interested in this topic. Two suggestions: Please describe briefly, why Banerjee and Das Gupta (1977) actually provide a 'better' formula for this purpose. Please provide a link to the code or some information how to obtain the code. – daniel.heydebreck Jun 27 '17 at 13:24
• Thanks, I will check it out. I will read the paper and ask you a question if necessary, – jaime garbanzo Jun 30 '17 at 19:14

I have recently had to deal with the same problem and tackled it in Matlab,

If it helps, here's a link to it on the file exchange

https://au.mathworks.com/matlabcentral/fileexchange/57349-nagyprism-x1-x2-y1-y2-h-rho-

Here is the code:

% A function to perform terrain corrections using the Nagy
% prism formula.
% To calculate the terrain correction for a block of topography
% with coordinates (x1,y1) (x2,y1) (x1,y2) (x2,y2) in the x,y plane and
% height h.
%......................................
%............___________
%.........../!........./!?.............
%........../_!________/.!?.............
%..........!.!.......!..!h ............
%..........!.!.......!..!?.............
%..........!.!.......!..!?.............
%..........!.!.......!..!?.............
%..........!.!(x1,y2)!..!(x2,y2).......
%..........!./.......!./...............
%(x1,y1)...!/________!/(x2,y1).........
%
% with density rho;
% the terrain correction T is given by
% T = nagyprism(x1,x2,y1,y2,h,rho)
%
% The function can handle vectors of coordinates so that correction
% for mutiple pieces of terrain can also be calculated
function [terrc] = nagyprism(x1,x2,y1,y2,h,rho)
twopiG = 0.0419; % for g in mgal, h in m and rho in Mg/m3
G=twopiG/(2*pi);
fac11=sqrt(x1.^2+y1.^2);
fac11h=sqrt(x1.^2+y1.^2+h.^2);
fac12=sqrt(x1.^2+y2.^2);
fac12h=sqrt(x1.^2+y2.^2+h.^2);
fac21=sqrt(x2.^2+y1.^2);
fac21h=sqrt(x2.^2+y1.^2+h.^2);
fac22=sqrt(x2.^2+y2.^2);
fac22h=sqrt(x2.^2+y2.^2+h.^2);
fac2h=sqrt(y2.^2+h.^2);
fac1h=sqrt(y1.^2+h.^2);
y2h=y2.^2+h.^2;
y1h=y1.^2+h.^2;
terrc=x2.*(log( (y2+fac22)./(y2+fac22h))-...
log( (y1+fac21)./(y1+fac21h))) -...
x1.*( log( (y2+fac12)./(y2+fac12h))-log( (y1+fac11)./(y1+fac11h))) +...
y2.*(log( (x2+fac22)./(x2+fac22h))-log( (x1+fac12)./(x1+fac12h))) -...
y1.*(log( (x2+fac21)./(x2+fac21h))-log( (x1+fac11)./(x1+fac11h))) +...
h.*(asin( (y2h +y2.*fac22h)./( (y2+fac22h).*fac2h))-...
asin( (y2h +y2.*fac12h)./( (y2+fac12h).*fac2h))-...
asin( (y1h +y1.*fac21h)./( (y1+fac21h).*fac1h))+...
asin( (y1h +y1.*fac11h)./( (y1+fac11h).*fac1h)));
terrc=sum(G.*rho.*terrc);
end


The licence is as follows:

Copyright (c) 2016, Jack

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are
met:

* Redistributions of source code must retain the above copyright
notice, this list of conditions and the following disclaimer.
* Redistributions in binary form must reproduce the above copyright
notice, this list of conditions and the following disclaimer in
the documentation and/or other materials provided with the distribution

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
POSSIBILITY OF SUCH DAMAGE.

• Welcome to Earth Science! Whilst this may theoretically answer the question, it would be preferable to include the essential parts of the answer here, and provide the link for reference. – bon May 26 '16 at 12:15
• If you could post the code here as well that would be very helpful. – hichris123 May 26 '16 at 22:55
• Since the Mathworks licence explicitly allows redistribution, I have taken the liberty of editing the code into the answer. – Pont Jan 12 '17 at 10:33
• This works perfectly, thanks you so much. – jaime garbanzo Feb 6 '17 at 15:49
• I actually was trying to tested with a huge raster, and some rare things happened. try these data: x1 = 0.1 x2 =47.536 y1= 99173.917 y2 =99295.0 h= 1111.769 rho = 2.67. It gives me some negative number which couldn't be possible. Any Idea? – jaime garbanzo Feb 9 '17 at 21:17