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 isotroc(a,rho,px,py,eval,evec)

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

c Finds e-vals and e-vecs for an isotropic layer.
c Based on modified (corrected) version of Sean Guest's routine isoeig.
c This version allows complex slowness and frequency.
c
      implicit real*8 (a-h,o-z)
      dimension a(3,3,3,3)
      complex*16 px,py
      complex*16 p(3),eval(6),evec(6,6)

      p(1)=px
      p(2)=py
      c33=a(3,3,3,3)*rho
      c44=a(2,3,2,3)*rho
      call isoeigc (rho,p,c33,c44,eval,evec)

      return
      end
c
c........|.........|.........|.........|.........|.........|.........|..
c
      subroutine isoeigc(den,p,c33,c44,wc,zc)

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

c
c Original version of isoeig: 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 Modified for complex lateral slowness, February 1995. CJT.
c Modified to test branches better, November 1996. CJT.
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
      complex*16 p2,p(3),wc(6),zc(6,6),norm,acomp,ccomp,bcomp
c
c-LK Correction for the normal incidence
      if(p(1).eq.0..and.p(2).eq.0.) then
        p(1)=1.E-20
      end if
c-LK
      p2=p(1)**2 +p(2)**2
c
c branch test parameter
c
      acc=1.d-12
      acc=4.d0*datan(1.d0)
c
c Define the eigenvalues
c
      wc(1)=cdsqrt(den/c33-p2)
c     test=datan2(dimag(wc(1)),dble(wc(1)))
c     if(test.lt.0.d0.and.test.ge.-acc)wc(1)=-wc(1)
c     if(test.eq.acc)wc(1)=-wc(1)

      if(dimag(wc(1)).lt.0.d0)wc(1)=-wc(1)

c       if(dble(wc(1)).eq.0.d0)then
c          if(dimag(wc(1)).eq.0.d0)then
c             print *,'*** error 1 in isotroc'
c             stop
c          endif
c          test=2.d0*acc
c       else
c          if(dimag(wc(1)).ne.0.d0)then
c             test=dabs(dimag(wc(1))/dble(wc(1)))
c          endif
c       endif
c     if(test.gt.acc)then
c         if(dimag(wc(1)).lt.0.d0)wc(1)=-wc(1)
c     else
c         if(dble(wc(1)).lt.0.d0)wc(1)=-wc(1)
c     endif

      wc(2)=cdsqrt(den/c44-p2)
c     test=datan2(dimag(wc(2)),dble(wc(2)))
c     if(test.lt.0.d0.and.test.ge.-acc)wc(2)=-wc(2)
c     if(test.eq.acc)wc(2)=-wc(2)

      if(dble(wc(2)).lt.0.d0)wc(2)=-wc(2)

c       if(dble(wc(2)).eq.0.d0)then
c          if(dimag(wc(2)).eq.0.d0)then
c             print *,'*** error 1 in isotroc'
c             stop
c          endif
c          test=2.d0*acc
c       else
c          if(dimag(wc(2)).ne.0.d0)then
c             test=dabs(dimag(wc(2))/dble(wc(2)))
c          endif
c       endif
c     if(test.gt.acc)then
c         if(dimag(wc(2)).lt.0.d0)wc(2)=-wc(2)
c     else
c         if(dble(wc(2)).lt.0.d0)wc(2)=-wc(2)
c     endif

      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)
        acomp=zc(1,j)*dconjg(zc(1,j))
        bcomp=zc(2,j)*dconjg(zc(2,j))
        ccomp=zc(3,j)*dconjg(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