Sources Esope | Cast3M

pdirec
C PDIREC    SOURCE    CB215821  17/07/21    21:15:23     9513           
c---------------------------------------------------------------------
c
      SUBROUTINE PDIREC (NDIME, NSTR1, SIGMA, SGPRI, BETAM)
c
c=====================================================================
c                                                                    =
c   This routine computes the eigenvectors of a structural tensor.   =
c                                                                    =
c    Note:   sigma = (Sxx, Syy, Szz, Sxy, Sxz, Syz)                  =
c            sgpri = (S11, S22, S33)         for NSTR1=6             =
c            sgpri = (S11, S22, Szz)         for NSTR1=4             =
c                                                                    =
c=====================================================================
      IMPLICIT INTEGER(I-N)
      integer   ndime, nstr1
      real*8    betam (ndime,*), sigma (*), sgpri (*)
c
      integer   i
      real*8    twopi, error, str11, str22, str12, denta,
     .          angle, ca   , ca2  , sa   , sa2  , value,
     .          dummy (3,3) , absnor
c
      data twopi, error / 6.283185307179586 d0, 1.0 d-12 /
c
      if (ndime .eq. 2)  then
          str11 = sigma (1)
          str22 = sigma (2)
          str12 = sigma (4)
c
          absnor=sqrt(str11**2 + str22**2 + str12**2)
c
          denta = str11 - str22
c
          if (abs (str12)/absnor .le. error)  then
              angle = 0.0 d0
              if (denta .lt. 0.0 d0)  angle = twopi * 0.25 d0
          else
              if (abs (denta)/absnor .le. error)  then
                  angle = twopi * 0.125 d0
                  if (str12 .lt. 0.0 d0)  angle = -angle
              else
                  angle = 0.5 d0 * atan (2.0d0*str12/abs(denta))
              end if
          end if
c
          ca    = cos (angle)
          sa    = sin (angle)
          ca2   = ca **2
          sa2   = sa **2
          value = 2.0 d0 * sa * ca * str12
c
          sgpri (1) = ca2 * str11 + sa2 * str22 + value
          sgpri (2) = sa2 * str11 + ca2 * str22 - value
          sgpri (3) = sigma (3)
c
          if (sgpri (1) .lt. sgpri (2))  then
c
              value =  ca
              ca    = -sa
              sa    =  value
c
              value = sgpri (1)
              sgpri (1) = sgpri (2)
              sgpri (2) = value
          end if
c=====----------------------------------------------------------------
c  build transformation matrix                                       |
c=====----------------------------------------------------------------
          betam (1,1) =  ca
          betam (2,2) =  ca
          betam (2,1) = -sa
          betam (1,2) =  sa
      else
          if (sgpri (1) .eq. 0.0 d0 .and. sgpri (2) .eq. 0.0 d0
     .                              .and. sgpri (3) .eq. 0.0 d0)  then
              call zero  (betam, 3, 3)
              do 10 i = 1, 3
                 betam (i,i) = 1.0 d0
   10         continue
              return
          end if
c
          call pdire3 (sigma, sgpri, dummy)
c
          betam (1,1) = dummy (1,1)
          betam (1,2) = dummy (2,1)
          betam (1,3) = dummy (3,1)
          betam (2,1) = dummy (1,2)
          betam (2,2) = dummy (2,2)
          betam (2,3) = dummy (3,2)
          betam (3,1) = dummy (1,3)
          betam (3,2) = dummy (2,3)
          betam (3,3) = dummy (3,3)
      end if
c
      return
      end