function [ ] = PlotMohrsCircle (Sx, Sy, Sz, Txy, Tyz, Tzx)
%
%  This function uses the applied stresses to plot a Mohr's circle diagram.
%  You must manually set the axises to "Auto" in the graph to see the
%  plots.  Options->Axis Limits (Auto, Auto)
%
%  Input
%    Sx, Sy, Sz        the applied normal stresses
%    Txy, Tyz, Tzx     the applied shear stresses
%

%  Compute the principal stresses

[ S1, S2, S3, T12, T23, T13 ] = PrnStresses (Sx, Sy, Sz, Txy, Tyz, Tzx);

%  draw circle 13

x = (S1 + S3) / 2;
r = (S1 - S3) / 2;
if (r > 0)
    Circle(x, 0, r);
end
hold on

%  Draw circle 12

x = (S1 + S2) / 2;
r = (S1 - S2) / 2;
if (r > 0)
    Circle(x, 0, r);
end
hold on

%  Draw Circle 23

x = (S2 + S3) / 2;
r = (S2 - S3) / 2;
if (r > 0)
    Circle(x, 0, r);
end

%  Draw in the Axes

Length = (S1 - S3) * .1;
X = [S3 - Length, S1 + Length];
Y = [0, 0];

plot(X, Y);               %  Plot stress axis

Length = 1.1 * (S1 - S3) / 2
X = [0, 0];
Y = [-Length, Length];
plot(X, Y);               %  Plot shear axis

%  Change the aspect ratio so that the circles are circular

set(gca,'DataAspectRatio',[1 1 1]);
hold off

end

function Circle(x,y,r)
%x and y are the coordinates of the center of the circle
%r is the radius of the circle
%0.01 is the angle step, bigger values will draw the circle faster but
%you might notice imperfections (not very smooth)
ang=0:0.1:2*pi; 
xp=r*cos(ang);
yp=r*sin(ang);
plot(x+xp,y+yp);
end


