Mohr's Circle

%Variables
sigx = -75;
sigy = 60;
tauxy = 45;

%Euation of a circle is: (x - xc)^2 + (y - yc)^2 = R^2

%Computing principal stress
R = sqrt(((sigx - sig y)/2)^2 + tauxy^2);
C = (sigx + sigy)/2;
sig1 = C + R; %Largest principal stress
Sig2 = C - R; %Smallest principal stress

%Define variable for normal stress
sigma = linspace(sig2, sig1, 1001); %Vectoring function

%Evaluate shear stress. Use + and - square root
tau1 = sqrt(R^2 - (sigma - C).^2);
tau2 = -tau1;
plot(sigma, tau1, 'b', sigma, tau2, 'b');
xlabel('Sigma');

%Fix the aspect ratio so the circle would look like a
%circle and not an ellipse
axis('square');

%Add in diamater for principal stresses
hold on;
plot([sig1, sig2], [0, 0], 'k');

%Add in diamater for initial stress state
plot([sigx, sigy], [-tauxy, tauxy], 'g');

%(Optional) Grid helps for visualization
grid on;