
%  stress_tensor_2D_angle.m  ver 1.0  by Tom Irvine

function[pstress,pangle]=stress_tensor_2D_angle(sx,sy,txy)

ss=size(sx);

NT=ss(1);

pstress=zeros(NT,2);
pangle=zeros(NT,1);

eps0 = 1e-9;
    
for i=1:NT
        
    q1=(sx(i)+sy(i))/2;
    q2=sqrt( ((sx(i)-sy(i))/2).^2 + txy(i).^2 );  
    
    s1= q1+q2;
    s2= q1-q2; 
    
    pstress(i,:)=[s1 s2];

    d=sx(i)-sy(i);

    if(abs(d)<eps0 && abs(txy(i))<eps0)
        % This is a hydrostatic or isotropic stress state 
        % (equal normal stresses, no shear). 
        % Every direction is a principal direction
        pangle(i)=0;    
    else    
        pangle(i)=atan2(2*txy(i),d)/2;
    end
end

pangle = rad2deg(unwrap(pangle));