c
c........|.........|.........|.........|.........|.........|.........|..
c........|.........|.........|.........|.........|.........|.........|..
c
c Rmatrix and associated programs may be used freely by University
c researchers and teachers. Commercial parties interested in using
c the program should contact C. J. Thomson.
c
c........|.........|.........|.........|.........|.........|.........|..
c........|.........|.........|.........|.........|.........|.........|..
c
c
      subroutine aniso6c1(a,rho,px,py,eval,evec)
      implicit real*8 (a-h,o-z)

c Copyright (c) 1996 C. J. Thomson.
c All rights reserved by the author.

c Finds e-vals and e-vecs for given layer.
c Based on Sean Guest's routines for finding vertical slowness given
c horizontal slownesses by solving the 6X6 matrix problem (using
c public domain library routines - Eispack? They look like it.).
c
      dimension a(3,3,3,3),p(3)
      complex*16 px,py
      complex*16 eval(6),evec(6,6),xnorm
      logical iso
      common/info/lay,id,iso

      iso=.false.

      p(1)=dble(px)
      p(2)=dble(py)
      call eig (rho,a,p,eval,evec)

c Normalise these evecs to unit displacement

      do 100 j=1,6
        xnorm=cdsqrt(evec(1,j)**2+evec(2,j)**2+evec(3,j)**2)
        do 101 i=1,6
          evec(i,j)=evec(i,j)/xnorm
  101   continue
  100 continue

      return
      end
c
c........|.........|.........|.........|.........|.........|.........|..
c
      subroutine eig (den,ae,p,wc,zc)

c Copyright (c) 1996 W. S. Guest and C. J. Thomson.
c All rights reserved by the authors.

c
c          W.S. Guest  October 1988
c
c This subroutine finds the eigenvalues (vertical slowness) and
c eigenvectors (displacement-stress) when given the 81 component
c elastic constants AE, the density DEN, and the horizontal slownesses
c p(1) and p(2). If the medium is isotropic, then a call is made
c to ISOEIG which provide a rapid analytic solution (Fryer & Frazer).
c THe anisotropic formulation is after Woodhouse (1974) and then the
c system matrix is reduced using the `numerical recipes' routine DEIGV.
c The vectors must then be sorted into
c                 Pu,SVu,SHu,Pd,SVd,SHd
c This is done in SORT which bases the S orientation on polarization
c or in SORT_AMP which sorts the eigenvalues from smallest to largest.
c
      implicit real*8(a-h,o-z)
c
      integer i,j,k,l,m,n,lay,id
      logical iso
      dimension c(3,3,3,3),ae(3,3,3,3),cm(3,3),ci(3,3),t(3,3),tt(3,3)
      dimension ct2(3,3,3),ct(3,3),qc(3,3,3,3),s(3,3),a(6,6),p(3)
      dimension zi(6,6),wr(6),wi(6),zr(6,6)
      complex*16 wc(6),zc(6,6)
c Common Block
      common/info/lay,id,iso
c
c
c     write (*,*) 'ineig',iso
      do 10 i=1,3
       do 20 j=1,3
        do 30 k=1,3
         do 40 l=1,3
           c(i,j,k,l) = ae(i,j,k,l) * den
 40      continue
 30     continue
 20    continue
 10   continue
c
      if (iso) then
c        write (*,*)'IsoEig used to solve "A"'
         c33=c(3,3,3,3)
         c44=c(2,3,2,3)
c        write (*,*) c33,c44
         call isoeig(den,p,c33,c44,wc,zc)
         goto 999
      endif
c
c Calculate the c33 inverse matrix
c
      cm(1,1)=c(2,3,2,3)*c(3,3,3,3) - c(2,3,3,3)**2.d0
      cm(2,2)=c(1,3,1,3)*c(3,3,3,3) - c(1,3,3,3)**2.d0
      cm(3,3)=c(1,3,1,3)*c(2,3,2,3) - c(1,3,2,3)**2.d0
      cm(1,2)=-1.d0*(c(2,3,1,3)*c(3,3,3,3) - c(2,3,3,3)*c(3,3,1,3))
      cm(1,3)=c(2,3,1,3)*c(3,3,2,3) - c(2,3,2,3)*c(3,3,1,3)
      cm(2,3)=-1.d0*(c(1,3,1,3)*c(3,3,2,3) - c(1,3,2,3)*c(3,3,1,3))
      cm(2,1)=cm(1,2)
      cm(3,1)=cm(1,3)
      cm(3,2)=cm(2,3)
      det=c(1,3,1,3)*cm(1,1) + c(2,3,1,3)*cm(2,1) + c(3,3,1,3)*cm(3,1)
c
      do 50 i=1,3
        do 55 j=1,3
           ci(i,j)=cm(i,j)/det
 55     continue
 50   continue
c
c Calculate the 'T' and 'T' transpose matrices
c
      do 60 m=1,3
        do 65 n=1,3
          c1=0.d0
          c2=0.d0
          do 70 k=1,3
            c1 = c1 + c(m,1,k,3)*ci(k,n)
            c2 = c2 + c(m,2,k,3)*ci(k,n)
 70       continue
          tt(n,m) = -1.d0*(p(1)*c1 + p(2)*c2)
          t(m,n) = tt(n,m)
 65     continue
 60   continue
c
c Calculate the 'Q' matrix
c
      do 75 i=1,2
        do 80 k=1,3
          do 85 l=1,3
            ct(k,l)=0.d0
            do 90 m=1,3
              ct(k,l) = ct(k,l) + c(k,i,m,3)*ci(m,l)
 90         continue
 85       continue
 80     continue
c
        do 95 j=1,2
          do 100 k=1,3
            do 105 l=1,3
              ct2(j,k,l)=0.d0
              do 110 m=1,3
                ct2(j,k,l)= ct2(j,k,l) + ct(k,m)*c(m,3,l,j)
 110          continue
              qc(i,j,k,l)=c(k,i,l,j) - ct2(j,k,l)
 105        continue
 100      continue
 95     continue
 75   continue
c
c Calculate the 'S' matrix
c
      do 115 k=1,3
        do 120 l=1,3
          s(k,l)=0.d0
          do 125 i=1,2
            do 130 j=1,2
              s(k,l) = s(k,l) + p(i)*p(j)*qc(i,j,k,l)
 130        continue
 125      continue
          s(k,l) = -1.d0 *s(k,l)
          if (k.eq.l) s(k,l)=s(k,l)+den
 120    continue
 115  continue
c
c Form the complete system matrix 'A'
c
      do 135 i=1,3
        j=i+3
        do 140 k=1,3
          l=k+3
          a(i,k) =  tt(i,k)
          a(i,l) =  ci(i,k)
          a(j,k) =  s(i,k)
          a(j,l) =  t(i,k)
 140    continue
 135  continue
c
c Find the eigenvalues and eigenvectors
c
      call deigv(6,a,wr,wi,zr,zi)
c
      call sort_pol(wr,wi,zr,zi,wc,zc,ae,p)
c     call sort_amp(wr,wi,zr,zi,wc,zc,ae,p)
c
 999  continue
c
c      write (19,*)
c      write (19,'(6e12.3)') wc(1),wc(2),wc(3),wc(4),wc(5),wc(6)
c      write (19,*)
c      write (19,'(12d10.2)')  zc(1,1),zc(1,2),zc(1,3),
c    *                       zc(1,4),zc(1,5),zc(1,6)
c      write (19,'(12d10.2)')  zc(2,1),zc(2,2),zc(2,3),
c    *                       zc(2,4),zc(2,5),zc(2,6)
c      write (19,'(12d10.2)')  zc(3,1),zc(3,2),zc(3,3),
c    *                       zc(3,4),zc(3,5),zc(3,6)
c
      return
      end
c
c........|.........|.........|.........|.........|.........|.........|..
c
      subroutine sort_amp (wr,wi,zr,zi,wc,zc,ae,p)

c Copyright (c) 1996 W. S. Guest and C. J. Thomson.
c All rights reserved by the authors.

c
c same as SORT except that the shear eigenvalue separation is based
c soley on magnitudes. If all are real then the eigenvalues are arranged
c as 1, the smallest magnitude, to 3, the largest. This means that 3 will
c be the last to go evanescant. When all three waves are imaginary then the
c magnitude of the imaginary components are monitored and they are arranged
c from 1, the largest, to 3, the smallest. If there are no intersections then
c this will guarentee smooth flowing curves.
c
      implicit real*8 (a-h,o-z)
c
      integer np(6),neg(3),pos(3),neg_ev(3),pos_ev(3),in_ev,ip_ev,
     +        nvmax,pvmax,nvmin,pvmin,in,ip,i,j,k
      dimension wr(6),wi(6),zr(6,6),zi(6,6),p(3),v(3),ae(3,3,3,3)
      complex*16 wc(6),zc(6,6)
c Common Blocks
      common/info/lay,id,iso
c
c
      numreal=0
      in=0
      ip=0
      in_ev=0
      ip_ev=0
      velmaxn=0.d0
      velmaxp=0.d0
      velminn=10000.d0
      velminp=10000.d0
      holdn=-10.d0
      holdp=10.d0
c
      do i=1,6
         np(i)=0
      enddo
c
c      write (19,'(6e10.3)') wr
c      write (19,'(6e10.3)') wi
c
c
c
      do i=1,6
        if (dabs(wi(i)) .lt. 1.d-9)then
          numreal=numreal+1
          p(3)=wr(i)
          call vtot(ae,p,v)
          velocity=dsqrt(v(1)**2. + v(2)**2. + v(3)**2.)
c          write (19,'(i4,4f10.1)')i,v,velocity
c
          if (v(3) .lt. 0.d0)then
            in=in+1
            neg(in)=i
            if (velocity .lt. velminn)then
               velminn=velocity
               nvmin=i
            endif
            if (velocity .gt. velmaxn)then
               velmaxn=velocity
               nvmax=i
            endif
          else
            ip=ip+1
            pos(ip)=i
            if (velocity .lt. velminp)then
               velminp=velocity
               pvmin=i
            endif
            if (velocity .gt. velmaxp)then
               velmaxp=velocity
               pvmax=i
            endif
          endif
        elseif (wi(i) .lt. 0.d0)then
          in=in+1
          in_ev=in_ev+1
          neg_ev(in_ev)=i
          neg(in)=i
        else
          ip=ip+1
          ip_ev=ip_ev+1
          pos_ev(ip_ev)=i
          pos(ip)=i
        endif
       enddo
c
c      write (19,*) ip,in,pos,neg
c      write (19,*) 'number of real eigenvalues',numreal
c      write (19,*) pos_ev,neg_ev
c      write (19,*) 'pvmax,pvmin',pvmax,pvmin
c      write (19,*) 'nvmax,nvmin',nvmax,nvmin
c
      j=0
      k=0
c
      if (numreal .ge. 4)then
c
          np(3)=pvmin
          np(6)=nvmin
c
          if (numreal .eq. 6)then
            np(1)=pvmax
            np(4)=nvmax
          else
            np(1)=pos_ev(1)
            np(4)=neg_ev(1)
          endif
c
          do i=1,3
           if (pos(i).ne.np(1) .and. pos(i).ne.np(3)) np(2)=pos(i)
           if (neg(i).ne.np(4) .and. neg(i).ne.np(6)) np(5)=neg(i)
          enddo
c
      elseif (numreal .ge. 2)then
c
          np(3)=pvmin
          np(6)=nvmin
c
          if (wi(pos_ev(1)) .gt. wi(pos_ev(2)))then
            np(1)=pos_ev(1)
            np(2)=pos_ev(2)
          else
            np(2)=pos_ev(1)
            np(1)=pos_ev(2)
          endif
          if (wi(neg_ev(1)) .lt. wi(neg_ev(2)))then
            np(4)=neg_ev(1)
            np(5)=neg_ev(2)
          else
            np(5)=neg_ev(1)
            np(4)=neg_ev(2)
          endif
c
      else
c
          do 200 i=1,3
            if (wi(pos_ev(1)) .gt. wi(pos_ev(2))) then
              j=j+1
              if (j .eq. 1)then
                holdp=wi(pos(i))
                np(2)=pos(i)
              elseif (wi(pos(i)) .lt. holdp) then
                np(3)=np(2)
                np(2)=pos(i)
              else
                np(3)= pos(i)
              endif
            endif
c
            if (neg(i) .ne. np(4)) then
              k=k+1
              if (j .eq. 1)then
                holdn=wi(neg(i))
                np(5)=neg(i)
              elseif (wi(neg(i)) .gt. holdn) then
                np(6)=np(5)
                np(5)=neg(i)
              else
                np(6)= neg(i)
              endif
            endif
 200      continue
c
      endif
c
c     write (19,*) 'np',np
c
c Switch the eigenvalues and eigenvectors to the order P,S1,S2
c
      do i=1,6
        wc(i)=dcmplx(wr(np(i)),wi(np(i)))
        do j=1,3
          zc(j,i)=dcmplx(zr(j,np(i)),zi(j,np(i)))
c         zc(j+3,i)=dcmplx(zr(j+3,np(i)),zi(j+3,np(i)))*(0.d0,1.d0)
          zc(j+3,i)=dcmplx(zr(j+3,np(i)),zi(j+3,np(i)))
        enddo
      enddo
c
      return
      end
c
c........|.........|.........|.........|.........|.........|.........|..
c
      subroutine sort_pol (wr,wi,zr,zi,wc,zc,ae,p)

c Copyright (c) 1996 W. S. Guest and C. J. Thomson.
c All rights reserved by the authors.

c
c              W.S. Guest   July 1989
c
c This subroutine takes the input of the eigenvalues and eigenvectors
c calculated in eig and orders them into pu,svu,shu,pd,svd,shd. a
c check is first made to see if any of the eigenvalues are complex
c representing evanescant waves. separation of the up and downgoing
c waves is based on th sign of the eigenvalue but a check is made to
c make sure there are three of each sign, if not then the total vel.
c are calculated to get energy flux direction.
c if all of the eigenvalues are real,
c then the p-wave is taken as the smallest eigenvalues, both in the
c up and downgoing direction. the s-waves are not sorted, but
c always come out as sv and then sh. s-waves can be separated based
c on displacement, but problem occur for vertically travelling waves
c
      implicit real*8 (a-h,o-z)
c
      integer i,j,k,np(6),neg(3),npos(3)
      dimension wr(6),wi(6),zr(6,6),zi(6,6),
     +          p(3),v(3),ae(3,3,3,3)
      complex*16 wc(6),zc(6,6)
c Common Blocks
      common/info/lay,id,iso
c Parameter defintions
      parameter (pi = 3.1415926535897932d0)
c
      j=0
      k=0
      holdn=-10.d0
      holdp=10.d0
      numev=0
      hldevn=0.d0
      hldevp=0.d0
      do 5 i=1,6
         np(i)=0
 5    continue
c
c Check to see is any of the waves are evanesant
c
      do 10 i=1,6
        if (wi(i) .ne. 0.d0)then
          if (wi(i) .lt. 0.d0)then
            j=j+1
            neg(j)=i
            numev=numev+1
            if (wi(i) .lt. hldevn) then
              np(4)=i
              hldevn=wi(i)
            endif
          else
            numev=numev+1
            k=k+1
            npos(k)=i
            if (wi(i) .gt. hldevp) then
              np(1)=i
              hldevp=wi(i)
            endif
          endif
        else
          if (wr(i) .lt. 0.d0)then
            j=j+1
            neg(j)=i
            if (numev .eq. 0) then
              if (wr(i) .gt. holdn)then
                holdn=wr(i)
                np(4)=i
              endif
            endif
          else
            k=k+1
            npos(k)=i
            if (numev .eq. 0) then
              if (wr(i) .lt. holdp)then
                holdp=wr(i)
                np(1)=i
              endif
            endif
          endif
        endif
 10   continue
c
c     write (19,*) 'numev',numev,'   j',j,'   k',k
c
      if(j .gt. 3 .or. k .gt. 3)then
        j=0
        k=0
        do 20 i=1,6
          if (dabs(wi(i)) .lt. 1.d-9)then
            p(3)=-wr(i)
            call vtot(ae,p,v)
            if (v(3) .gt. 0.d0)then
              j=j+1
              neg(j)=i
            else
              k=k+1
              npos(k)=i
            endif
          elseif (wi(i) .lt. 0.d0)then
            j=j+1
            neg(j)=i
          else
            k=k+1
            npos(k)=i
          endif
 20     continue
      endif
c
      j=0
      k=0
      if (numev .le. 2)then
        do 755 i=1,3
          if (npos(i) .ne. np(1)) then
            j=j+1
            hfrac=dsqrt(dble(zr(1,npos(i))**2.d0)+dble
     *    (zr(2,npos(i)))**2.d0)/dsqrt(dble(zr(1,npos(i)))**2.d0
     *    +dble(zr(2,npos(i)))**2.d0+dble(zr(3,npos(i)))**2.d0)
            if (j .eq. 1) then
              hfrac1=hfrac
              ihld=i
            elseif (dabs(hfrac1) .gt. dabs(hfrac))then
              np(3)= npos(ihld)
              np(2)= npos(i)
            else
              np(3)= npos(i)
              np(2)= npos(ihld)
            endif
          endif
          if (neg(i) .ne. np(4)) then
            k=k+1
            hfracn=dsqrt(dble(zr(1,neg(i)))**2.d0+
     *    dble(zr(2,neg(i)))**2.d0)/dsqrt(dble(zr(1,neg(i)))**2.d0
     *    +dble(zr(2,neg(i)))**2.d0+dble(zr(3,neg(i)))**2.d0)
            if (k .eq. 1) then
              hfracn1=hfracn
              ihld=i
            elseif (dabs(hfracn1) .gt. dabs(hfracn))then
              np(6)= neg(ihld)
              np(5)= neg(i)
            else
              np(6)= neg(i)
              np(5)= neg(ihld)
            endif
          endif
 755    continue
c
      elseif (numev .le. 4)then
c
       do 756 i=1,3
c       if (npos(i) .ne. np(1)) then
c          j=j+1
c          if (j .eq. 1) then
c            if (wi(npos(i)) .gt. 0.d0)then
c              ihld=i
c            else
c              if (dabs(dble(zr(2,npos(i)))) .gt.
c     *            dabs(dble(zr(3,npos(i)))))then
c                 np(3)=npos(i)
c              else
c                 np(2)=npos(i)
c              endif
c            endif
c          else
c            if (wi(npos(i)) .gt. 0.d0)then
c              if (np(3) .eq. 0)then
c                np(3)=npos(i)
c              else
c                np(2)=npos(i)
c              endif
c            else
c              if (dabs(dble(zr(2,npos(i)))) .gt.
c     *            dabs(dble(zr(3,npos(i)))))then
c                 np(2)=npos(ihld)
c                 np(3)=npos(i)
c              else
c c                np(2)=npos(i)
c                 np(3)=npos(ihld)
c              endif
c            endif
c          endif
c         endif
        if (npos(i) .ne. np(1)) then
           if (wi(npos(i)) .gt. 0.d0)then
              np(5-id)=npos(i)
c              np(2)=neg(i)
           else
              np(id)=npos(i)
c              np(3)=neg(i)
           endif
        endif
c         if (neg(i) .ne. np(4)) then
c          k=k+1
c          if (k .eq. 1) then
c            if (wi(neg(i)) .gt. 0.d0)then
c              jhld=i
c            else
c              if (dabs(dble(zr(2,neg(i)))) .gt.
c     *            dabs(dble(zr(3,neg(i)))))then
c                 np(6)=neg(i)
c              else
c c                np(5)=neg(i)
c              endif
c            endif
c          else
c           if (wi(neg(i)) .gt. 0.d0)then
c              if (np(6) .eq. 0)then
c                np(6)=neg(i)
c              else
c                np(5)=neg(i)
c              endif
c           else
c              if (dabs(dble(zr(2,neg(i)))) .gt.
c     *            dabs(dble(zr(3,neg(i)))))then
c                 np(5)=neg(jhld)
c                 np(6)=neg(i)
c              else
c                 np(5)=neg(i)
c                 np(6)=neg(jhld)
c              endif
c            endif
c           endif
c          endif
c
       if (neg(i) .ne. np(4)) then
           if (wi(neg(i)) .lt. 0.d0)then
              np(8-id)=neg(i)
c              np(5)=neg(i)
           else
c              np(6)=neg(i)
              np(3+id)=neg(i)
           endif
        endif
756   continue
c
      elseif (numev .le. 6)then
c
       do 757 i=1,3
        if (npos(i) .ne. np(1)) then
           if (np(2) .eq. 0.d0)then
              np(2)=npos(i)
           else
              np(3)=npos(i)
           endif
        endif
        if (neg(i) .ne. np(4)) then
           if (np(5) .eq. 0.d0)then
              np(5)=neg(i)
           else
              np(6)=neg(i)
           endif
        endif
 757   continue
      endif
c
c     write (19,*) 'np',np
c
c Switch the eigenvalues and eigenvectors to the order P,S1,S2
c
      do i=1,6
        wc(i)=dcmplx(wr(np(i)),wi(np(i)))
        do j=1,3
          zc(j,i)=dcmplx(zr(j,np(i)),zi(j,np(i)))
c         zc(j+3,i)=dcmplx(zr(j+3,np(i)),zi(j+3,np(i)))*(0.d0,1.d0)
          zc(j+3,i)=dcmplx(zr(j+3,np(i)),zi(j+3,np(i)))
        enddo
      enddo
c
      return
      end
c
c........|.........|.........|.........|.........|.........|.........|..
c
      SUBROUTINE VTOT(CC,P,V)
C
      IMPLICIT REAL*8 (A-H,O-Z)
      integer lay,id
      logical iso
      DIMENSION CC(3,3,3,3),A(3,3),CP(3,3),DCP(3,3,3),V(3),P(3),DET(3)
      COMMON/info/lay,id,iso
C
      if (iso)then
c
        pmag=p(1)**2.d0+p(2)**2.d0+p(3)**2.d0
        do 10 i=1,3
           v(i)=p(i)/pmag
 10     continue
c
      else
c
        DO 28 I=1,3
          DO 27 L=1,3
             PSUM = 0.0D0
             DO 26 J=1,3
               DO 25 K=1,3
                  PSUM = PSUM+CC(I,J,K,L)*P(J)*P(K)
 25            CONTINUE
 26          CONTINUE
             CP(I,L) = PSUM
             CP(L,I) = PSUM
 27       CONTINUE
 28     CONTINUE
C
C
        DO 34 M=1,3
           do 33 i=1,3
              DO 32 L=1,3
                X = 0.0D0
                Y = 0.0D0
              DO 30 J=1,3
                 X = X+CC(I,J,M,L)*P(J)
 30           CONTINUE
              DO 31 K=1,3
                 Y = Y+CC(I,M,K,L)*P(K)
 31           CONTINUE
              DCP(M,I,L) = X+Y
 32        CONTINUE
 33      CONTINUE
 34    CONTINUE
C
C
       DO 45 I=1,3
          DET(I)=0.0D0
 45    CONTINUE
C
C
       DO 55 N=1,3
          DO 54 M=1,3
            DO 51 I=1,3
              DO 50 J=1,3
                 A(I,J)=CP(I,J)
 50           CONTINUE
 51         CONTINUE
            DO 52 I=1,3
               A(I,I) = CP(I,I)-1.D0
 52         CONTINUE
            DO 53 I=1,3
               A(I,M) = DCP(N,I,M)
 53         CONTINUE
            Z1 = A(1,1)*A(2,2)*A(3,3)-A(1,1)*A(2,3)*A(3,2)
            Z2 = -A(1,2)*A(2,1)*A(3,3)+A(1,2)*A(2,3)*A(3,1)
            Z3 = A(1,3)*A(2,1)*A(3,2)-A(1,3)*A(2,2)*A(3,1)
            DET(N) = DET(N)+Z1+Z2+Z3
 54      CONTINUE
 55    CONTINUE
C
       DEN = 0.0D0
       DO 65 I=1,3
         DEN = DEN + P(I)*DET(I)
 65    CONTINUE
C
       DO 71 I=1,3
         V(I) = DET(I)/DEN
 71    CONTINUE
c
      endif
C
       RETURN
       END
c
c........|.........|.........|.........|.........|.........|.........|..
c
      subroutine isoeig(den,p,c33,c44,wc,zc)

c Copyright (c) 1996 W. S. Guest and C. J. Thomson.
c All rights reserved by the authors.

c
c Original version: W.S. Guest   November 1988
c
c Modified/Corrected: C. Thomson   October 1993
c                     Sean's original lines commented out for
c                     comparison.
c  NOTE: c33 = lambda + 2mu  and  c44 = mu   in Lame notation
c
c
c Sean's comments follow...
c
c ``This subroutine calculates the displacement-stress eigenvectors
c and the eigenvalues in the case of an isotropic medium. The coding
c follows Fryer and Frazer (1987) with the obvious errors corrected.
c The vectors are choosen to be in a right handed system. P is along
c the slowness vector, SV is pointed upwards, and SH is horizontal
c but its direction depends on the value of p(x) (It goes to the
c left of the p(x) direction). All of the vectors are normalized so
c that the displacements are of unit magnitude.
c
c     Orientations rechecked   07/92     Sean Guest''
c
      implicit real*8(a-h,o-z)
c
      integer i
      dimension p(3)
      complex*16 wc(6),zc(6,6),norm,acomp,ccomp,bcomp
c
      p2=p(1)**2.d0 +p(2)**2.d0
      fac1=den/c33-p2
      fac2=den/c44-p2
c     write (19,*) p2,fac1,fac2
c
c Define the eigenvalues
c
      wc(1)=cdsqrt(dcmplx(fac1,0.d0))
      wc(2)=cdsqrt(dcmplx(fac2,0.d0))
      wc(3)=wc(2)
      wc(4)=-wc(1)
      wc(5)=-wc(2)
      wc(6)=-wc(3)
c
c     write (19,'(12e10.3)') wc
c
c Evaluation of the Eigenvectors
c Solution of the Eigenvectors
c
      do 10 i=0,3,3
c
       zc(1,1+i)=p(1)
       zc(2,1+i)=p(2)
       zc(3,1+i)=wc(1+i)
       zc(4,1+i)=2.d0*p(1)*c44*wc(1+i)
       zc(5,1+i)=2.d0*p(2)*c44*wc(1+i)
       zc(6,1+i)=den-2.d0*p2*c44
c
c      zc(4,1+i)=-2.d0*p(1)*c44*wc(1+i)
c      zc(5,1+i)=-2.d0*p(2)*c44*wc(1+i)
c      zc(6,1+i)=2.d0*p2*c44-den
c
       zc(1,2+i)=p(1)
       zc(2,2+i)=p(2)
       zc(3,2+i)=-p2/wc(2+i)
       zc(4,2+i)=p(1)*zc(6,1+i)/wc(2+i)
       zc(5,2+i)=p(2)*zc(6,1+i)/wc(2+i)
       zc(6,2+i)=-2.d0*p2*c44
c
c      zc(4,2+i)=p(1)*(2.d0*p2*c44-den)/wc(2+i)
c      zc(5,2+i)=p(2)*(2.d0*p2*c44-den)/wc(2+i)
c      zc(6,2+i)=2.d0*p2*c44
c      if (dble(wc(2+i)) .gt. 0.d0)then
c        do j=1,6
c          zc(j,2+i) = -zc(j,2+i)
c        enddo
c      endif
c
       zc(1,3+i)=-p(2)
       zc(2,3+i)=p(1)
       zc(3,3+i)=0.d0
       zc(4,3+i)=-p(2)*c44*wc(3+i)
       zc(5,3+i)=p(1)*c44*wc(3+i)
       zc(6,3+i)=0.d0
c
 10   continue
c
cs      do 40 j=1,6
cs        do 30 i=4,6
cs          zc(i,j)=zc(i,j)*(0.d0,-1.d0)
cs 30    continue
cs 40   continue
c
c Normalize the vectors based on displacement
c
      do 50 j=1,6
        acomp=zc(1,j)*zc(1,j)
        bcomp=zc(2,j)*zc(2,j)
        ccomp=zc(3,j)*zc(3,j)
        norm=cdsqrt(acomp+bcomp+ccomp)
        do 60 i=1,6
          zc(i,j)=zc(i,j)/norm
 60     continue
 50   continue
c
      return
      end
c
c........|.........|.........|.........|.........|.........|.........|..
c
