% 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