

%mu=matrl(2);
%nu=matrl(3);
nu2=2*(1-0.25);
xn=[0.001:100];

nis=2;        % number of sub-intervals on [0,1] on which integration is done
             % using a 16-point Gauss quadrature (i.e., total of nis*16 points)
eps=1e-6;    % solution accuracy for Fredholm integral equations (stop 
             % iterating if relative change is less than eps)

[ix,iy]=size(xn);
U3=zeros(ix,iy);
U1=zeros(ix,iy);

depth=19;	%	depth of source, m
h=0.5; 	%	dimensionless crack depth (depth/radius)
P=0.5;	%	effective elastic strain (excess pressure/rigidity)
% explicitly assign source parameters
% Solve a coupled system of Fredholm eqs. of 2nd kind for basis functions fi,psi
% t and Wt are nodes and weights of the num. integration quadrature
 [fi,psi,t,Wt]=fredholm(h,nis,eps);
 
% r=(xn.^2+yn.^2).^(1/2);      % radial distance from x0,y0 to observ. points
 r=xn;    

 [U3,Ur]=intgr(r./depth*h,fi,psi,h,Wt,t);  %calculate vertical and radial displ.
 U3=-U3*P*depth/h*nu2;     % in my solutions, vertical axis is positive downward; 
                % positive P should result in positive Uz and Ur.
 U1=U1+U3;


depth=35;	%	depth of source, m
h=1; 		%dimensionless crack depth (depth/radius)
P=0;		%effective elastic strain (excess pressure/rigidity)
%P=-1.5;		%effective elastic strain (excess pressure/rigidity)
% explicitly assign source parameters
% Solve a coupled system of Fredholm eqs. of 2nd kind for basis functions fi,psi
% t and Wt are nodes and weights of the num. integration quadrature
 [fi,psi,t,Wt]=fredholm(h,nis,eps);
 
% r=(xn.^2+yn.^2).^(1/2);      % radial distance from x0,y0 to observ. points
 r=xn;    

 [U3,Ur]=intgr(r./depth*h,fi,psi,h,Wt,t);  %calculate vertical and radial displ.
 U3=-U3*P*depth/h*nu2;     % in my solutions, vertical axis is positive downward; 
 U1=U1+U3;
plot(xn,U1);
