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 on observation points (Xp,Yp,Zp) placed on the surface z=0 of a half-space, 
% according to the equations (27)-(29) 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, depth l_b and height l_c of the source;
%
% The thermo-poro-elastic strain variation 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;
    
    for ind_i=1:2
        for ind_j=1:2
            for ind_k=1:2 
                    
                % distance 
                R=sqrt((Xobs(i)-x_extremes(ind_i)).^2+(Yobs(i)-y_extremes(ind_j)).^2+(Zobs(i)-z_extremes(ind_k)).^2); 
                mu=(-1)^ind_i*(-1)^ind_j*(-1)^ind_k; 

                %% x component
                arg_x=mu*((Xobs(i)-x_extremes(ind_i))*atan(((Zobs(i)-z_extremes(ind_k))*(Yobs(i)-y_extremes(ind_j)))/((Xobs(i)-x_extremes(ind_i))*R))-(Zobs(i)-z_extremes(ind_k))*log(R+(Yobs(i)-y_extremes(ind_j)))-(Yobs(i)-y_extremes(ind_j))*log(R+(Zobs(i)-z_extremes(ind_k))));
                sum_x=sum_x+arg_x;
                U1=(1+ni)/(3*pi)*eps_0*sum_x;  

                %% y component 
                arg_y=mu*((Yobs(i)-y_extremes(ind_j))*atan(((Zobs(i)-z_extremes(ind_k))*(Xobs(i)-x_extremes(ind_i)))/((Yobs(i)-y_extremes(ind_j))*R))-(Zobs(i)-z_extremes(ind_k))*log(R+(Xobs(i)-x_extremes(ind_i)))-(Xobs(i)-x_extremes(ind_i))*log(R+(Zobs(i)-z_extremes(ind_k))));
                sum_y=sum_y+arg_y;
                U2=(1+ni)/(3*pi)*eps_0*sum_y;               
                    
                %% z component 
                arg_z=mu*((Zobs(i)-z_extremes(ind_k))*atan(((Xobs(i)-x_extremes(ind_i))*(Yobs(i)-y_extremes(ind_j)))/((Zobs(i)-z_extremes(ind_k))*R))-(Xobs(i)-x_extremes(ind_i))*log(R+(Yobs(i)-y_extremes(ind_j)))-(Yobs(i)-y_extremes(ind_j))*log(R+(Xobs(i)-x_extremes(ind_i))));
                sum_z=sum_z+arg_z;
		        U3=(1+ni)/(3*pi)*eps_0*sum_z;                                     
                    
            end
        end
    end
		        
    %% 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