function dg = mLocal_prismEffect(Input,Xi,Yi,Hi,roi,Height,grid,sum_switch)
%MLOCAL_PRISMEFFECT Function for the compuation of gravitational effect of prism
% Function serves for the computation of gravitational effect of simple 
% prism in cartesian coordinate system, i.e. in x,y+height and constant 
% density. The prism height is automaticaly corrected for Earth's 
% curvature using 
% radius = 6371008 m.
% dg = mLocal_prismEffect(Input,Xi,Yi,Hi,roi,Height,grid,1)
% Input:
%   bodyvypoctu     ...     [x,y,height] of computation point
%                           vector
%   Xi              ...     x coordinate of the centre of prism upper side
%                           scalar or matrix
%   Yi              ...     y coordinate of the centre of prism upper side
%                           scalar or matrix
%   Hi              ...     height of the centre of prism upper side
%                           scalar or matrix
%   roi             ...     prism difference density
%                           scalar or matrix
%   Height           ...     prism height
%                           scalar or matrix
%   grid            ...     dimensions of prism in x/y direction [dx,dy]
%                           scalar or vector
%   sum_switch      ...     ==1 => result is the sum of all grid cells
%                           ==0 => result is a matrix
% 
% Output:
%   dg              ...     gravitational effect of prism(s), i.e. 
%                           scalar/matrix with all grid cells in m/s^2
% 
% Check: 24.10.2012, to Mod3D - OK (for fine grid/resolution),
%                    to Karcolom - OK (1 prism 2x2x1 s const.density).



%% COMPUTATION
R = 6371008;
d = sqrt((Input(1)-Xi).^2 + (Input(2)-Yi).^2);
Hi = Hi - ((sqrt(d.^2 + R^2) - R).*R)./(sqrt(d.^2 + R^2));                  % correct height for Earths curvature
clear d R;
Hi = -Hi+Input(:,3);
Xi = Xi-Input(:,1);                                                         % Create new grid = zero in point of compuation
Yi = Yi-Input(:,2);

if size(grid) == 1;
    grid(2) = grid(1);                                                      % Create y resolution if not given
end
z1 = Hi;
z2 = Hi+Height;
clear Hi Height
x1 = Yi-grid(2)/2;                                                          % Yi => xi on purpose, see rows 43-44
x2 = Yi+grid(2)/2;
clear Yi
y1 = Xi-grid(1)/2;
y2 = Xi+grid(1)/2;
clear Xi
% Create temp variables
x22py22 = x2.^2+y2.^2;
x12py22 = x1.^2+y2.^2;
x22py12 = x2.^2+y1.^2;
x12py12 = x1.^2+y1.^2;
z12 = z1.^2;
z22 = z2.^2;

G = 6.674215*10^-11;                                                        % Gravity constant

dg = -G*roi.*(x2.*(log((y2+sqrt(x22py22+z22))./(y2+sqrt(x22py22+z12)))-log((y1+sqrt(x22py12+z22))./(y1+sqrt(x22py12+z12))))-...
             x1.*(log((y2+sqrt(x12py22+z22))./(y2+sqrt(x12py22+z12)))-log((y1+sqrt(x12py12+z22))./(y1+sqrt(x12py12+z12))))+...
             y2.*(log((x2+sqrt(x22py22+z22))./(x2+sqrt(x22py22+z12)))-log((x1+sqrt(x12py22+z22))./(x1+sqrt(x12py22+z12))))-...
             y1.*(log((x2+sqrt(x22py12+z22))./(x2+sqrt(x22py12+z12)))-log((x1+sqrt(x12py12+z22))./(x1+sqrt(x12py12+z12))))+...
             z2.*(atan2((z2.*sqrt(x22py22+z22)),(x2.*y2))-atan2((z2.*sqrt(x12py22+z22)),(x1.*y2))-atan2((z2.*sqrt(x22py12+z22)),(x2.*y1))+atan2((z2.*sqrt(x12py12+z22)),(x1.*y1)))-...
             z1.*(atan2((z1.*sqrt(x12py12+z12)),(x1.*y1))-atan2((z1.*sqrt(x12py22+z12)),(x1.*y2))-atan2((z1.*sqrt(x22py12+z12)),(x2.*y1))+atan2((z1.*sqrt(x22py22+z12)),(x2.*y2))));
if sum_switch == 1
    dg = sum(sum(dg));                                                 % final effect = sum of all grid cells in uGal (=10^-8 m/s^2 = 10 nm/s^2)
end