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


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%% FORWARD MODEL %%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% ANALYTICAL SOLUTION for the TRIAXIAL 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 triaxial ellipsoid on observation points (Xp,Yp,Zp) placed on the surface z=0 of a half-space, 
% according to the equations (12)-(14) 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, b and c of the semi-axes along the x, y and z 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)] 



% 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;
% 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; % Zc



%% SEMI-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);  
    
    fA = @(x)(1./((a.^2+x).*sqrt((a.^2+x).*(b.^2+x).*(c.^2+x))));     
    A_lambda = quadgk(fA,lambda,Inf,'RelTol',1e-12,'AbsTol',1e-15);
        
    fB = @(x)(1./((b.^2+x).*sqrt((a.^2+x).*(b.^2+x).*(c.^2+x))));     
    B_lambda = quadgk(fB,lambda,Inf,'RelTol',1e-12,'AbsTol',1e-15);
        
    fC = @(x)(1./((c.^2+x).*sqrt((a.^2+x).*(b.^2+x).*(c.^2+x))));     
    C_lambda = quadgk(fC,lambda,Inf,'RelTol',1e-12,'AbsTol',1e-15);
        
    U1=-2*(1+ni)/3*a*b*c*eps_0*x1*A_lambda;
    U2=-2*(1+ni)/3*a*b*c*eps_0*x2*B_lambda;
    U3=-2*(1+ni)/3*a*b*c*eps_0*x3*C_lambda;
            
        
    %% 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;