[Key,Smax,Smin,Savg,Tmxmn, ...

function [Key,Smax,Smin,Savg,Tmxmn, ...
  Smax_ang,Smin_ang,Tmax_ang,Tmin_ang]= ...
  prstress(Sigma_x,Sigma_y,Tau_xy,epsilon)
%
% [Key,Smax,Smin,Savg,Tmxmn,Smax_ang, ...
%  Smin_ang,Tmax_ang,Tmin_ang]= ...
%  prstress(Sigma_x,Sigma_y,Tau_xy,epsilon)
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% This function calculates the principal normal 
% and shear stresses for an element.
%
% Sigma_x  - sigma in the x direction
% Sigma_y  - sigma in the y direction
% Tau_xy   - tau in the xy plane
% epsilon  - tolerance for determining
%            if two numbers are equal
%            (optional)
%
% Key      - =0, normal exit
%            =1, all axes are principal
% Smax     - Sigma maximum
% Smin     - Sigma minimum
% Savg     - Sigma on planes of max/min
%            shear stress
% Tmxmn    - Tau maximum. Tau minimum is 
%            negative of Tau maximum
% Smax_ang - angle to Sigma_max
% Smin_ang - angle to Sigma_min
% Tmax_ang - angle to Tau_max
% Tmin_ang - angle to Tau_min
%
% NOTES:  1)  All angles measured CCW from
%             positive x-axis in degrees.
%         2)  Normal stresses are + for
%             tension, - for compression.
%         3)  Shear stresses are + for
%             up on right face.
%
% User m functions called:  strangle
%----------------------------------------------

%...Initialize
raddeg=180/pi; Key=0;

%...Max/min values
Savg=(Sigma_x+Sigma_y)/2;
Sdiff=(Sigma_x-Sigma_y)/2;
Tmxmn=sqrt(Sdiff^2+Tau_xy^2);
Smax=Savg+Tmxmn; Smin=Savg-Tmxmn;

if nargin ~= 4
  if Smax == Smin
    epsilon=1e-3;
  else
    epsilon=abs(Smax-Smin)*1e-3;
  end
end

if abs(Smax-Smin) < epsilon
  %...Hydrostate case, all axes are principal
  Key=1;
  Smax=Sigma_x; Smin=Sigma_x; Tmxmn=0;
  Smax_ang=0; Smin_ang=0;
  Tmax_ang=0; Tmin_ang=0;
elseif abs(Sigma_x-Sigma_y) < epsilon
  %...In state of max/min shear
  if Tau_xy < 0
    Smax_ang=45; Smin_ang=-45;
    Tmax_ang=90; Tmin_ang=0;
  else
    Smax_ang=-45; Smin_ang=45;
    Tmax_ang=0; Tmin_ang=90;
  end
elseif abs(Tau_xy) < epsilon
  %...In state of max/min normal stress
  if Sigma_x > Sigma_y
    %...sigma-x is max
    Smax_ang=0; Smin_ang=90;
    Tmax_ang=45; Tmin_ang=-45;
  else
    %...sigma-y is max
    Smax_ang=90; Smin_ang=0;
    Tmax_ang=-45; Tmin_ang=45;
  end
else
  %...Get angle to Smax plane
  angle=0.5*atan2(-2*Tau_xy,Sigma_x-Sigma_y);
  Smax_ang=angle*raddeg;
  Smin_ang=(angle+pi/2)*raddeg;
  Tmax_ang=Smax_ang+45;
  Tmin_ang=Tmax_ang+90;

  %...Make angles between -90 and +90 
  Smax_ang=flpang(Smax_ang,0);
  Smin_ang=flpang(Smin_ang,0);
  Tmax_ang=flpang(Tmax_ang,0);
  Tmin_ang=flpang(Tmin_ang,0);
end