function [Ur,Uy,Uz] = deformation_triaxial_ellipsoid_radial(rho_radial_obs,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=[Ur,Uz] 
% in the radial and vertical components Ur and Uz induced by a triaxial ellipsoid on the surface z=0 of a half-space.
% The displacements were computed over a ground surface profile rho_radial_obs along the x-axis and with respect to the centre of the source 


Yobs=zeros(1,length(rho_radial_obs));
Zobs=Yobs;



% 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(rho_radial_obs)
    
    x1=rho_radial_obs(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_r(i)=U_rot(1);
    U_y(i)=U_rot(2);
    U_z(i)=U_rot(3); 
            
end

Ur=-U_r; 
Uy=-U_y;
Uz=U_z;