C PVECTO    SOURCE    CHAT      05/01/13    02:39:59     5004           
c---------------------------------------------------------------------
c
      SUBROUTINE PVECTO (S, XL, V, COM)
c
c=====================================================================
c                                                                    =
c   This routine calculates one eigenvector.                         =
c                                                                    =
c    Input:    xl     eigenvalue                                     =
c              s (6)  original matrix                                =
c    Output:   v (3)  eigenvector                                    =
c                                                                    =
c    Note:   s = (Sxx, Syy, Szz, Sxy, Sxz, Syz)                      =
c                                                                    =
c=====================================================================
      IMPLICIT INTEGER(I-N)
      real*8    s (6), v (3), xl, com
c
      real*8    tol, x1, x2, x3, x4, x5, x6, xn, det, deti, v1, v2
c
      data tol  / 1.0 d-07 / 
c
      x1 = s (1) - xl
      x2 = s (2) - xl
      x3 = s (3) - xl
      x4 = s (4)
      x5 = s (6)
      x6 = s (5)
c
      det = x1 * x2 - x4 **2
      if (dabs (det/com) .gt. tol)  then
        deti = 1.0 d0 / det
        v1 = (x4 * x5 - x6 * x2) * deti
        v2 = (x6 * x4 - x5 * x1) * deti
        xn = 1.0 d0 / dsqrt (v1**2+v2**2+1.0d0)
        v (1) = v1 * xn
        v (2) = v2 * xn
        v (3) = xn
        return
      end if
c
      det = x4 * x5 - x6 * x2
      if (dabs (det/com) .gt. tol)  then
        v1 = (x6 * x4 - x5 * x1) / det
        xn = 1.0 d0 / dsqrt (v1**2+1.0d0)
        v (1) = xn
        v (2) = v1 * xn
        v (3) = 0.0 d0
        return
      end if
c
      det = x1 * x5 - x4 * x6
      if (dabs (det/com) .gt. tol)  then
        v (1) = 0.0 d0
        v (2) = 1.0 d0
        v (3) = 0.0 d0
        return
      end if
c
      det = x4 * x3 - x5 * x6
      if (dabs (det/com) .gt. tol)  then
        v1 = (x5 **2 - x2 * x3) / det
        xn = 1.0 d0 / dsqrt (v1**2+1.0d0)
        v (1) = v1 * xn
        v (2) = xn
        v (3) = 0.0 d0
        return
      end if
c
      det = x2 * x3 - x5 **2
      if (dabs (det/com) .gt. tol)  then
        v (1) = 1.0 d0
        v (2) = 0.0 d0
        v (3) = 0.0 d0
        return
      end if
c
      v (1) = 0.0 d0
      v (2) = 1.0 d0
      v (3) = 0.0 d0
c
      end