Esope source files

   1 : C PVECTO    SOURCE    CHAT      05/01/13    02:39:59     5004           
   2 : c---------------------------------------------------------------------
   3 : c
   4 :       SUBROUTINE PVECTO (S, XL, V, COM)
   5 : c
   6 : c=====================================================================
   7 : c                                                                    =
   8 : c   This routine calculates one eigenvector.                         =
   9 : c                                                                    =
  10 : c    Input:    xl     eigenvalue                                     =
  11 : c              s (6)  original matrix                                =
  12 : c    Output:   v (3)  eigenvector                                    =
  13 : c                                                                    =
  14 : c    Note:   s = (Sxx, Syy, Szz, Sxy, Sxz, Syz)                      =
  15 : c                                                                    =
  16 : c=====================================================================
  17 :       IMPLICIT INTEGER(I-N)
  18 :       real*8    s (6), v (3), xl, com
  19 : c
  20 :       real*8    tol, x1, x2, x3, x4, x5, x6, xn, det, deti, v1, v2
  21 : c
  22 :       data tol  / 1.0 d-07 / 
  23 : c
  24 :       x1 = s (1) - xl
  25 :       x2 = s (2) - xl
  26 :       x3 = s (3) - xl
  27 :       x4 = s (4)
  28 :       x5 = s (6)
  29 :       x6 = s (5)
  30 : c
  31 :       det = x1 * x2 - x4 **2
  32 :       if (dabs (det/com) .gt. tol)  then
  33 :         deti = 1.0 d0 / det
  34 :         v1 = (x4 * x5 - x6 * x2) * deti
  35 :         v2 = (x6 * x4 - x5 * x1) * deti
  36 :         xn = 1.0 d0 / dsqrt (v1**2+v2**2+1.0d0)
  37 :         v (1) = v1 * xn
  38 :         v (2) = v2 * xn
  39 :         v (3) = xn
  40 :         return
  41 :       end if
  42 : c
  43 :       det = x4 * x5 - x6 * x2
  44 :       if (dabs (det/com) .gt. tol)  then
  45 :         v1 = (x6 * x4 - x5 * x1) / det
  46 :         xn = 1.0 d0 / dsqrt (v1**2+1.0d0)
  47 :         v (1) = xn
  48 :         v (2) = v1 * xn
  49 :         v (3) = 0.0 d0
  50 :         return
  51 :       end if
  52 : c
  53 :       det = x1 * x5 - x4 * x6
  54 :       if (dabs (det/com) .gt. tol)  then
  55 :         v (1) = 0.0 d0
  56 :         v (2) = 1.0 d0
  57 :         v (3) = 0.0 d0
  58 :         return
  59 :       end if
  60 : c
  61 :       det = x4 * x3 - x5 * x6
  62 :       if (dabs (det/com) .gt. tol)  then
  63 :         v1 = (x5 **2 - x2 * x3) / det
  64 :         xn = 1.0 d0 / dsqrt (v1**2+1.0d0)
  65 :         v (1) = v1 * xn
  66 :         v (2) = xn
  67 :         v (3) = 0.0 d0
  68 :         return
  69 :       end if
  70 : c
  71 :       det = x2 * x3 - x5 **2
  72 :       if (dabs (det/com) .gt. tol)  then
  73 :         v (1) = 1.0 d0
  74 :         v (2) = 0.0 d0
  75 :         v (3) = 0.0 d0
  76 :         return
  77 :       end if
  78 : c
  79 :       v (1) = 0.0 d0
  80 :       v (2) = 1.0 d0
  81 :       v (3) = 0.0 d0
  82 : c
  83 :       end
  84 :  
  85 :