function [Ux,Uy,Uz] = deformation_prolate_oblate_ellipsoid(Xobs,Yobs,Zobs,Xc,Yc,depth,a,b,alpha,delta,gamma,eps_0,ni)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%% FORWARD MODEL %%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% ANALYTICAL SOLUTION for the PROLATE ELLIPSOID %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



% The subroutine computes the thermo-poro-elastic ground deformation U_calc=[Ux,Uy,Uz], 
% in the components Ux, Uy and Uz, induced by a prolate/oblate ellipsoid on observation points (Xp,Yp,Zp) placed on the surface z=0 of a half-space, 
% according to the equations (18)-(20) 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;
%
% The length a and b of the semi-axes along the x and y direction
%
% 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)] 





c=b; 


% traslation of the observation points along x and y directions
Xobs=Xobs-Xc;
Yobs=Yobs-Yc;
% compute Zc knowing d (depth of the top)
[R_1]=rotation(alpha,delta,gamma);
alpha2=0;        
delta2=180; 
gamma2=0; 
[R_2] = rotation(alpha2,delta2,gamma2);                  
R=R_2*R_1;
Rot=R';

P0 = [0 0 0]'; % The centroid (Xc,Yc,Zc)
P1 = P0+b*R(:,2);
P2 = P0-b*R(:,2);    
Q1 = P0+c*R(:,3);
Q2 = P0-c*R(:,3);
R1 = P0+a*R(:,1);
R2 = P0-a*R(:,1);   
x=max([P1(3),P2(3),Q1(3),Q2(3),R1(3),R2(3)]);
Zc=-x+depth+Zobs;


%% ANALYTICAL FORMULATION

for i=1:length(Xobs) 

    x1=Xobs(i);
    x2=Yobs(i);

    %% Resolution in the rotated coordinate system

    %% ROTATION OF THE CARTESIAN COORDINATE SYSTEM
    % coordinates of the points (Xobs,Yobs,Zc) in the body coordinate system
    P_rot=Rot*[x1;x2;Zc(i)];
    x1=P_rot(1);
    x2=P_rot(2);
    x3=P_rot(3);


    
    p2=a^2+b^2+c^2-x1^2-x2^2-x3^2; 
    p1=a^2*b^2+b^2*c^2+c^2*a^2-(b^2+c^2)*x1^2-(c^2+a^2)*x2^2-(b^2+a^2)*x3^2; 
    p0=a^2*b^2*c^2-b^2*c^2*x1^2-c^2*a^2*x2^2-a^2*b^2*x3^2; 
    poly=[1 p2 p1 p0];
    lambdaT=roots(poly);
    lambda=max(lambdaT);  
    
    
    f1=sqrt((a^2-b^2)/(a^2+lambda))-log((sqrt((a^2-b^2))+sqrt((a^2+lambda)))/sqrt(b^2+lambda));
    f2=log((sqrt((a^2-b^2))+sqrt((a^2+lambda)))/sqrt(b^2+lambda))-sqrt(a^2-b^2)*sqrt(a^2+lambda)/(b^2+lambda);  
    
    
    U1=4*(1+ni)/3*a*b^2*eps_0*x1/((a^2-b^2)^1.5)*f1;     
    U2=2*(1+ni)/3*a*b^2*eps_0*x2/((a^2-b^2)^1.5)*f2;
    U3=2*(1+ni)/3*a*b^2*eps_0*x3/((a^2-b^2)^1.5)*f2;
    

    %% Deformation components in the Cartesian coordinate system       
    U_rot = Rot'*[U1;U2;U3];
    U_x(i)=U_rot(1);
    U_y(i)=U_rot(2);
    U_z(i)=U_rot(3);
    
end

Ux=-U_x; 
Uy=-U_y;
Uz=U_z;