clear
sig_x=10
sig_y=-10
tau_xy=-20

% sigma_average, R
sigma_ave=(sig_x+sig_y)/2
R=sqrt((0.5*(sig_x-sig_y))^2+tau_xy^2)
%**************************************************
% principal stress calculations
%**************************************************
thetap=0.5*atan((2*tau_xy)/(sig_x-sig_y));
thetap_in_deg=thetap*180/pi


sig_1=sigma_ave+R
sig_2=sigma_ave-R

 % theta_p corresponds to which principal stress plane (sigma_1 or sigma_2?
sigma_xp=sig_x*cos(thetap)^2+sig_y*sin(thetap)^2+2*tau_xy*sin(thetap)*cos(thetap)

% if theta_p corresponds to sigma_1 then call it theta_p1
% if theta_p corresponds to sigma_2 then call it theta_p2

%**************************************************
% Maxm shearing stress calculations
%**************************************************
thetas=0.5*atan(-(sig_x-sig_y)/(2*tau_xy));
thetas_in_deg=thetas*180/pi
taumax=+R
taumin=-R

 % theta_s corresponds to which shearing stress plane (max or min)?
tau_xyp=tau_xy*(cos(thetas)^2-sin(thetas)^2)+(sig_y-sig_x)*sin(thetas)*cos(thetas)
% if theta_s corresponds to taumax then call it theta_s1
% if theta_s corresponds to taumin then call it theta_s2