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
      subroutine aniso6cg(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
c Finds e-vals and e-vecs for given layer.
c Based on EISPACK routines for finding vertical slowness given
c horizontal slownesses by solving the 6X6 matrix problem.
c
c COMPLEX SYSTEM MATRIX ALLOWED -- this version is for complex
c horizontal slowness.
c
      dimension a(3,3,3,3)
      complex*16 px,py,p(3),eval(6),evec(6,6),xnorm
      logical iso
      common/info/lay,id,iso

      iso=.false.

      p(1)=px
      p(2)=py
      call eigc (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 eigc (den,ae,p,wc,zc)

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

c
c          Based on subroutine eig by:
c          W.S. Guest,  October 1988.
c
c          Modified by C J Thomson, February 1995.
c          This version allows for complex horizontal slowness.
c
c Sean's comments from the original follow (but do not relate exactly
c to the new version) ...
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
c This new version calls isoeigc or Eispack routine cg and the sorting
c is done by routine sortslwc.
c
c C. Thomson February 1995.
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)
      dimension ct2(3,3,3),ct(3,3),qc(3,3,3,3)
      dimension amr(6,6),ami(6,6)
      dimension fv1(6),fv2(6),fv3(6)
      dimension zi(6,6),wr(6),wi(6),zr(6,6)
      complex*16 p(3),wc(6),zc(6,6)
      complex*16 t(3,3),tt(3,3),s(3,3)
c Common Block
      common/info/lay,id,iso
c
c
c     write (*,*) 'in eigc',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 (*,*)'isoeigc used to solve "A"'
         c33=c(3,3,3,3)
         c44=c(2,3,2,3)
c        write (*,*) c33,c44
         call isoeigc(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,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
          amr(i,k) =  dreal(tt(i,k))
          ami(i,k) =  dimag(tt(i,k))
          amr(i,l) =  ci(i,k)
          ami(i,l) =  0.d0
          amr(j,k) =  dreal(s(i,k))
          ami(j,k) =  dimag(s(i,k))
          amr(j,l) =  dreal(t(i,k))
          ami(j,l) =  dimag(t(i,k))
 140    continue
 135  continue
c
c Find the eigenvalues and eigenvectors
c
      n6=6
      call cg(n6,n6,amr,ami,wr,wi,1,zr,zi,fv1,fv2,fv3,ierr)

      if(ierr.ne.0)print *,'**IERR.ne.0 afer CG -- in eigc'
c
      call sortslwc(wr,wi,zr,zi,wc,zc)
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 sortslwc(wr,wi,zr,zi,wc,zc)

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

      implicit real*8 (a-h,o-z)
c
c Sorts eigenvalues/vectors into separate up/down sets.
c
c This version sorts real and imaginary parts into positive and
c negative groups. Does not sort within a group (unecessary?).
c
c C. Thomson, February 1995.
c
      dimension wr(6),wi(6),zr(6,6),zi(6,6)
      dimension wwrkr(6),wwrki(6),zwrkr(6,6),zwrki(6,6)
      complex*16 wc(6),zc(6,6)
c
c Sort according to size of p3 (vertical slowness)
c
        ipos=0
        ineg=0
        acc=1.d-12
        do 100 l=1,6
c
          test=0.d0
          if(wr(l).eq.0.d0)then
                 if(wi(l).eq.0.d0)then
                     print *,'*** error 1 in sortslwc'
                     stop
                 endif
                 test=2.d0*acc
          else
              if(wi(l).ne.0.d0)then
                  test=dabs(wi(l)/wr(l))
              endif
          endif

          if(test.gt.acc)then
              if(wi(l).gt.0.d0)then
                 ipos=ipos+1
                 wwrkr(ipos)=wr(l)
                 wwrki(ipos)=wi(l)
                 call trnsfv6(zr(1,l),zi(1,l),
     .                    zwrkr(1,ipos),zwrki(1,ipos))
              else
                 ineg=ineg+1
                 wwrkr(ineg+3)=wr(l)
                 wwrki(ineg+3)=wi(l)
                 call trnsfv6(zr(1,l),zi(1,l),
     .                    zwrkr(1,ineg+3),zwrki(1,ineg+3))
              endif
          else
              if(wr(l).gt.0.d0)then
                 ipos=ipos+1
                 wwrkr(ipos)=wr(l)
                 wwrki(ipos)=wi(l)
                 call trnsfv6(zr(1,l),zi(1,l),
     .                    zwrkr(1,ipos),zwrki(1,ipos))
              else
                 ineg=ineg+1
                 wwrkr(ineg+3)=wr(l)
                 wwrki(ineg+3)=wi(l)
                 call trnsfv6(zr(1,l),zi(1,l),
     .                    zwrkr(1,ineg+3),zwrki(1,ineg+3))
              endif
          endif

100     continue

         if(ipos.ne.3.or.ineg.ne.3)print *,'**error 2 in sortslwc'

         do 101 i=1,6
               wc(i)=dcmplx(wwrkr(i),wwrki(i))
         do 102 l=1,6
               zc(l,i)=dcmplx(zwrkr(l,i),zwrki(l,i))
  102    continue
  101    continue

      return
      end
c
c........|.........|.........|.........|.........|.........|.........|..
c
       subroutine trnsfv6(vin1,vin2,vout1,vout2)

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

       implicit real*8 (a-h,o-z)

       dimension vin1(6),vin2(6),vout1(6),vout2(6)

       do 100 l=1,6
          vout1(l)=vin1(l)
          vout2(l)=vin2(l)
100    continue

       return
       end
c
c........|.........|.........|.........|.........|.........|.........|..
c
