C
```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
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
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.

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.

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
C                                                                 ```