function [Ux,Uy,Uz] = deformation_sphere(Xobs,Yobs,Zobs,Xc,Yc,Zc,DeltaV,ni)



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%% FORWARD MODEL %%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% ANALYTICAL SOLUTION for the SPHERE %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



% The subroutine computes the thermo-poro-elastic ground deformation U_calc=[Ux,Uy,Uz], 
% in the components Ux, Uy and Uz, induced by a sphere on observation points (Xp,Yp,Zp) placed on the surface z=0 of
% a half-space (eq. (6) 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,Zc) of the center of the source;
%
% The volume variation DeltaV volume variation undergoes by the source;
%
% 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)] 



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



%% ANALYTICAL FORMULATION

R3=((Xobs-Xc).^2 + (Yobs-Yc).^2 + (Zobs-Zc).^2).^(1.5);  
Ux=(1+ni)/(3*pi)*DeltaV*(Xobs-Xc).*1./R3;
Uy=(1+ni)/(3*pi)*DeltaV*(Yobs-Yc).*1./R3;
Uz=(1+ni)/(3*pi)*DeltaV*(Zobs-Zc).*1./R3;