function [Ux,Uy,Uz] = deformation_prism(Xobs,Yobs,Zobs,Xc,Yc,depth,l_a,l_b,l_c,alpha,delta,gamma,eps_0,ni)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%% FORWARD MODEL %%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% ANALYTICAL SOLUTION for the RECTANGULAR PRISM %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



% The subroutine computes the thermo-poro-elastic ground deformation U_calc=[Ux,Uy,Uz], 
% in the components Ux, Uy and Uz, induced by a rectangular prism at observation points (Xp,Yp,Zp) placed on the surface z=0 of a half-space, 
% according to the equations (26)-(28) in the paper 



%%%%%%%%%%
%% INPUTS:
%%%%%%%%%%

% Xobs=[Xobs(1),...,Xobs(n)], Yobs=[Yobs(1),...,Yobs(n)], Zobs=[Zobs(1),...,Zobs(n)], with n the number of the observation points
% and P=(Xobs(i),Yobs(i),Zobs(i)) the coordinates of the i-th observation point in the 3D Cartesian coordinate system;
%
% The coordinates (Xc,Yc) of the center of the source;
%
% The upper limit depth of the source;
%
% Length l_a, width l_b and height l_c of the source;
%
% The stress-free strain eps_0 at the source;
%
% The angles alpha, delta and gamma of rotation about the z, y and x axes, respectively;
%
% The Poisson's ratio ni=0.25



%%%%%%%%%%%
%% OUTPUTS:
%%%%%%%%%%%

% The ground deformation at the observation points P in the components
% Ux=[Ux(1),...,Ux(n)], Uy=[Uy(1),...,Uy(n)] and Uz=[Uz(1),...,Uz(n)] 



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



%% TRASLATION OF THE OBSERVATION POINTS
% traslation of the observation points along x and y directions in such a way that (Xc,Yc)=(0,0)
Xobs=Xobs-Xc;
Yobs=Yobs-Yc;
% traslation of the observation points along z direction knowing the depth d of the top in such a way that Zc=0
[R] = rotation(alpha,delta,gamma);    
Rot=R';
P0 = [0 0 0]'; % The centroid (Xc,Yc,Zc)
P1 = P0+l_a*R(:,1)/2+l_c*R(:,3)/2-l_b/2*R(:,2);
P2 = P1-l_c*R(:,3);
P3 = P2-l_a*R(:,1);
P4 = P1-l_a*R(:,1);
Q1 = P0+l_a*R(:,1)/2+l_c*R(:,3)/2+l_b/2*R(:,2);
Q2 = Q1-l_c*R(:,3);
Q3 = Q2-l_a*R(:,1);
Q4 = Q1-l_a*R(:,1);
x=max([P1(3),P2(3),P3(3),P4(3),Q1(3),Q2(3),Q3(3),Q4(3)]);
Zobs=Zobs+x-depth;


%% RECTANGULAR PRISM EXTREMES
% along x direction with x1<x2
x1=-l_a/2; % [m]
x2=l_a/2;
% along y direction with y1<y2
y1=-l_b/2;
y2=l_b/2;
% along z direction with z1>z2
z1=l_c/2; 
z2=-z1; 
% vectors of extremes
x_extremes=[x1,x2];
y_extremes=[y1,y2];
z_extremes=[z1,z2];
    
	

	

%% ANALYTICAL FORMULATION
	
for i=1:length(Xobs)


    %% Resolution in the rotated coordinate system

    %% ROTATION OF THE CARTESIAN COORDINATE SYSTEM
    P_rot=Rot*[Xobs(i);Yobs(i);Zobs(i)]; 
    Xobs_rot=P_rot(1);
    Yobs_rot=P_rot(2);
    Zobs_rot=P_rot(3);
    % coordinates of the observation points P in the body coordinate system
    Xobs(i)=Xobs_rot;
    Yobs(i)=Yobs_rot;
    Zobs(i)=Zobs_rot;
    
    sum_z=0;
    sum_x=0;
    sum_y=0;
    
    ind_i=[1:2];
    ind_j=ind_i;
    ind_k=ind_i;

    [ind_i,ind_j,ind_k]=meshgrid(ind_i,ind_j,ind_k);
    ind_i=ind_i(:);
    ind_j=ind_j(:);
    ind_k=ind_k(:);
    
    for j=1:length(ind_i)

                    
        % distance 
        R=sqrt((Xobs(i)-x_extremes(ind_i(j))).^2+(Yobs(i)-y_extremes(ind_j(j))).^2+(Zobs(i)-z_extremes(ind_k(j))).^2); 
        mu=(-1)^ind_i(j)*(-1)^ind_j(j)*(-1)^ind_k(j); 

        %% x component
        arg_x=mu*((Xobs(i)-x_extremes(ind_i(j)))*atan(((Zobs(i)-z_extremes(ind_k(j)))*(Yobs(i)-y_extremes(ind_j(j))))/((Xobs(i)-x_extremes(ind_i(j)))*R))-(Zobs(i)-z_extremes(ind_k(j)))*log(R+(Yobs(i)-y_extremes(ind_j(j))))-(Yobs(i)-y_extremes(ind_j(j)))*log(R+(Zobs(i)-z_extremes(ind_k(j)))));
        sum_x=sum_x+arg_x;

        %% y component 
        arg_y=mu*((Yobs(i)-y_extremes(ind_j(j)))*atan(((Zobs(i)-z_extremes(ind_k(j)))*(Xobs(i)-x_extremes(ind_i(j))))/((Yobs(i)-y_extremes(ind_j(j)))*R))-(Zobs(i)-z_extremes(ind_k(j)))*log(R+(Xobs(i)-x_extremes(ind_i(j))))-(Xobs(i)-x_extremes(ind_i(j)))*log(R+(Zobs(i)-z_extremes(ind_k(j)))));
        sum_y=sum_y+arg_y;               
                    
        %% z component 
        arg_z=mu*((Zobs(i)-z_extremes(ind_k(j)))*atan(((Xobs(i)-x_extremes(ind_i(j)))*(Yobs(i)-y_extremes(ind_j(j))))/((Zobs(i)-z_extremes(ind_k(j)))*R))-(Xobs(i)-x_extremes(ind_i(j)))*log(R+(Yobs(i)-y_extremes(ind_j(j))))-(Yobs(i)-y_extremes(ind_j(j)))*log(R+(Xobs(i)-x_extremes(ind_i(j)))));
        sum_z=sum_z+arg_z;                                     
                    

    end

    U1=sum_x;
    U2=sum_y;
    U3=sum_z;
		        
    %% Deformation components in the Cartesian coordinate system 
    U_rot = Rot'*[U1;U2;U3];
    Ux(i)=U_rot(1);
    Uy(i)=U_rot(2);
    Uz(i)=U_rot(3);
            
end  

Ux=(1+ni)/(3*pi)*eps_0*Ux;
Uy=(1+ni)/(3*pi)*eps_0*Uy;
Uz=(1+ni)/(3*pi)*eps_0*Uz;