nagyprism(x1,x2,y1,​y2,h,rho)

A function file to compute Nagy prism terrain corrections

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 all in metres.
.....................................................
..................__________...............
................/!.................../!↑.............
............../_!________./..!↑.............
..............!..!................!...!h ............
..............!..!................!...!↓.............
..............!..!................!...!↓.............
..............!..!................!...!↓.............
..............!..!(x1,y2).....!_.!(x2,y2)....
..............!./.................!./................
(x1,y1)...!/_________!/(x2,y1).....

with density rho in Mg/m3, the terrain correction T is given by
T = nagyprism(x1,x2,y1,y2,h,rho)

The function can also handle vectors of coordinates, x1,x2,y1,y2,h so that the total correction for multiple pieces of terrain can also be calculated at once.