!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!                                                                              |
!              ||===\\                                                         | 
!              ||    \\                                                        |
!              ||     ||   //==\\   ||  ||   //==||  ||/==\\                   |
!              ||     ||  ||    ||  ||  ||  ||   ||  ||    ||                  |
!              ||    //   ||    ||  ||  ||  ||   ||  ||                        |
!              ||===//     \\==//    \\==\\  \\==\\  ||                        |
!                                                                              |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!                                                                              |
!              ROUTINE    Nov. 2006                                            |
!                                                                              |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|

subroutine PS (cc11,cc22,cc33,cc12,cc23,cc31,l1,l2,l3, &
               n11,n12,n13, &
               n21,n22,n23, &
               n31,n32,n33)

implicit none

!------------------------------------------------------------------------------|
!(((((((((((((((( Purpose of the routine  ))))))))))))))))))))))))))))))))))))))
!------------------------------------------------------------------------------|
! calculates principal directions and eigenvalues of symmetrical 
! tensor cc11,cc22,cc33,cc12,cc23,cc31

!------------------------------------------------------------------------------|
!((((((((((((((((  declaration of the subroutine arguments  ))))))))))))))))))))
!------------------------------------------------------------------------------|

double precision cc11,cc22,cc33,cc12,cc23,cc31
double precision l1,l2,l3
double precision n11,n12,n13,n21,n22,n23,n31,n32,n33

!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
!------------------------------------------------------------------------------|

double precision a1,a2,a3,q
double precision r, theta, a13, xnorm
double precision q1,q2,q3 
double precision x1,x2,x3 
double precision pi3,pi23,pi43
double precision c11,c22,c33,c12,c23,c13
double precision mi(3),mimax
integer im

!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------

pi3=3.141592654d0/3.d0
pi23=pi3*2.d0
pi43=pi3*4.d0

c11=cc11
c22=cc22
c33=cc33
c12=cc12
c23=cc23
c13=cc31

! find the stress tensor eigenvalues (principal stresses)
  
a3=-c11*c22*c33-c12*c23*c13-c12*c23*c13 &
   +c13*c22*c13+c12*c12*c33+c11*c23*c23
a2=c11*c22+c22*c33+c33*c11-c13*c13-c23*c23-c12*c12
a1=-c11-c22-c33
  
q=(a1*a1-3.d0*a2)/9.d0
r=(2.d0*a1*a1*a1-9.d0*a1*a2+27.d0*a3)/54.d0
  
q3=q*q*q

! check for multiple eigenvalues
                                 
if (q3-r*r.le.0.d0) then

   l1=c11
   l2=c22
   l3=c33

else
  
   theta=acos(r/sqrt(q3))/3.d0
   q2=sqrt(q)
   a13=a1/3.d0
   l1=-2.d0*q2*cos(theta)-a13
   l2=-2.d0*q2*cos(theta+pi23)-a13
   l3=-2.d0*q2*cos(theta+pi43)-a13

endif

! first eigenvector

c11=c11-l1
c22=c22-l1
c33=c33-l1
mi(1)=c22*c33-c23*c23
mi(2)=c33*c11-c13*c13
mi(3)=c11*c22-c12*c12

im=2
if (abs(mi(1)).gt.abs(mi(2))) im=1
if (abs(mi(im)).lt.abs(mi(3))) im=3
                               
mimax=mi(im)

if (im.eq.1) then
   x2=(c13*c23-c12*c33)/mimax
   x3=(c12*c23-c13*c22)/mimax
   xnorm=sqrt(1.d0+x2*x2+x3*x3)
   n11=1.d0/xnorm
   n12=x2/xnorm
   n13=x3/xnorm
elseif (im.eq.2) then
   x3=(c12*c13-c23*c11)/mimax
   x1=(c13*c23-c12*c33)/mimax
   xnorm=sqrt(1.d0+x3*x3+x1*x1)
   n11=x1/xnorm
   n12=1.d0/xnorm
   n13=x3/xnorm
elseif (im.eq.3) then
   x1=(c23*c12-c13*c22)/mimax
   x2=(c13*c12-c23*c11)/mimax
   xnorm=sqrt(1.d0+x1*x1+x2*x2)
   n11=x1/xnorm
   n12=x2/xnorm
   n13=1.d0/xnorm
endif

! second eigenvector

              c11=c11+l1-l2
              c22=c22+l1-l2
              c33=c33+l1-l2
              mi(1)=c22*c33-c23*c23
              mi(2)=c33*c11-c13*c13
              mi(3)=c11*c22-c12*c12
              
              im=2
              if (abs(mi(1)).gt.abs(mi(2))) im=1
              if (abs(mi(im)).lt.abs(mi(3))) im=3
                                     
              mimax=mi(im)
                if (im.eq.1) then
                x2=(c13*c23-c12*c33)/mimax
                x3=(c12*c23-c13*c22)/mimax
                xnorm=sqrt(1.d0+x2*x2+x3*x3)
                n21=1.d0/xnorm
                n22=x2/xnorm
                n23=x3/xnorm
                  elseif (im.eq.2) then
                  x3=(c12*c13-c23*c11)/mimax
                  x1=(c13*c23-c12*c33)/mimax
                  xnorm=sqrt(1.d0+x3*x3+x1*x1)
                  n21=x1/xnorm
                  n22=1.d0/xnorm
                  n23=x3/xnorm
                    elseif (im.eq.3) then
                    x1=(c23*c12-c13*c22)/mimax
                    x2=(c13*c12-c23*c11)/mimax
                    xnorm=sqrt(1.d0+x1*x1+x2*x2)
                    n21=x1/xnorm
                    n22=x2/xnorm
                    n23=1.d0/xnorm
                    endif
  
! third eigenvector

c11=c11+l2-l3
c22=c22+l2-l3
c33=c33+l2-l3
mi(1)=c22*c33-c23*c23
mi(2)=c33*c11-c13*c13
mi(3)=c11*c22-c12*c12
              
im=2
if (abs(mi(1)).gt.abs(mi(2))) im=1
if (abs(mi(im)).lt.abs(mi(3))) im=3
                                     
mimax=mi(im)

if (im.eq.1) then
   x2=(c13*c23-c12*c33)/mimax
   x3=(c12*c23-c13*c22)/mimax
   xnorm=sqrt(1.+x2*x2+x3*x3)
   n31=1.d0/xnorm
   n32=x2/xnorm
   n33=x3/xnorm
elseif (im.eq.2) then
   x3=(c12*c13-c23*c11)/mimax
   x1=(c13*c23-c12*c33)/mimax
   xnorm=sqrt(1.d0+x3*x3+x1*x1)
   n31=x1/xnorm
   n32=1.d0/xnorm
   n33=x3/xnorm
elseif (im.eq.3) then
   x1=(c23*c12-c13*c22)/mimax
   x2=(c13*c12-c23*c11)/mimax
   xnorm=sqrt(1.d0+x1*x1+x2*x2)
   n31=x1/xnorm
   n32=x2/xnorm
   n33=1.d0/xnorm
endif

return
end

!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------