Terrain correction for gravity data using Nagy, 1973
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