function [U1,U2,U3] = calc_penny(xn,yn,depth,h,P);

%depth		depth of source, m
%h 		dimensionless crack depth (depth/radius)
%xn,yn		calculation points, m
%P		effective elastic strain (excess pressure/rigidity)

%mu=matrl(2);
%nu=matrl(3);
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);
U1=zeros(ix,iy);
U2=zeros(ix,iy);
U3=zeros(ix,iy);

% 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

 [U3,Ur]=intgr(r./depth*h,fi,psi,h,Wt,t);  %calculate vertical and radial displ.
 theta=atan2(yn,xn);
 U1=Ur.*cos(theta)*P*depth/h;  % decompose radial displ.into x and y components
 U2=Ur.*sin(theta)*P*depth/h;
 U3=-U3*P*depth/h;     % in my solutions, vertical axis is positive downward; 
                % positive P should result in positive Uz and Ur.
