function [Ur,Uz] = deformation_cylinder(rho_radial_obs,Zobs,depth,R,h,eps_0,ni)



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%% FORWARD MODEL %%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% SEMI-ANALYTICAL SOLUTION for the CYLINDER %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



% 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 cylinder on the surface z=0 of a half-space
% according to the formulas (7)-(8) in the paper. 



%%%%%%%%%%
%% INPUTS:
%%%%%%%%%%

% rho_radial_obs=[rho_radial_obs(1),...,rho_radial_obs(n)], Z_obs=[Z_obs(1),...,Z_obs(n)], with n the number of the observation points
% and P=(rho_radial_obs(i),Z_obs(i)) the coordinates of the i-th observation point in the 2D cylindrical coordinate system;
%
% The coordinates (Xc,Yc) of the center of the source;
%
% The upper limit depth of the source;
%
% The radius R of the source;
%
% The height h of the source;
%
% The thermo-poro-elastic strain variation eps_0 at the source;
%
% The Poisson's ratio ni=0.25



%%%%%%%%%%%
%% OUTPUTS:
%%%%%%%%%%%

% The ground deformation at the observation points P in the components
% Ur=[Ur(1),...,Ur(n)] and Uz=[Uz(1),...,Uz(n)] 



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


lim_sup=depth; % upper limit of the source
lim_inf=depth-h; % lower limit of the source


% we compute the deformation as the sum of the deformations induced by 100 smaller cylindrical sources
m=10; % number of block inside the source along the r direction
l=10; % number of block inside the source along the z direction


amp_r=R/m; % width of each block along the r direction
dr=amp_r*ones(m,1); % vector
dr_tot=[dr';dr';dr';dr';dr';dr';dr';dr';dr';dr']; % matrix (to avoid for cicle)
rho_radial=[amp_r/2:amp_r:R-amp_r/2]; % vector of the coordinates of the center of each small cylinder along the r direction
% rho_matrix=[];
% for i=1:m
%     rho_matrix=[rho_matrix; rho_radial];
% end
rho_matrix=[rho_radial;rho_radial;rho_radial;rho_radial;rho_radial;rho_radial;rho_radial;rho_radial;rho_radial;rho_radial]; % matrix (to avoid for cicle)


dz=h/l; % width of each block along the z direction
Z=[lim_sup-dz/2:-dz:lim_inf+dz/2]'; % vector of the coordinates of the center of each small cylinder along the z direction (from top to bottom)
% Z_matrix=[];
% for i=1:l
%     Z_matrix=[Z_matrix Z];
% end
Z_matrix=[Z Z Z Z Z Z Z Z Z Z]; % matrix (to avoid for cicle)


rho_2=rho_matrix(:);
Z_2=Z_matrix(:);


for j=1:length(rho_radial_obs)
    parfor i=1:prod(size(rho_matrix))  

       fr = @(x)(  (rho_radial_obs(j)-rho_2(i)*cos(x)).*(rho_2(i).^2-2*rho_2(i)*rho_radial_obs(j)*cos(x)+(Zobs(j)-Z_2(i)).^2+rho_radial_obs(j).^2).^(-1.5) );
       fz = @(x)((rho_2(i).^2-2*rho_2(i)*rho_radial_obs(j)*cos(x)+(Zobs(j)-Z_2(i)).^2+rho_radial_obs(j).^2).^(-1.5));     
       
       integralz = quad(fz,0,2*pi); 
       integralr = quad(fr,0,2*pi);

       facz(i,j)=rho_2(i)*(Zobs(j)-Z_2(i))*integralz; 
       facr(i,j)=rho_2(i)*integralr;

    end
end



Ur=(1+ni)/(3*pi)*facr'*(eps_0.*dr_tot(:))*dz; 
Uz=(1+ni)/(3*pi)*facz'*(eps_0.*dr_tot(:))*dz; 

Ur=Ur'; % row vector;
Uz=Uz'; % row vector;