<p>Instead of the Nagy prism formula, I suggest you to use the formula quoted in the following paper:</p>

<blockquote>
  <p>B. Banerjee, S.P. Das Gupta (1977): <br>
  "Gravitational attraction of a rectangular parallelepiped" <br>
  Geophysics, vol. 42, n. 5, pp. 1053-1055 <br>
  doi: <a href="http://dx.doi.org/10.1190/1.1440766" rel="nofollow noreferrer">10.1190/1.1440766</a></p>
</blockquote>

<p>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:</p>

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

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

<p>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".</p>

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

<p>Kind regards and best wishes to all, Francesco</p>

CC BY-SA 3.0
2017-06-27T13:24:52.710Z
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
CC BY-SA 3.0
2017-06-30T19:14:50.523Z
Thanks, I will check it out. I will read the paper and ask you a question if necessary,
jaime garbanzo
CC BY-SA 3.0
2017-06-27T12:27:48.633Z
2017-06-28T11:04:33.757Z
2017-06-28T11:04:33.757Z
francesco palmieri
<p>I have recently had to deal with the same problem and tackled it in Matlab,</p>

<p>If it helps, here's a link to it on the file exchange</p>

<p><a href="https://au.mathworks.com/matlabcentral/fileexchange/57349-nagyprism-x1-x2-y1-y2-h-rho-" rel="nofollow noreferrer">https://au.mathworks.com/matlabcentral/fileexchange/57349-nagyprism-x1-x2-y1-y2-h-rho-</a></p>

<p>Here is the code:</p>

<pre><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
</code></pre>

<p>The licence is as follows:</p>

<pre><code>Copyright (c) 2016, Jack
All rights reserved.

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.
</code></pre>

CC BY-SA 3.0
2016-05-26T12:15:12.050Z
Welcome to Earth Science! Whilst this may theoretically answer the question, [it would be preferable](//meta.stackoverflow.com/q/8259) to include the essential parts of the answer here, and provide the link for reference.
bon
CC BY-SA 3.0
2016-05-26T22:55:16.583Z
If you could post the code here as well that would be very helpful.
hichris123
CC BY-SA 3.0
2017-01-12T10:33:02.870Z
Since the Mathworks licence explicitly allows redistribution, I have taken the liberty of editing the code into the answer.
Pont
CC BY-SA 3.0
2017-02-06T15:49:53.353Z
This works perfectly, thanks you so much.
jaime garbanzo
CC BY-SA 3.0
2017-02-09T21:17:48.217Z
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
CC BY-SA 3.0
2016-05-26T11:49:15.483Z
2017-01-12T10:32:14.963Z
2017-01-12T10:32:14.963Z
Jack
<p>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.</p>

<p>So, this is the problem. I have this function
<a href="https://i.sstatic.net/gvbMP.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/gvbMP.png" alt="enter image description here"></a></p>

<ol>
<li>where x,y,z are the coordinates of a prism corner.</li>
<li>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$ </li>
<li>The equation 3 gives the $z$ component of the gravitational attraction </li>
</ol>

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

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

<p>I have the following coordinates</p>

<ul>
<li>$(x1,x2) = (-500,500)$</li>
<li>$(y1,y2) = (500,1500)$</li>
<li>$(z1,z2) = (500,1500)$</li>
</ul>

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

<p>I attach the code in python, Anyone with an answer?</p>

<pre><code>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
</code></pre>

<p>Call the function:</p>

<pre><code>x= (-500, 500)
y = (500, 1500)
z = (500, 1500)
FZ1 = nagy1973(x, y, z)
print FZ1 
</code></pre>

CC BY-SA 3.0
2016-01-28T00:13:35.320Z
Despite it being related to earth science, this question is about coding. It may be better suited to another site.
Fred
CC BY-SA 3.0
2016-01-28T00:14:58.583Z
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 (http://planetfacts.org/density-and-mass-of-earth/)
Fred
CC BY-SA 3.0
2016-01-28T02:14:33.573Z
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
CC BY-SA 3.0
2016-03-01T06:56:22.580Z
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
<geophysics><geodesy><python><gravity>
Terrain correction for gravity data using Nagy, 1973
jaime garbanzo
# Terrain correction for gravity data using Nagy, 1973

(Asked by: username_0 on 2016-01-27)

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
[![enter image description here](https://i.sstatic.net/gvbMP.png)](https://i.sstatic.net/gvbMP.png)

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 

```

**Tags:** geophysics, geodesy, python, gravity

> **Comments:**
> **username_1**: Despite it being related to earth science, this question is about coding. It may be better suited to another site.  
> **username_1**: 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 (http://planetfacts.org/density-and-mass-of-earth/)  
> **username_0**: 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.  
> **username_2**: I humbly disagree with username_1. 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!  

---

## 2 answers

### Answer by username_3 on 2017-06-27. Score: 3

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](http://dx.doi.org/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 username_3.gravimetry@libero.it

Kind regards and best wishes to all, username_3

> **Comments:**
> **username_4**: 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.  
> **username_0**: Thanks, I will check it out. I will read the paper and ask you a question if necessary,  

---

### Answer by username_5 on 2016-05-26. Score: 2

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, username_5
All rights reserved.

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.

```

> **Comments:**
> **username_6**: Welcome to Earth Science! Whilst this may theoretically answer the question, [it would be preferable](//meta.stackoverflow.com/q/8259) to include the essential parts of the answer here, and provide the link for reference.  
> **username_7**: If you could post the code here as well that would be very helpful.  
> **username_8**: Since the Mathworks licence explicitly allows redistribution, I have taken the liberty of editing the code into the answer.  
> **username_0**: This works perfectly, thanks you so much.  
> **username_0**: 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?  

---

