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........|.........|.........|.........|.........|.........|.........|..
c
      subroutine aniso3c(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 Exploits Georg R"umpker's routines for finding vertical slowness
c given real horizontal slownesses.
c Requires Georg's routines in slsubx.f.
c
      dimension a(3,3,3,3)
      complex*16 px,py
      complex*16 eval(6),evec(6,6),pvec(3),pol(3),sigma(3)
c
c vertsl2 is ct's modified version of gr's vertsl.
      pxr=dble(px)
      pyr=dble(py)
      call vertsl2(a,pxr,pyr,eval)
c
      do 100 i=1,6
       pvec(1)=px
       pvec(2)=py
       pvec(3)=eval(i)
       call polarc(a,pvec,pol,i)


c
c Construct stress terms
c
       call stress(a,rho,pvec,pol,sigma)
c
c
       do 200 j=1,3
         evec(j,i)=pol(j)
         evec(j+3,i)=sigma(j)
  200  continue
c
  100 continue
c
      return
      end
c
c........|.........|.........|.........|.........|.........|.........|..
c
      subroutine stress(c,rho,p,pol,sigma)
      implicit real*8 (a-h,o-z)

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

c
c Compute stress vector for given displacement vector.
c
      dimension c(3,3,3,3)
      dimension cm(3,3),c33i(3,3),t(3,3),tt(3,3)
      complex*16 p(3),pol(3),sigma(3),temp
c
c Calculate the c33 inverse matrix -- first the minors
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
           c33i(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)*c33i(k,n)
            c2 = c2 + c(m,2,k,3)*c33i(k,n)
 70       continue
          tt(n,m) = -1.d0*(dreal(p(1))*c1 + dreal(p(2))*c2)
          t(m,n) = tt(n,m)
 65     continue
 60   continue
c
c Lastly the stress
c
      do 99 i=1,3
        sigma(i)=dcmplx(0.d0,0.d0)
        do 100 j=1,3
          do 101 k=1,3
            temp=dcmplx(-tt(j,k),0.d0)
            if(j.eq.k)temp=p(3)+temp
            sigma(i) = sigma(i)+
     .          c(i,3,j,3)*temp*pol(k)
  101     continue
  100   continue
c Scale by rho since a(i,j,k,l) actually passed in, not c(i,j,k,l)
        sigma(i)=rho*sigma(i)
   99 continue
c
      return
      end
c
c........|.........|.........|.........|.........|.........|.........|..
c
      subroutine vertsl2(a,p1,p2,p3)
      implicit real*8 (a-h,o-z)

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

c
c Finds vertical slowness for given horizontal slowness
c
c This version sorts real and imaginary parts into positve and
c negative groups. Does not sort within a group (unecessary?).
c CJT's version of GR's routine vertsl.
c
      dimension a(3,3,3,3),p3r(6),p3i(6)
      dimension cotr(3,3),pvec(3),vg(3,6)
      complex*16 p3(6)
c
      call eigp3(a,p1,p2,p3r,p3i)
c
c Sort according to size of p3
c
        ipos=0
        ineg=0
        acc=1.d-12
        do 100 l=1,6
c
          test=0.d0
          if(p3r(l).eq.0.d0)then
                 if(p3i(l).eq.0.d0)then
                     print *,'*** error 1 in vertsl2'
                     stop
                 endif
                 test=2.d0*acc
          else
              if(p3i(l).ne.0.d0)then
                  test=dabs(p3i(l)/p3r(l))
              endif
          endif

          if(test.gt.acc)then
              if(p3i(l).gt.0.d0)then
                 ipos=ipos+1
                 p3(ipos)=dcmplx(p3r(l),p3i(l))
              else
                 ineg=ineg+1
                 p3(ineg+3)=dcmplx(p3r(l),p3i(l))
              endif
          else
              if(p3r(l).gt.0.d0)then
                 ipos=ipos+1
                 p3(ipos)=dcmplx(p3r(l),p3i(l))
              else
                 ineg=ineg+1
                 p3(ineg+3)=dcmplx(p3r(l),p3i(l))
              endif
          endif

100     continue

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

      return
      end
c
c........|.........|.........|.........|.........|.........|.........|..
c
        subroutine polarc(a,p,pol,ivec)

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

c
c Generates polarisation vectors of
c
c   a_jk=(a_ijkl p_i p_l - delta_jk) pol_k = 0
c
c   input: slowness vector and elasticity
c   output: polarisation vectors
c
c Original version of 'polar' by G. R"umpker ---
c significanlty modified for complex values and more general
c orientations by C. Thomson. Sept. 1993
c
        implicit real*8(a-h,o-z)
        dimension a(3,3,3,3)
        complex*16 p(3),g(3,3),pol(3),de11,de22,de33,xnorm
c
c****forming gamma****
c   generating g_jk
c
         do 65 j=1,3
           do 55 k=1,3
             g(j,k) = (0.0d0,0.d0)
             do 45 i=1,3
               do 35 l=1,3
                 g(j,k) = g(j,k) + p(l)*p(i)*a(i,j,k,l)
 35            continue
 45          continue
 55        continue
 65      continue
c
c
        g(1,1)=g(1,1)-(1.d0,0.d0)
        g(2,2)=g(2,2)-(1.d0,0.d0)
        g(3,3)=g(3,3)-(1.d0,0.d0)

c        write(11,*)ivec
c        write(11,*)((g(i,j),i=1,3),j=1,3)
c        write(11,*)

         de11=g(2,2)*g(3,3)-g(2,3)*g(3,2)
         de22=g(1,1)*g(3,3)-g(1,3)*g(3,1)
         de33=g(1,1)*g(2,2)-g(1,2)*g(2,1)

c        write(11,*)ivec,de11,de22,de33
c        write(11,*)

       testv=1.d-12
       if(cdabs(de11).gt.cdabs(de22).and.cdabs(de11).gt.cdabs(de33))then
           test=cdabs(g(2,2)*g(3,3))+cdabs(g(2,3)*g(3,2))
           if(test.eq.0.d0)print *,'*** error 1 in polarc'
           if(cdabs(de11)/test.lt.testv)
     .       write(11,*)'test failed in polar ... de11=',de11
           pol(1)=de11
           pol(2)=-g(2,1)*g(3,3)+g(2,3)*g(3,1)
           pol(3)=g(2,1)*g(3,2)-g(2,2)*g(3,1)
       elseif(cdabs(de22).gt.cdabs(de33))then
           test=cdabs(g(1,1)*g(3,3))+cdabs(g(1,3)*g(3,1))
           if(test.eq.0.d0)print *,'*** error 2 in polarc'
           if(cdabs(de22)/test.lt.testv)
     .       write(11,*)'test failed in polar ... de22=',de22
           pol(1)=-g(1,2)*g(3,3)+g(1,3)*g(3,2)
           pol(2)=de22
           pol(3)=-g(1,1)*g(3,2)+g(1,2)*g(3,1)
       else
           test=cdabs(g(1,1)*g(2,2))+cdabs(g(1,2)*g(2,1))
           if(test.eq.0.d0)print *,'*** error 3 in polarc'
           if(cdabs(de33)/test.lt.testv)
     .       write(11,*)'test failed in polar ... de33=',de33
           pol(1)=g(1,2)*g(2,3)-g(1,3)*g(2,2)
           pol(2)=-g(1,1)*g(2,3)+g(1,3)*g(2,1)
           pol(3)=de33
       endif
c       xnorm=cdsqrt(pol(1)**2+pol(2)**2+pol(3)**2)
c       if(xnorm.eq.(0.d0,0.d0))print *,'*** error 4 in polarc'
        rnorm=
     .    dsqrt(cdabs(pol(1))**2+cdabs(pol(2))**2+cdabs(pol(3))**2)
        pol(1)=pol(1)/rnorm
        pol(2)=pol(2)/rnorm
        pol(3)=pol(3)/rnorm

       iwrite=0
       if(iwrite.eq.1)then
        write(11,*)
        write(11,*)'in polarc...Fedorov polarisation pol(',ivec,') = '
        write(11,*)pol
        write(11,*)
       endif

        return
        end
c
c........|.........|.........|.........|.........|.........|.........|..
c
      subroutine eigp3(a,p1,p2,p3r,p3i)

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

      implicit real*8 (a-h,o-z)
      dimension a(3,3,3,3)
      dimension alpha(3,3), beta(3,3), gamma(3,3)
      dimension axx(6),bxx(6),cxx(6),dxx(6),exx(6),fxx(6),gxx(6)
      dimension coe(7),p3r(6),p3i(6),xcoe(7)
c
c  input: matrix of elastic constants, devided by density a
c         two slowness components p1 and p2
c  output: six solutions for the slowness component p3
c          splittet into real and imaginary part
c  experience shows that
c  as soon as the eigenvalues are not elements of the
c  slowness surface the imaginary part becomes unequal 0
c

      do 1 j=1,3
       do 2 k=1,3
        alpha(j,k)=a(1,j,k,1)*p1*p1 + a(2,j,k,1)*p1*p2
     &            +a(1,j,k,2)*p2*p1 + a(2,j,k,2)*p2*p2
        if(j.eq.k) alpha(j,k)=alpha(j,k)-1.d0
        beta(j,k)=a(3,j,k,1)*p1 + a(3,j,k,2)*p2
     &           +a(1,j,k,3)*p1 + a(2,j,k,3)*p2
        gamma(j,k)=a(3,j,k,3)
2      continue
1     continue

      do 5 i=1,6
       if(i.eq.1) then
        j1=2
        k1=2
        j2=3
        k2=3
        j3=1
        k3=1
       else if(i.eq.2) then
        j1=2
        k1=3
        j2=3
        k2=2
        j3=1
        k3=1
       else if(i.eq.3) then
        j1=2
        k1=3
        j2=3
        k2=1
        j3=1
        k3=2
       else if(i.eq.4) then
        j1=2
        k1=1
        j2=3
        k2=3
        j3=1
        k3=2
       else if(i.eq.5) then
        j1=2
        k1=1
        j2=3
        k2=2
        j3=1
        k3=3
       else
        j1=2
        k1=2
        j2=3
        k2=1
        j3=1
        k3=3
       end if

       ax=alpha(j1,k1)*alpha(j2,k2)
       bx=alpha(j1,k1)*beta(j2,k2)+beta(j1,k1)*alpha(j2,k2)
       cx=alpha(j1,k1)*gamma(j2,k2)
     &    +beta(j1,k1)*beta(j2,k2)
     &    +gamma(j1,k1)*alpha(j2,k2)
       dx=beta(j1,k1)*gamma(j2,k2)+gamma(j1,k1)*beta(j2,k2)
       ex=gamma(j1,k1)*gamma(j2,k2)

       axx(i)=alpha(j3,k3)*ax
       bxx(i)=alpha(j3,k3)*bx+beta(j3,k3)*ax
       cxx(i)=alpha(j3,k3)*cx+beta(j3,k3)*bx+gamma(j3,k3)*ax
       dxx(i)=alpha(j3,k3)*dx+beta(j3,k3)*cx+gamma(j3,k3)*bx
       exx(i)=alpha(j3,k3)*ex+beta(j3,k3)*dx+gamma(j3,k3)*cx
       fxx(i)=beta(j3,k3)*ex+gamma(j3,k3)*dx
       gxx(i)=gamma(j3,k3)*ex

5     continue

      coe(1)=axx(1)-axx(2)+axx(3)-axx(4)+axx(5)-axx(6)
      coe(2)=bxx(1)-bxx(2)+bxx(3)-bxx(4)+bxx(5)-bxx(6)
      coe(3)=cxx(1)-cxx(2)+cxx(3)-cxx(4)+cxx(5)-cxx(6)
      coe(4)=dxx(1)-dxx(2)+dxx(3)-dxx(4)+dxx(5)-dxx(6)
      coe(5)=exx(1)-exx(2)+exx(3)-exx(4)+exx(5)-exx(6)
      coe(6)=fxx(1)-fxx(2)+fxx(3)-fxx(4)+fxx(5)-fxx(6)
      coe(7)=gxx(1)-gxx(2)+gxx(3)-gxx(4)+gxx(5)-gxx(6)

      coe(1)=coe(1)/coe(7)
      coe(2)=coe(2)/coe(7)
      coe(3)=coe(3)/coe(7)
      coe(4)=coe(4)/coe(7)
      coe(5)=coe(5)/coe(7)
      coe(6)=coe(6)/coe(7)
      coe(7)=1.d0

      call dpoly(coe,xcoe,6,p3r,p3i,ier)

      return
      end
c
c........|.........|.........|.........|.........|.........|.........|..
c
      subroutine dpoly(xcof,cof,m,rootr,rooti,ier)
c=======================================================================
c  Extract roots of a polynomial/double precision
c                                       Last revision: February 29, 1988
c  xcof: coeffizienten des polynoms, maximal 37
c         dcoef(1) ist das konstante Glied
c  cof:hilfsvector, der nicht benoetigt wird, aber deklarieren
c  m:   grad des polynoms, maximal 36
c  rootr: vector mit realteilen der wurzeln
c  rooti: vector mit imaginaerteilen der Wurzeln
c  ier:   fehlerangabe
c=======================================================================
      dimension cof(37),rooti(36),rootr(36),xcof(37)
      double precision cof,rooti,rootr,xcof
      double precision alpha,dabs,dx,dy,sumsq,temp,u,ux,uy
      double precision v,x,xo,xpr,xt,xt2,y,yo,ypr,yt,yt2
c=======================================================================
c  Check for input data errors
c=============================
      ifit=0
      n=m
      ier=0
      if(xcof(n+1).eq.0.) then
        ier=4

      elseif(n.le.0) then
        ier=1
      elseif(n.gt.36) then
        ier=2
      else
c================================
c  No input data errors, so begin
c================================
        nx=n
        nxx=n+1
        n2=1
        kj1=n+1
        do 10 l=1,kj1
          mt=kj1-l+1
          cof(mt)=xcof(l)
 10     continue
c====================
c  Set initial values
c====================
c 20     xo=.00500101d0    feb.92 sometimes wrong ?
 20     xo=-.00500101d0
        yo=.01000101d0
c============================
c  Zero initial value counter
c============================
        in=0
 30     x=xo
c======================================
c  Increment initial values and counter
c======================================
        xo=-10.d0*yo
        yo=-10.d0*x
c==============================
c  Set x and y to current value
c==============================
        x=xo
        y=yo
        in=in+1
        go to 50
 40     ifit=1
        xpr=x
        ypr=y
c=====================================
c  Evaluate polynomial and derivatives
c=====================================
 50     ict=0
 60     ux=0.d0
        uy=0.d0
        v =0.d0
        yt=0.d0
        xt=1.d0
        u=cof(n+1)
        if(u.eq.0.d0) go to 120
        do 70 i=1,n
          l=n-i+1
          temp=cof(l)
          xt2=x*xt-y*yt
          yt2=x*yt+y*xt
          u=u+temp*xt2
          v=v+temp*yt2
          fi=i
          ux=ux+fi*xt*temp
          uy=uy-fi*yt*temp
          xt=xt2
          yt=yt2
 70     continue
        sumsq=ux*ux+uy*uy
        if(sumsq.eq.0.d0) go to 100
        dx=(v*uy-u*ux)/sumsq
        x=x+dx
        dy=-(u*uy+v*ux)/sumsq
        y=y+dy
        if((dabs(dy)+dabs(dx)).lt.1.d-12) go to 80
c========================
c  Step iteration counter
c========================
        ict=ict+1
        if(ict.lt.500) go to 60
        if(ifit.ne.0) go to 80
        if(in.eq.5) then
          ier=3
          go to 170
        else
          go to 30
        endif
 80     do 90 l=1,nxx
          mt=kj1-l+1
          temp=xcof(mt)
          xcof(mt)=cof(l)
          cof(l)=temp
 90     continue
        itemp=n
        n=nx
        nx=itemp
        if(ifit.eq.0) then
          go to 40
        else
          go to 110
        endif
 100    if(ifit.eq.0) go to 30
        x=xpr
        y=ypr
 110    ifit=0
        if(dabs(y).lt.1.d-10*dabs(x)) go to 130
        alpha=x+x
        sumsq=x*x+y*y
        n=n-2
        go to 140
 120    x=0.d0
        nx=nx-1
        nxx=nxx-1
 130    y=0.d0
        sumsq=0.d0
        alpha=x
        n=n-1
 140    cof(2)=cof(2)+alpha*cof(1)
        do 150 l=2,n
          cof(l+1)=cof(l+1)+alpha*cof(l)-sumsq*cof(l-1)
 150    continue
 160    rooti(n2)=y
        rootr(n2)=x
        n2=n2+1
        if(sumsq.ne.0.d0) then
          y=-y
          sumsq=0.d0
          go to 160
        endif
        if(n.gt.0) go to 20
      endif
 170  return
      end
c
c........|.........|.........|.........|.........|.........|.........|..
c