function [f,g,h,i]=tau_principal(tau_xx,tau_yy,tau_xy)
%TAU_PRINCIPAL - compute principal component of stress tensor
%
%   Usage:
%      [f,g,h,i]=tau_principal(tau_xx,tau_yy,tau_xy)

nel=length(tau_xx);
tau_xx=tau_xx/(10^5);
tau_yy=tau_yy/(10^5);
tau_xy=tau_xy/(10^5);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%determination of principal components of stress tensor%%%%%%%%%%%%
teta1=zeros(nel,1);
teta2=zeros(nel,1);
ind=find((tau_xx-tau_yy)~=0);
teta1(ind)=1/2*atan(tau_xy(ind)./(tau_xx(ind)-tau_yy(ind)));
ind=find((tau_xx-tau_yy)==0);
teta1(ind)=pi/2;
teta2=teta1+pi/2;

tmp1=tau_xx.*cos(teta1).^2+tau_yy.*sin(teta1).^2+2*tau_xy.*sin(teta1).*cos(teta1);
tmp2=tau_xx.*cos(teta2).^2+tau_yy.*sin(teta2).^2+2*tau_xy.*sin(teta2).*cos(teta2);

a=tmp1<=tmp2;
b=tmp1>=tmp2;
tau_1=tmp1.*b+a.*tmp2;
tau_2=tmp1.*a+b.*tmp2;

teta1_pp=teta1.*b+a.*teta2;
teta2_pp=teta1.*a+b.*teta2;
%teta1_pp=mod(teta1_pp,2*pi);
%teta2_pp=mod(teta2_pp,2*pi);
f=tau_1;
g=tau_2;
h=teta1_pp;
i=teta2_pp;

