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

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