function [Srr,Srt,Stt,Sxx,Sxy,Syy,taumax] = uniaxialhole2(S,X,Y)
% Stresses around a hole of unit radius in a plate under a uniaxial tension
% The solution is from Barber (1992, p. 94).
% S is the value of the uniaxial far-field tension along the x-axis
% X is the matrix (or vector) of x-coordinates of observation points
% Y is the matrix (or vector) of y-coordinates of observation points
% For a grid of points, X and Y can be found easily using
% Matlab's meshgrid command: [X,Y] = meshgrid(x,y)
% Type
% help meshgrid
% for how to use the meshgrid command
% Example: [Srr,Srt,Stt,Sxx,Syx,Syy,taumax] = uniaxialhole2(S,X,Y)
% Calculate the radial distance to the first, second, and fourth powers
% for each point on the observation grid
r = sqrt(X.*X + Y.*Y);
r2 = X.*X + Y.*Y;
r4 = r2.*r2;
% Calculate theta for each point on the observation grid
theta = atan2(Y,X);
% Calculate the trig functions of theta for each point on the grid
cost = cos(theta); % cosine of angle between x and r directions
sint = ;
cos2t = ;
sin2t = ;
cost2 = ;
sint2 = ;
axr = ; % cosine of angle between x and r directions
ayr = ; % cosine of angle between y and r directions
axt = ; % cosine of angle between x and t directions
ayt = ; % cosine of angle between y and t directions
% Calculate the stresses in a polar reference frame at the grid points
Srr = ;
Srt = ;
Str = ;
Stt = ;
% Because some grid points might be inside the circle, which is
% traction-free, check for points inside the circle and set the
% stresses at those points to zero
index = find(r<1);
Srr(index) = zeros(size(index));
Srt(index) = zeros(size(index));
Str(index) = zeros(size(index));
Stt(index) = zeros(size(index));
% Calculate the stresses in an x,y reference frame at the grid points
% See equation 5.30 in Chou and Pagano, or eqs. 1.7-1.9 of Barber,
% or equations 4.3-4.9 of course notes for GG711
Sxx = ;
Sxy = ;
Syy = ;
% Calculate the maximum shear stress
% See equation 1.15 of Chou and Pagano
taumax = ;