% fig_09_27
% calculate stresses around an elliptcal hole
% under biaxial stress with intermal pressure
% Compression is positive.
% Pollard, 1973, Mathematical Geology, Vol.5, No.1, 10, 973. 
% see also Mushkelishvili 1953, and Stevenson, 1945
% M-file by JOK 21.3.2005 modified by DDP 26.07.2005

clear all, clf reset; % clear memory and figures
a = 1; % semi-major axis length 
b = .5; % semi-minor axis length, b < a
beta = 0; % angle (degrees) of s3 to x-axis
s3 = 1; % remote stress acting at beta to x-axis
s1 = 0; % remote stress acting at beta + pi/2 to x-axis
S = [s3,s1]; % vector of remote principal stresses
p = 0; % pressure in the hole
pr = 0.25; % Poisson's ratio

x=-3*a+eps:.05*a:3*a;   % define x-coords. of grid
y=-3*a+eps:.05*a:3*a;   % define y-coords. of grid
[X,Y] = meshgrid(x,y); % define Cartesian grid

SXX = p*ones(size(X)); SYY = SXX; % hydrostatic stress p
SXY = 0*SXX;

f = sqrt(a^2-b^2); % focal length f
RT = sqrt((X+f).^2+Y.^2)+sqrt((X-f).^2+Y.^2);

XI = acosh(RT/(2*f)); % coordinates xi of grid points
ET = atan2(Y,(X.*tanh(XI))); % coordinates eta of grid points

xio = atanh(b/a); % xi coordinate of hole boundary
ET(find(XI<xio)) = nan; % set coords within hole to nan
XI(find(XI<xio)) = nan; % set coords within hole to nan

% Shorthhands
s2xio = sinh(2*xio); c2xio = cosh(2*xio);
e2xio = s2xio + c2xio;

S2XI = sinh(2*XI); S2ET = sin(2*ET);
DEN = cosh(2*XI)-cos(2*ET); DEN3 = DEN.^3;

I = sinh(XI).*sin(ET);
J = cosh(XI).*cos(ET);
K = sinh(XI).*cos(ET);
L = cosh(XI).*sin(ET);

A = J.*K.*K.*K-3*J.*K.*L.*L-3*I.*K.*K.*L+I.*L.*L.*L;
B = I.*K.*K.*K-3*I.*K.*L.*L+3*J.*K.*K.*L-J.*L.*L.*L; MB = -B;
C = J.*K.*K.*K-3*J.*K.*L.*L+3*I.*K.*K.*L-I.*L.*L.*L;
D = -I.*K.*K.*K+3*I.*K.*L.*L+3*J.*K.*K.*L-J.*L.*L.*L;
E = I.*K.*K.*K-3*J.*K.*K.*L-3*I.*K.*L.*L+J.*L.*L.*L;

for k = 1:2
    be = beta*pi/180 + (k-1)*pi/2;
    s2b = sin(2*be); c2b = cos(2*be);
    smp = S(k)-p;
      
% Derivatives of the stress functions
REPSIP = .25*smp*(e2xio*c2b + (S2XI-e2xio*(S2XI*c2b+S2ET*s2b))./DEN);

REZPSIPP = -2*smp*(((1 - e2xio*c2b))*A - e2xio*s2b*B)./DEN3;         

IMZPSIPP = -2*smp*(((1 - e2xio*c2b))*MB - e2xio*s2b*A)./DEN3;

RECHIPP = -2*smp*(.25*e2xio*(c2xio*c2b - (s2xio*c2b*S2XI + c2xio*s2b*S2ET)./DEN)...
                     +((c2b-c2xio+e2xio*s2xio*c2b)*C + e2xio*c2xio*s2b*D)./DEN3);

IMCHIPP = -2*smp*(.25*e2xio*(s2xio*s2b - (c2xio*s2b*S2XI + s2xio*c2b*S2ET)./DEN)...
                     +((c2b-c2xio+e2xio*s2xio*c2b)*E + e2xio*c2xio*s2b*C)./DEN3);

% Stresses
SXX = SXX + 2*REPSIP - REZPSIPP - RECHIPP;
SYY = SYY + 2*REPSIP + REZPSIPP + RECHIPP;
SXY = SXY + IMZPSIPP + IMCHIPP;
end;

%SZZ = pr*(SXX+SYY); % out of plane normal stress
SM = (SXX+SYY)/2; % mean normal stress in xy-plane
SS = sqrt(((SXX - SYY)./2).^2 + SXY.^2); % maximum shear stress

% Plotting stresses
contourf(X,Y,(SXX),15), axis equal tight, colorbar
title('\sigma_{xx}'), ylabel('y/|a|'), xlabel('x/|a|')
figure, contourf(X,Y,(SYY),15), axis equal tight, colorbar
title('\sigma_{yy}'), ylabel('y/|a|'), xlabel('x/|a|')
figure, contourf(X,Y,(SXY),15), axis equal tight, colorbar
title('\sigma_{xy}'), ylabel('y/|a|'), xlabel('x/|a|')
%figure, contourf(X,Y,(SZZ),15), axis equal tight, colorbar
%title('\sigma_{zz}'), ylabel('y/|a|'), xlabel('x/|a|')
figure, contourf(X,Y,(SM),15), axis equal tight, colorbar
title('\sigma_{m}'), ylabel('y/|a|'), xlabel('x/|a|')
figure, contourf(X,Y,(SS),15), axis equal tight, colorbar
title('\sigma_{s}'), ylabel('y/|a|'), xlabel('x/|a|')