C
```c
c........|.........|.........|.........|.........|.........|.........|..
c
c Subroutines used by program rmatrix1.f, which
c computes r/t matrices for a stack of homogeneous
c layers. Uses matrix algorithm of B. L. N. Kennett.
c General anisotropy permitted.
c
c C. J. Thomson.
c
c This version allows complex px,py and omega. February 1995.
c
c See subroutine evset (next) for important comments regarding
c eigenvalue/eigenvector calculations.
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 evset(px,py,omega,ilay,eval,evec)

c Copyright (c) 1996 C. J. Thomson.

c Sets up the e-values and e-vectors for a given layer.
c
c Calls routine isotroc f the layer is isotropic.
c
c For anisotropic layers variable itype=1,2,3 controls the method
c used to find the e-vals/e-vecs.
c
c itype=1 : calls Eispack routine cg for a complex general system
c           matrix A -- i.e. for the case of complex slowness or
c           elastic parameters (note complex a(i,j,k,l) not allowed
c           at present but code easily extended).
c
c itype=2 : for the case of a real system matrix A; calls aniso6c1
c           which in turn calls Sean Guest's routine eig which uses
c           public domain e-value/e-vector routines that look
c           similar to Eispack routines.
c
c itype=3 : again for real A matrices; this time the vertical
c           slownesses are found by a 3x3 matrix method rather
c           like ray theory. Calls aniso3c which uses Georg
c           R"umpker's routine eigp3 which uses a public-domain
c           routine for the roots of a sextic.
c
c See routines aniso6cg, aniso6c1 and polarc to check  normalisation
c of the eigenvectors (amplitude of displacement part set to
c unity).
c
c IMPORTANT NOTE CONCERNING EIGENVECTOR ORDERING:
c The first three e-vals/e-vecs correspond to downgoing or
c downwardly decaying waves and the last three to upgoing or
c upwardly decaying waves. This is all that the R/T matrix
c calculation really requires. In the isotropic case there is
c easy sorting within these two groups into P, SV and SH.
c In the anisotropic case the sorting is generally from largest
c to smallest vertical slowness. Polarisations in anisotropic
c media can behave in interesting ways and so it is assumed the
c user will take care to understand which wave function should
c be called, say, qSH and which column it occupies (simply
c by looking at the eigenvectors).
c

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

complex*16 px,py,omega

complex*16 eval(6), evec(6,6)
c
c COMMON BLOCKS to be INCLUDEd
c
c max number of layers ---    1 = half-space above or top layer
c    in the stack                 just below free surface
c                          nlay = half-space below
c
parameter (nlaymx=30)
c
c Earth model
c
common /elast/isoflg(nlaymx),thickn(nlaymx),a(3,3,3,3,nlaymx),
.              rho(nlaymx)
c

complex*16 evals,evecs
common /evprob/evals(6,nlaymx),evecs(6,6,nlaymx)

complex*16 q11i,q12i,q21i,q22i
complex*16 tui,rui,tdi,rdi
common /intfce/q11i(3,3,nlaymx),q12i(3,3,nlaymx),
.               q21i(3,3,nlaymx),q22i(3,3,nlaymx),
.               tui(3,3,nlaymx),rui(3,3,nlaymx),
.               tdi(3,3,nlaymx),rdi(3,3,nlaymx)

complex*16 tusi,rusi,tdsi,rdsi
common /stackrt/tusi(3,3,nlaymx),rusi(3,3,nlaymx),
.                tdsi(3,3,nlaymx),rdsi(3,3,nlaymx)

c
c symmetry testing data
c
complex*16 diag1,j2cst1,j2cstn
common /symts/diag1(6,nlaymx)
common /symts2/j2cst1(3,3),j2cstn(3,3)
c
c END OF COMMON BLOCKS
c
if(isoflg(ilay).eq.1)then
call isotroc(a(1,1,1,1,ilay),rho(ilay),
.                     px,py,eval,evec)
c-LK       call write1(ilay,isoflg(ilay),px,py,eval,evec)
else
itype=1
if(itype.eq.1)then
c find e-vals + e-vecs of complex general matrix ...
call aniso6cg(a(1,1,1,1,ilay),rho(ilay),
.                      px,py,eval,evec)
c-LK         call write1(ilay,isoflg(ilay),px,py,eval,evec)
elseif(itype.eq.2)then
c find e-vals + e-vecs for case of real horizontal slownesses ...
call aniso6c1(a(1,1,1,1,ilay),rho(ilay),
.                      px,py,eval,evec)
c-LK         call write1(ilay,isoflg(ilay),px,py,eval,evec)
elseif(itype.eq.3)then
c find e-vals + e-vecs for case of real horizontal slownesses
c using 3x3 matrix method (like ray theory) ...
call aniso3c(a(1,1,1,1,ilay),rho(ilay),
.                     px,py,eval,evec)
c-LK         call write1(ilay,isoflg(ilay),px,py,eval,evec)
endif
endif

return
end
c
c........|.........|.........|.........|.........|.........|.........|..
c
subroutine estore(eval,evec,evals,evecs)

c Copyright (c) 1996 C. J. Thomson.

complex*16 eval(6),evec(6,6),evals(6),evecs(6,6)

do 100 i=1,6

evals(i)=eval(i)

do 200 j=1,6

evecs(i,j)=evec(i,j)

200       continue

100   continue

return
end
c
c........|.........|.........|.........|.........|.........|.........|..
c
subroutine qstore(q11,q12,q21,q22,q11i,q12i,q21i,q22i)

c Copyright (c) 1996 C. J. Thomson.

complex*16 q11(3,3),q12(3,3),q21(3,3),q22(3,3)
complex*16 q11i(3,3),q12i(3,3),q21i(3,3),q22i(3,3)

do 100 i=1,3

do 200 j=1,3

q11i(i,j)=q11(i,j)
q12i(i,j)=q12(i,j)
q21i(i,j)=q21(i,j)
q22i(i,j)=q22(i,j)

200       continue

100   continue

return
end
c
c........|.........|.........|.........|.........|.........|.........|..
c
subroutine rtstore(tu,ru,td,rd,tui,rui,tdi,rdi)

c Copyright (c) 1996 C. J. Thomson.

complex*16 tu(3,3),ru(3,3),td(3,3),rd(3,3)
complex*16 tui(3,3),rui(3,3),tdi(3,3),rdi(3,3)

do 100 i=1,3

do 200 j=1,3

tui(i,j)=tu(i,j)
rui(i,j)=ru(i,j)
rdi(i,j)=rd(i,j)
tdi(i,j)=td(i,j)

200       continue

100   continue

return
end
c
c........|.........|.........|.........|.........|.........|.........|..
c
subroutine rtget(tu,ru,td,rd,tui,rui,tdi,rdi)

c Copyright (c) 1996 C. J. Thomson.

complex*16 tu(3,3),ru(3,3),td(3,3),rd(3,3)
complex*16 tui(3,3),rui(3,3),tdi(3,3),rdi(3,3)

do 100 i=1,3

do 200 j=1,3

tu(i,j)=tui(i,j)
ru(i,j)=rui(i,j)
rd(i,j)=rdi(i,j)
td(i,j)=tdi(i,j)

200       continue

100   continue

return
end
c
c........|.........|.........|.........|.........|.........|.........|..
c
subroutine xmult6(a,b,c)

c Copyright (c) 1996 C. J. Thomson.

c
c complex matrix multiplication a*b=c
c dimension 6 ... allows overwriting
c
complex*16 a(6,6), b(6,6), c(6,6), wrk(6,6)

do 11 i=1,6
do 13 j=1,6
wrk(i,j)=(0.d0,0.d0)
do 15 k=1,6
wrk(i,j)=wrk(i,j) + a(i,k)*b(k,j)
15        continue
13      continue
11    continue

do 17 i=1,6
do 19 j=1,6
c(i,j)=wrk(i,j)
19      continue
17    continue

return
end
c
c........|.........|.........|.........|.........|.........|.........|..
c
subroutine xeveci(nlay,ilay,evec,eveci)

c Copyright (c) 1996 C. J. Thomson.

implicit real*8(a-h,o-z)
c
c Inverts the eigenvector matrix, and sets up the diagonal 'constant'
c of J-primed symmetry in the process (see theory notes, section 6,
c equations (6.5,6.6)).
c
complex*16 evec(6,6), eveci(6,6), wrk(6,6)

c
c COMMON BLOCKS to be INCLUDEd
c
c max number of layers ---    1 = half-space above or top layer
c    in the stack                 just below free surface
c                          nlay = half-space below
c
parameter (nlaymx=30)
c
c Earth model
c
common /elast/isoflg(nlaymx),thickn(nlaymx),a(3,3,3,3,nlaymx),
.              rho(nlaymx)

complex*16 evals,evecs
common /evprob/evals(6,nlaymx),evecs(6,6,nlaymx)

complex*16 q11i,q12i,q21i,q22i
complex*16 tui,rui,tdi,rdi
common /intfce/q11i(3,3,nlaymx),q12i(3,3,nlaymx),
.               q21i(3,3,nlaymx),q22i(3,3,nlaymx),
.               tui(3,3,nlaymx),rui(3,3,nlaymx),
.               tdi(3,3,nlaymx),rdi(3,3,nlaymx)

complex*16 tusi,rusi,tdsi,rdsi
common /stackrt/tusi(3,3,nlaymx),rusi(3,3,nlaymx),
.                tdsi(3,3,nlaymx),rdsi(3,3,nlaymx)

c
c symmetry testing data
c
complex*16 diag1,j2cst1,j2cstn
common /symts/diag1(6,nlaymx)
common /symts2/j2cst1(3,3),j2cstn(3,3)
c
c END OF COMMON BLOCKS
c

do 11 i=1,3
do 15 j=1,3
eveci(i,j)=evec(j+3,i)
eveci(i,j+3)=evec(j,i)
eveci(i+3,j)=evec(j+3,i+3)
eveci(i+3,j+3)=evec(j,i+3)
15      continue
11    continue

call xmult6(eveci,evec,wrk)

iwrite=0
if(iwrite.eq.1)then
write(11,*)
write(11,*)'------------------------------------------'
write(11,*)
write(11,*)'in xeveci...diag constant, layer ',ilay,'...'
write(11,*)(wrk(i,i),i=1,6)
endif

do 17 i=1,6
diag1(i,nlay+1-ilay)=wrk(i,i)
do 19 j=1,6
eveci(i,j)=eveci(i,j)/wrk(i,i)
19      continue
17    continue

if(iwrite.eq.1)then
write(11,*)
write(11,*)'------------------------------------------'
write(11,*)
endif

return
end
c
c........|.........|.........|.........|.........|.........|.........|..
c
subroutine trnsf6(e1,e2,x1,x2)

c Copyright (c) 1996 C. J. Thomson.

complex*16 e1(6),e2(6),x1(6,6),x2(6,6)

do 100 i=1,6
e1(i)=e2(i)
do 200 j=1,6
x1(i,j)=x2(i,j)
200   continue
100 continue

return
end
c
c........|.........|.........|.........|.........|.........|.........|..
c
subroutine retra(q11,q12,q21,q22,tu,ru,td,rd)

c Copyright (c) 1996 C. J. Thomson.

c Extracts the R/T coefficients from the Q-matrix for a given single
c interface (see theory notes, equation (3.4)).

implicit real*8(a-h,o-z)
complex*16 q11(3,3),q12(3,3),q21(3,3),q22(3,3)

complex*16 tu(3,3),ru(3,3),td(3,3),rd(3,3)
complex*16 tuneg(3,3),xhelp(3,3)

call xinv3(q22,tu)

call xmult3(q12,tu,ru)

call xneg3(tu,tuneg)

call xmult3(tuneg,q21,rd)

call xmult3(q12,rd,xhelp)

return
end
c
c........|.........|.........|.........|.........|.........|.........|..
c
subroutine qmatrix(eveci1,evec2,q11,q12,q21,q22)

c Copyright (c) 1996 C. J. Thomson.

c Sets up the Q matrix for a single interface from the e-vectors
c defined on either side (see theory notes, equation (3.3)).

implicit real*8(a-h,o-z)
complex*16 eveci1(6,6),evec2(6,6),qq(6,6)
complex*16 q11(3,3),q12(3,3),q21(3,3),q22(3,3)

call xmult6(eveci1,evec2,qq)

do 17 i=1,3
do 19 j=1,3
q11(i,j)=qq(i,j)
q12(i,j)=qq(i,j+3)
q21(i,j)=qq(i+3,j)
q22(i,j)=qq(i+3,j+3)
19      continue
17    continue

return
end
c
c........|.........|.........|.........|.........|.........|.........|..
c
subroutine xmult3(a,b,c)

c Copyright (c) 1996 C. J. Thomson.

c
c Complex matrix multiplication a*b=c
c dimension 3 ... allows overwriting
c
complex*16 a(3,3), b(3,3), c(3,3), wrk(3,3)

do 11 i=1,3
do 13 j=1,3
wrk(i,j)=(0.d0,0.d0)
do 15 k=1,3
wrk(i,j)=wrk(i,j) + a(i,k)*b(k,j)
15        continue
13      continue
11    continue

do 17 i=1,3
do 19 j=1,3
c(i,j)=wrk(i,j)
19      continue
17    continue

return
end
c
c........|.........|.........|.........|.........|.........|.........|..
c
subroutine xneg3(a,aneg)

c Copyright (c) 1996 C. J. Thomson.

c
c Complex matrix multiplication with (-1)
c
complex*16 a(3,3), aneg(3,3)

do 11 i=1,3
do 13 j=1,3
aneg(i,j)=-a(i,j)
13      continue
11    continue

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

c Copyright (c) 1996 C. J. Thomson.

c
c
complex*16 a(3,3), b(3,3), c(3,3)

do 11 i=1,3
do 13 j=1,3
c(i,j)=a(i,j)+b(i,j)
13      continue
11    continue

return
end
c
c........|.........|.........|.........|.........|.........|.........|..
c
subroutine xinv3(a,ainv)

c Copyright (c) 1996 C. J. Thomson.

complex*16 a(3,3),ainv(3,3)
complex*16 co(3,3), deta, detb, detc
c
c Inversion of 3x3 matrix
c

co(1,1)=(a(2,2)*a(3,3)-a(2,3)*a(3,2))
co(1,2)=-(a(2,1)*a(3,3)-a(2,3)*a(3,1))
co(1,3)=(a(2,1)*a(3,2)-a(2,2)*a(3,1))

co(2,1)=-(a(1,2)*a(3,3)-a(1,3)*a(3,2))
co(2,2)=(a(1,1)*a(3,3)-a(1,3)*a(3,1))
co(2,3)=-(a(1,1)*a(3,2)-a(1,2)*a(3,1))

co(3,1)=(a(1,2)*a(2,3)-a(1,3)*a(2,2))
co(3,2)=-(a(1,1)*a(2,3)-a(1,3)*a(2,1))
co(3,3)=(a(1,1)*a(2,2)-a(1,2)*a(2,1))

deta=a(1,1)*co(1,1)+a(1,2)*co(1,2)+a(1,3)*co(1,3)
detb=a(2,1)*co(2,1)+a(2,2)*co(2,2)+a(2,3)*co(2,3)
detc=a(3,1)*co(3,1)+a(3,2)*co(3,2)+a(3,3)*co(3,3)

c     write(11,*)
c     write(11,*)'det check in xinv3 ... deta, detb, detc ='
c     write(11,*)deta,detb,detc

do 10 i=1,3
do 11 j=1,3
ainv(i,j)=co(j,i)/deta
11      continue
10    continue

return
end
c
c........|.........|.........|.........|.........|.........|.........|..
c
.                  tu,ru,td,rd)

c Copyright (c) 1996 C. J. Thomson.

c
c Addition/recursion formulas (see theory notes, section 4,
c equations (4.8)).
c

complex*16 tus(3,3),rus(3,3),tds(3,3),rds(3,3)
complex*16 tu(3,3),ru(3,3),td(3,3),rd(3,3)

c workspaces

complex*16 tusw(3,3),rusw(3,3),tdsw(3,3),rdsw(3,3)
complex*16 ximat(3,3),reverb(3,3),reverbi(3,3)
save ximat

if(ilay.eq.nlay-1)then
c first interface -- trivial
call trnsf3(tus,tu)
call trnsf3(rus,ru)
call trnsf3(tds,td)
call trnsf3(rds,rd)
do 100 i=1,3
do 101 j=1,3
ximat(i,j)=(0.d0,0.d0)
101       continue
ximat(i,i)=(1.d0,0.d0)
100     continue
return
endif

iwrite=0
if(iwrite.eq.1)then
write(11,*)
write(11,*)'--------------------------------------------------'
write(11,*)

do 200 i=1,3
write(11,1000)(ximat(i,j),j=1,3)
200     continue
endif
1000  format(3(' (',d10.4,',',d10.4,')'))

c symbolic notation...    reverb = ximat - rds * ru
call xmult3(rds,ru,reverb)
call xneg3(reverb,reverb)
call xinv3(reverb,reverbi)

if(iwrite.eq.1)then
write(11,*)
write(11,*)'reverb and reverbi...'
do 201 i=1,3
write(11,1000)(reverb(i,j),j=1,3)
201     continue
do 202 i=1,3
write(11,1000)(reverbi(i,j),j=1,3)
202     continue
endif

c     tusw = tu * reverbi * tus
call xmult3(reverbi,tus,tusw)
call xmult3(tu,tusw,tusw)

c     rdsw = rd + tu * reverbi * rds * td
call xmult3(rds,td,rdsw)
call xmult3(reverbi,rdsw,rdsw)
call xmult3(tu,rdsw,rdsw)

c     rusw = rus + tds * ru * reverbi * tus
call xmult3(reverbi,tus,rusw)
call xmult3(ru,rusw,rusw)
call xmult3(tds,rusw,rusw)

c     tdsw = tds * [ ximat + ru*reverbi*rds ] * td
call xmult3(reverbi,rds,tdsw)
call xmult3(ru,tdsw,tdsw)
call xmult3(tdsw,td,tdsw)
call xmult3(tds,tdsw,tdsw)

c copy over

call trnsf3(tus,tusw)
call trnsf3(rus,rusw)
call trnsf3(tds,tdsw)
call trnsf3(rds,rdsw)

return
end
c
c........|.........|.........|.........|.........|.........|.........|..
c
subroutine layer(ilay,omega,xh,eval,tu,ru,td,rd)

c Copyright (c) 1996 C. J. Thomson.

c
c Includes the phase shift across a layer (theory notes, section 4,
c equations (4.6)).
c

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

complex*16 omega
complex*16 xi
complex*16 eval(6)
complex*16 tu(3,3),ru(3,3),td(3,3),rd(3,3)
complex*16 ed(3,3),eui(3,3)

iwrite=0
if(iwrite.eq.1)then
write(11,*)
write(11,*)'.................................................'
write(11,*)
write(11,*)'crossing layer ',ilay,'   omega, xh = ',omega,xh
write(11,*)
endif

xi=(0.d0,1.d0)

do 11 i=1,3
do 13 j=1,3
ed(i,j)=(0.d0,0.d0)
eui(i,j)=(0.d0,0.d0)
13      continue
11    continue

c
c elements of layer matrix for downgoing waves
c
ed(1,1)=cdexp(xi*omega*eval(1)*xh)
if(cdabs(ed(1,1)).gt.1.d0+1.d-12)
.print *,'**error ed11 in layer |ed11|=',cdabs(ed(1,1))
ed(2,2)=cdexp(xi*omega*eval(2)*xh)
if(cdabs(ed(2,2)).gt.1.d0+1.d-12)
.print *,'**error ed22 in layer |ed22|=',cdabs(ed(2,2))
ed(3,3)=cdexp(xi*omega*eval(3)*xh)
if(cdabs(ed(3,3)).gt.1.d0+1.d-12)
.print *,'**error ed33 in layer |ed33|=',cdabs(ed(3,3))
c
c elements of inverse layer matrix for upgoing waves
c
eui(1,1)=cdexp(-xi*omega*eval(4)*xh)
if(cdabs(eui(1,1)).gt.1.d0+1.d-12)
.print *,'**error eui11 in layer |eui11|=',cdabs(eui(1,1))
eui(2,2)=cdexp(-xi*omega*eval(5)*xh)
if(cdabs(eui(2,2)).gt.1.d0+1.d-12)
.print *,'**error eui22 in layer |eui22|=',cdabs(eui(2,2))
eui(3,3)=cdexp(-xi*omega*eval(6)*xh)
if(cdabs(eui(3,3)).gt.1.d0+1.d-12)
.print *,'**error eui33 in layer |eui33|=',cdabs(eui(3,3))

call xmult3(td,ed,td)

call xmult3(rd,ed,rd)
call xmult3(eui,rd,rd)

call xmult3(eui,tu,tu)

if(iwrite.eq.1)then
write(11,*)
write(11,*)'.................................................'
write(11,*)
endif

return
end
c
c........|.........|.........|.........|.........|.........|.........|..
c
subroutine trnsf3(x1,x2)

c Copyright (c) 1996 C. J. Thomson.

complex*16 x1(3,3),x2(3,3)

do 100 i=1,3
do 200 j=1,3
x1(i,j)=x2(i,j)
200   continue
100 continue

return
end
c
c........|.........|.........|.........|.........|.........|.........|..
c
subroutine j1symtst(nlay,ilay,tus,rus,tds,rds)

c Copyright (c) 1996 C. J. Thomson.

c Performs the J_1 symmetry test. See section 6 of the
c theory notes (equation (6.21)). Symmetries here satisfied
c by perfectly elastic medium, real slownesses/frequency
c (propagating waves) and general anisotropy.

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

complex*16 tus(3,3),rus(3,3),tds(3,3),rds(3,3)
complex*16 wrk1(3,3),wrk2(3,3),wrk3(3,3),wrk4(3,3)
complex*16 wrk5(3,3),wrk6(3,3),wrk7(3,3),wrk8(3,3)

c
c COMMON BLOCKS to be INCLUDEd
c
c max number of layers ---    1 = half-space above or top layer
c    in the stack                 just below free surface
c                          nlay = half-space below
c
parameter (nlaymx=30)
c
c Earth model
c
common /elast/isoflg(nlaymx),thickn(nlaymx),a(3,3,3,3,nlaymx),
.              rho(nlaymx)

complex*16 evals,evecs
common /evprob/evals(6,nlaymx),evecs(6,6,nlaymx)

complex*16 q11i,q12i,q21i,q22i
complex*16 tui,rui,tdi,rdi
common /intfce/q11i(3,3,nlaymx),q12i(3,3,nlaymx),
.               q21i(3,3,nlaymx),q22i(3,3,nlaymx),
.               tui(3,3,nlaymx),rui(3,3,nlaymx),
.               tdi(3,3,nlaymx),rdi(3,3,nlaymx)

complex*16 tusi,rusi,tdsi,rdsi
common /stackrt/tusi(3,3,nlaymx),rusi(3,3,nlaymx),
.                tdsi(3,3,nlaymx),rdsi(3,3,nlaymx)

c
c symmetry testing data
c
complex*16 diag1,j2cst1,j2cstn
common /symts/diag1(6,nlaymx)
common /symts2/j2cst1(3,3),j2cstn(3,3)
c
c END OF COMMON BLOCKS
c

do 100 i=1,3
do 101 j=1,3

wrk1(i,j)= diag1(i,1)*rus(i,j)
wrk2(i,j)= diag1(i+3,nlay+1-ilay)*tus(i,j)
wrk3(i,j)= diag1(i,1)*tds(i,j)
wrk4(i,j)= diag1(i+3,nlay+1-ilay)*rds(i,j)

101     continue
100   continue

tst1=(0.d0,0.d0)
tst2=(0.d0,0.d0)

do 200 i=1,3
do 201 j=1,3

wrk5(i,j)=(0.d0,0.d0)
wrk6(i,j)=(0.d0,0.d0)
wrk7(i,j)=(0.d0,0.d0)
wrk8(i,j)=(0.d0,0.d0)

do 202 k=1,3

wrk5(i,j)= wrk5(i,j)+dconjg(rus(k,i))*wrk1(k,j)
.                     - dconjg(tus(k,i))*wrk2(k,j)
wrk6(i,j)= wrk6(i,j)+dconjg(rus(k,i))*wrk3(k,j)
.                     - dconjg(tus(k,i))*wrk4(k,j)
wrk7(i,j)= wrk7(i,j)+dconjg(tds(k,i))*wrk1(k,j)
.                     - dconjg(rds(k,i))*wrk2(k,j)
wrk8(i,j)= wrk8(i,j)+dconjg(tds(k,i))*wrk3(k,j)
.                     - dconjg(rds(k,i))*wrk4(k,j)

202       continue

c form norms for accuracy/confidence measure

tst1=tst1+cdabs(wrk6(i,j))+cdabs(wrk7(i,j))
if(i.ne.j)then
tst1=tst1+cdabs(wrk5(i,j))+cdabs(wrk8(i,j))
endif
if(i.eq.j)then
tst2=tst2+cdabs(wrk5(i,j))+cdabs(wrk8(i,j))
endif

201     continue
200   continue

write(11,*)
write(11,*)'------------------------------------------'
write(11,*)
write(11,*)'test of stack r/t j1 symmetries...'
write(11,*)'(propagating waves only)...'

write(11,*)'diagonal constants check...'
write(11,*)(diag1(i+3,1),i=1,3)
write(11,*)(diag1(i,nlay+1-ilay),i=1,3)
c1000    format(3(' (',d9.3,',',d9.3,')'))

write(11,*)
write(11,*)'11 block...'
write(11,*)wrk5(1,1)
do 300 i=1,3
write(11,1001)(wrk5(i,j),j=1,3)
300     continue
1001    format(3(' (',d9.3,',',d9.3,')'))

write(11,*)
write(11,*)'12 block...'
write(11,*)wrk6(1,1)
do 301 i=1,3
write(11,1001)(wrk6(i,j),j=1,3)
301     continue

write(11,*)
write(11,*)'21 block...'
write(11,*)wrk7(1,1)
do 302 i=1,3
write(11,1001)(wrk7(i,j),j=1,3)
302     continue

write(11,*)
write(11,*)'22 block...'
write(11,*)wrk8(1,1)
do 303 i=1,3
write(11,1001)(wrk8(i,j),j=1,3)
303     continue

write(11,*)
write(11,*)'in F format...'
write(11,*)'11 block...'
write(11,*)wrk5(1,1)
do 400 i=1,3
write(11,1002)(wrk5(i,j),j=1,3)
400     continue
1002    format(3(' (',f8.1,',',f8.1,')'))

write(11,*)
write(11,*)'12 block...'
write(11,*)wrk6(1,1)
do 401 i=1,3
write(11,1002)(wrk6(i,j),j=1,3)
401     continue

write(11,*)
write(11,*)'21 block...'
write(11,*)wrk7(1,1)
do 402 i=1,3
write(11,1002)(wrk7(i,j),j=1,3)
402     continue

write(11,*)
write(11,*)'22 block...'
write(11,*)wrk8(1,1)
do 403 i=1,3
write(11,1002)(wrk8(i,j),j=1,3)
403     continue

symacc=tst1/tst2
write(11,*)
write(11,*)'measure of j1 accuracy... symacc = ',symacc
write(11,*)
write(11,*)'------------------------------------------'
write(11,*)

return
end
c
c........|.........|.........|.........|.........|.........|.........|..
c
subroutine j2sym(nlay,ilay,evec,j2cnst)

c Copyright (c) 1996 C. J. Thomson.

c Sets up 'constants' for use in J_2 symmetry test
c (horizontal mirror plane media, see theory section (6)).

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

complex*16 evec(6,6),j2cnst(3,3)
complex*16 wrk1(6,6),wrk2(6,6)

iwrite=0
if(iwrite.eq.1)then
write(11,*)
write(11,*)'------------------------------------------'
write(11,*)
write(11,*)'in j2sym, nlay=',nlay,'  ilay=',ilay
endif

do 100 i=1,3
do 101 j=1,3

if(i.lt.3)then
wrk1(i,j)= evec(i+3,j)
wrk1(i+3,j)= -evec(i,j)
wrk1(i,j+3)= evec(i+3,j+3)
wrk1(i+3,j+3)= -evec(i,j+3)
else
wrk1(i,j)= -evec(i+3,j)
wrk1(i+3,j)= evec(i,j)
wrk1(i,j+3)= -evec(i+3,j+3)
wrk1(i+3,j+3)= evec(i,j+3)
endif

101     continue
100   continue

do 200 i=1,6
do 201 j=1,6

wrk2(i,j)=(0.d0,0.d0)

do 202 k=1,6

wrk2(i,j)=wrk2(i,j)+evec(k,i)*wrk1(k,j)

202       continue

201     continue
200   continue

if(iwrite.eq.1)then
write(11,*)
write(11,*)'value of j2 constant...'

do 300 i=1,6

write(11,3000)(wrk2(i,j),j=1,6)

300     continue

write(11,*)'and in F format...'

do 301 i=1,6

write(11,3001)(wrk2(i,j),j=1,6)

301     continue
endif

3000  format(6(' (',d10.4,',',d10.4,')'))
3001  format(6(' (',f8.1,',',f8.1,')'))

do 400 i=1,3
do 401 j=1,3

j2cnst(i,j)= wrk2(i,j+3)

401     continue
400   continue

if(iwrite.eq.1)then
write(11,*)
write(11,*)'------------------------------------------'
write(11,*)
endif

return
end
c
c........|.........|.........|.........|.........|.........|.........|..
c
subroutine j2symtst(nlay,ilay,tus,rus,tds,rds)

c Copyright (c) 1996 C. J. Thomson.

c Performs J_2 symmtery test using previously set up
c 'constants' (theory section (6)). Media with horizontal
c mirror plane, but general anisotropy, anelasticity, allowed.

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

complex*16 tus(3,3),rus(3,3),tds(3,3),rds(3,3)
complex*16 wrk1(3,3),wrk2(3,3),wrk3(3,3),wrk4(3,3)

c
c COMMON BLOCKS to be INCLUDEd
c
c max number of layers ---    1 = half-space above or top layer
c    in the stack                 just below free surface
c                          nlay = half-space below
c
parameter (nlaymx=30)
c
c Earth model
c
common /elast/isoflg(nlaymx),thickn(nlaymx),a(3,3,3,3,nlaymx),
.              rho(nlaymx)

complex*16 evals,evecs
common /evprob/evals(6,nlaymx),evecs(6,6,nlaymx)

complex*16 q11i,q12i,q21i,q22i
complex*16 tui,rui,tdi,rdi
common /intfce/q11i(3,3,nlaymx),q12i(3,3,nlaymx),
.               q21i(3,3,nlaymx),q22i(3,3,nlaymx),
.               tui(3,3,nlaymx),rui(3,3,nlaymx),
.               tdi(3,3,nlaymx),rdi(3,3,nlaymx)

complex*16 tusi,rusi,tdsi,rdsi
common /stackrt/tusi(3,3,nlaymx),rusi(3,3,nlaymx),
.                tdsi(3,3,nlaymx),rdsi(3,3,nlaymx)

c
c symmetry testing data
c
complex*16 diag1,j2cst1,j2cstn
common /symts/diag1(6,nlaymx)
common /symts2/j2cst1(3,3),j2cstn(3,3)
c
c END OF COMMON BLOCKS
c

do 100 i=1,3
do 101 j=1,3

wrk1(i,j)=(0.d0,0.d0)
wrk2(i,j)=(0.d0,0.d0)
wrk3(i,j)=(0.d0,0.d0)
wrk4(i,j)=(0.d0,0.d0)

do 102 k=1,3

wrk1(i,j)= wrk1(i,j)+rus(k,i)*j2cstn(k,j)
wrk2(i,j)= wrk2(i,j)+rds(k,i)*j2cst1(j,k)
wrk3(i,j)= wrk3(i,j)+j2cstn(k,i)*tds(k,j)
wrk4(i,j)= wrk4(i,j)+tus(k,i)*j2cst1(j,k)

102       continue

101     continue
100   continue

tst1a=0.d0
tst1b=0.d0
tst2a=0.d0
tst2b=0.d0
tst3a=0.d0
tst3b=0.d0

do 103 i=1,3
do 104 j=1,3

tst1a=tst1a+cdabs(wrk1(i,j))
tst1b=tst1b+cdabs(wrk1(i,j)-wrk1(j,i))
tst2a=tst2a+cdabs(wrk2(i,j))
tst2b=tst2b+cdabs(wrk2(i,j)-wrk2(j,i))
tst3a=tst3a+cdabs(wrk3(i,j))+cdabs(wrk4(i,j))
tst3b=tst3b+cdabs(wrk3(i,j)-wrk4(i,j))

104     continue
103   continue

write(11,*)
write(11,*)'------------------------------------------'
write(11,*)
write(11,*)'test of stack r/t j2 symmetries...'
write(11,*)'(TI media only)...'

c......................................................................|

write(11,*)'wrk1...'
do 201 i=1,3
write(11,3000)(wrk1(i,j),j=1,3)
201       continue

write(11,*)'wrk2...'
do 202 i=1,3
write(11,3000)(wrk2(i,j),j=1,3)
202       continue

write(11,*)'wrk3...'
do 203 i=1,3
write(11,3000)(wrk3(i,j),j=1,3)
203       continue

write(11,*)'wrk4...'
do 204 i=1,3
write(11,3000)(wrk4(i,j),j=1,3)
204       continue

write(11,*)
write(11,*)'and in F format...'
write(11,*)'wrk1...'
do 301 i=1,3
write(11,3001)(wrk1(i,j),j=1,3)
301       continue

write(11,*)'wrk2...'
do 302 i=1,3
write(11,3001)(wrk2(i,j),j=1,3)
302       continue

write(11,*)'wrk3...'
do 303 i=1,3
write(11,3001)(wrk3(i,j),j=1,3)
303       continue

write(11,*)'wrk4...'
do 304 i=1,3
write(11,3001)(wrk4(i,j),j=1,3)
304       continue
3000  format(3(' (',d10.4,',',d10.4,')'))
3001  format(3(' (',f8.1,',',f8.1,')'))

symac1=tst1b/tst1a
symac2=tst2b/tst2a
symac3=2.d0*tst3b/tst3a
write(11,*)
write(11,*)'measure of j2 accuracy... symac1 = ',symac1
write(11,*)'measure of j2 accuracy... symac2 = ',symac2
write(11,*)'measure of j2 accuracy... symac3 = ',symac3
write(11,*)
write(11,*)'------------------------------------------'
write(11,*)

return
end
c
c........|.........|.........|.........|.........|.........|.........|..
c
subroutine jprtst(evalsf,evecsf,nlay,
.               tus,rus,tds,rds,tuspr,ruspr,tdspr,rdspr)

c Copyright (c) 1996 C. J. Thomson.

c Performs the J-primed symmetry test (theory section (6), eqn. (6.16)).
c This is the most general symmetry. See comments below regarding the
c effect of column ordering in the e-vector matrices.

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

complex*16 tus(3,3),rus(3,3),tds(3,3),rds(3,3)
complex*16 ttus(3,3),trus(3,3),ttds(3,3),trds(3,3)

complex*16 tuspr(3,3),ruspr(3,3),tdspr(3,3),rdspr(3,3)
complex*16 ttuspr(3,3),truspr(3,3),ttdspr(3,3),trdspr(3,3)
complex*16 evalsf(6,2),evecsf(6,6,2)

complex*16 diagt(6,6),diagb(6,6)
complex*16 tstru1,tsttd1,tsttu1,tstrd1

c
c COMMON BLOCKS to be INCLUDEd
c
c max number of layers ---    1 = half-space above or top layer
c    in the stack                 just below free surface
c                          nlay = half-space below
c
parameter (nlaymx=30)
c
c Earth model
c
common /elast/isoflg(nlaymx),thickn(nlaymx),a(3,3,3,3,nlaymx),
.              rho(nlaymx)
c
complex*16 evals,evecs
common /evprob/evals(6,nlaymx),evecs(6,6,nlaymx)

complex*16 q11i,q12i,q21i,q22i
complex*16 tui,rui,tdi,rdi
common /intfce/q11i(3,3,nlaymx),q12i(3,3,nlaymx),
.               q21i(3,3,nlaymx),q22i(3,3,nlaymx),
.               tui(3,3,nlaymx),rui(3,3,nlaymx),
.               tdi(3,3,nlaymx),rdi(3,3,nlaymx)

complex*16 tusi,rusi,tdsi,rdsi
common /stackrt/tusi(3,3,nlaymx),rusi(3,3,nlaymx),
.                tdsi(3,3,nlaymx),rdsi(3,3,nlaymx)

c
c symmetry testing data
c
complex*16 diag1,j2cst1,j2cstn
common /symts/diag1(6,nlaymx)
common /symts2/j2cst1(3,3),j2cstn(3,3)
c
c END OF COMMON BLOCKS
c

iwrite=1
if(iwrite.eq.1)then
write(11,*)
write(11,*)'-----------------------------------------'
write(11,*)
write(11,*)'in jprtst -- testing reversal symmetries'
endif

c set up matrices which should come out to be diagonal -- or ....

c Important Note. The up/down columns of the e-vector matrix
c are interchanged for the reversed system, at least in the
c theoretical description (papers, notes). However, in the code
c the columns of the e-vector matrix are always down first, then
c up (set eg in subroutine vertsl2) irrespective of the sign of the
c slowness. Therefore, in this code the so-called diagonal
c matrices actually come out to be trailing diagonal. This causes
c no particular problem so long as one is aware of it.
c
c Furthermore, the sorting in, eg, vertsl2 does not ensure that
c the columns within the up/down blocks are always in the order
c P, S1, S2. They may be interchanged and this means that the
c matrix d-primed in the notes may be block diagonal (with 1's and
c zeros) not truly diagonal.

do 100 i=1,3
do 200 j=1,3

diagt(i,j)=dcmplx(0.d0,0.d0)
diagt(i,j+3)=dcmplx(0.d0,0.d0)
diagt(i+3,j)=dcmplx(0.d0,0.d0)
diagt(i+3,j+3)=dcmplx(0.d0,0.d0)

diagb(i,j)=dcmplx(0.d0,0.d0)
diagb(i,j+3)=dcmplx(0.d0,0.d0)
diagb(i+3,j)=dcmplx(0.d0,0.d0)
diagb(i+3,j+3)=dcmplx(0.d0,0.d0)

do 300 l=1,3

diagt(i,j) = diagt(i,j)
.                      +evecs(l,i,1)*evecsf(l+3,j,1)
.                        -evecs(l+3,i,1)*evecsf(l,j,1)
diagt(i,j+3) = diagt(i,j+3)
.                      +evecs(l,i,1)*evecsf(l+3,j+3,1)
.                        -evecs(l+3,i,1)*evecsf(l,j+3,1)
diagt(i+3,j) = diagt(i+3,j)
.                      +evecs(l,i+3,1)*evecsf(l+3,j,1)
.                        -evecs(l+3,i+3,1)*evecsf(l,j,1)
diagt(i+3,j+3) = diagt(i+3,j+3)
.                      +evecs(l,i+3,1)*evecsf(l+3,j+3,1)
.                        -evecs(l+3,i+3,1)*evecsf(l,j+3,1)
diagb(i,j) = diagb(i,j)
.                      +evecs(l,i,nlay)*evecsf(l+3,j,2)
.                        -evecs(l+3,i,nlay)*evecsf(l,j,2)
diagb(i,j+3) = diagb(i,j+3)
.                      +evecs(l,i,nlay)*evecsf(l+3,j+3,2)
.                        -evecs(l+3,i,nlay)*evecsf(l,j+3,2)
diagb(i+3,j) = diagb(i+3,j)
.                      +evecs(l,i+3,nlay)*evecsf(l+3,j,2)
.                        -evecs(l+3,i+3,nlay)*evecsf(l,j,2)
diagb(i+3,j+3) = diagb(i+3,j+3)
.                      +evecs(l,i+3,nlay)*evecsf(l+3,j+3,2)
.                        -evecs(l+3,i+3,nlay)*evecsf(l,j+3,2)

300    continue

200  continue
100 continue

if(iwrite.eq.1)then

write(11,*)'top diagonal(?) constant'
do 400 i=1,6
write(11,3000)(diagt(i,j),j=1,6)
3000 format(3('(',d10.4,',',d10.4,')'),/,
.        5X,3('(',d10.4,',',d10.4,')'))
400   continue

write(11,*)
write(11,*)'bottom diagonal(?) constant'
do 500 i=1,6
write(11,3000)(diagb(i,j),j=1,6)
500   continue

endif

c r/t coefficient symmetries

tstru1=dcmplx(0.d0,0.d0)
tsttd1=dcmplx(0.d0,0.d0)
tsttu1=dcmplx(0.d0,0.d0)
tstrd1=dcmplx(0.d0,0.d0)

tstru2=0.d0
tsttd2=0.d0
tsttu2=0.d0
tstrd2=0.d0

do 600 i=1,3
do 700 j=1,3

trus(i,j) = dcmplx(0.d0,0.d0)
ttds(i,j) = dcmplx(0.d0,0.d0)
ttus(i,j) = dcmplx(0.d0,0.d0)
trds(i,j) = dcmplx(0.d0,0.d0)

truspr(i,j) = dcmplx(0.d0,0.d0)
ttdspr(i,j) = dcmplx(0.d0,0.d0)
ttuspr(i,j) = dcmplx(0.d0,0.d0)
trdspr(i,j) = dcmplx(0.d0,0.d0)

do 701 k=1,3

trus(i,j) = trus(i,j) + diagb(i+3,k)*rus(k,j)
ttds(i,j) = ttds(i,j) + diagb(i+3,k)*tds(k,j)
ttus(i,j) = ttus(i,j) + diagt(i,k+3)*tus(k,j)
trds(i,j) = trds(i,j) + diagt(i,k+3)*rds(k,j)

truspr(i,j) = truspr(i,j) + ruspr(k,i)*diagb(k,j+3)
ttdspr(i,j) = ttdspr(i,j) + tuspr(k,i)*diagt(k+3,j)
ttuspr(i,j) = ttuspr(i,j) + tdspr(k,i)*diagb(k,j+3)
trdspr(i,j) = trdspr(i,j) + rdspr(k,i)*diagt(k+3,j)

701   continue

tstru1= tstru1+truspr(i,j)+trus(i,j)
tstru2= tstru2+cdabs(truspr(i,j))+cdabs(trus(i,j))

tsttd1= tsttd1+ttdspr(i,j)-ttds(i,j)
tsttd2= tsttd2+cdabs(ttdspr(i,j))+cdabs(ttds(i,j))

tsttu1= tsttu1+ttuspr(i,j)-ttus(i,j)
tsttu2= tsttu2+cdabs(ttuspr(i,j))+cdabs(ttus(i,j))

tstrd1= tstrd1+trdspr(i,j)+trds(i,j)
tstrd2= tstrd2+cdabs(trdspr(i,j))+cdabs(trds(i,j))

700  continue
600 continue

if(iwrite.eq.1)then

write(11,*)'reversed and scaled r/t matrices'
write(11,*)'ru-primed'
write(11,3001)((truspr(i,j),j=1,3),i=1,3)
3001 format(3('(',d10.4,',',d10.4,')'))
write(11,*)'ru'
write(11,3001)((trus(i,j),j=1,3),i=1,3)
write(11,*)'td-primed'
write(11,3001)((ttdspr(i,j),j=1,3),i=1,3)
write(11,*)'td'
write(11,3001)((ttds(i,j),j=1,3),i=1,3)
write(11,*)'tu-primed'
write(11,3001)((ttuspr(i,j),j=1,3),i=1,3)
write(11,*)'tu'
write(11,3001)((ttus(i,j),j=1,3),i=1,3)
write(11,*)'rd-primed'
write(11,3001)((trdspr(i,j),j=1,3),i=1,3)
write(11,*)'rd'
write(11,3001)((trds(i,j),j=1,3),i=1,3)

accru=cdabs(tstru1)/tstru2
acctd=cdabs(tsttd1)/tsttd2
acctu=cdabs(tsttu1)/tsttu2
accrd=cdabs(tstrd1)/tstrd2

write(11,*)
write(11,*)'accuracy measures:'
write(11,*)'acc-ru =',accru
write(11,*)'acc-td =',acctd
write(11,*)'acc-tu =',acctu
write(11,*)'acc-rd =',accrd

endif

return
end
c
c........|.........|.........|.........|.........|.........|.........|..
c
subroutine frees(evec,tds,rds,det1,det2)

c Copyright (c) 1996 C. J. Thomson.

c Computes the determinants associated with the free--surface
c calculation (see theory notes section (5)). Not a required
c subroutine.

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

complex*16 evec(6,6)
complex*16 tds(3,3),rds(3,3)
complex*16 wrk1(3,3),det1,det2

do 100 i=1,3
do 200 j=1,3

wrk1(i,j)=(0.d0,0.d0)
do 300 k=1,3
wrk1(i,j)=wrk1(i,j)+evec(i+3,k+3)*rds(k,j)
300       continue

wrk1(i,j)=evec(i+3,j)+wrk1(i,j)

200     continue
100   continue

det1=wrk1(1,1)*(wrk1(2,2)*wrk1(3,3)-wrk1(2,3)*wrk1(3,2))
.     -wrk1(1,2)*(wrk1(2,1)*wrk1(3,3)-wrk1(2,3)*wrk1(3,1))
.     +wrk1(1,3)*(wrk1(2,1)*wrk1(3,2)-wrk1(2,2)*wrk1(3,1))

det2=tds(1,1)*(tds(2,2)*tds(3,3)-tds(2,3)*tds(3,2))
.     -tds(1,2)*(tds(2,1)*tds(3,3)-tds(2,3)*tds(3,1))
.     +tds(1,3)*(tds(2,1)*tds(3,2)-tds(2,2)*tds(3,1))

return
end
c
c........|.........|.........|.........|.........|.........|.........|..
c
subroutine write1(ilay,isoflg,px,py,eval,evec)

c Copyright (c) 1996 C. J. Thomson.

implicit real*8(a-h,o-z)
complex*16 px,py
complex*16 eval(6),evec(6,6)

write(11,*)
write(11,*)'---------------------------------------------'
write(11,*)
write(11,*)
write(11,*)'layer ',ilay,'     px=',px,'   py=',py
if(isoflg.eq.0)write(11,*)'anisotropic layer'
if(isoflg.eq.1)write(11,*)'isotropic layer'
write(11,*)
write(11,*)'eigenvalues or vertical slownesses...'
write(11,*)eval
write(11,*)'and eigenvectors...'
do 1000 j=1,6
write(11,*)(evec(i,j),i=1,6)
write(11,*)
1000   continue
write(11,*)
write(11,*)'---------------------------------------------'
write(11,*)

return
end
c
c........|.........|.........|.........|.........|.........|.........|..
c
subroutine write2(ilay,evec,eveci)

c Copyright (c) 1996 C. J. Thomson.

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

complex*16 evec(6,6),eveci(6,6),testm(6,6)

call xmult6(evec,eveci,testm)

write(11,*)
write(11,*)'---------------------------------------------'
write(11,*)
write(11,*)
write(11,*)'test of inverse e-vector matrix, layer ',ilay
do 55 i=1,6
write(11,3000)(testm(i,j),j=1,6)
55  continue
write(11,*)'and in F format...'
do 56 i=1,6
write(11,3001)(testm(i,j),j=1,6)
56  continue

call xmult6(eveci,evec,testm)

write(11,*)
write(11,*)'second test of inverse...'
do 57 i=1,6
write(11,3000)(testm(i,j),j=1,6)
57  continue
write(11,*)'and in F format...'
do 58 i=1,6
write(11,3001)(testm(i,j),j=1,6)
58  continue
3000 format(6(' (',d10.4,',',d10.4,')'))
3001 format(6(' (',f5.2,',',f5.2,')'))
write(11,*)
write(11,*)'---------------------------------------------'
write(11,*)

return
end
c
c........|.........|.........|.........|.........|.........|.........|..
c
subroutine write3(ilay,tu,ru,td,rd,evec1,evec2)

c Copyright (c) 1996 C. J. Thomson.

implicit real*8(a-h,o-z)
complex*16 tu(3,3),ru(3,3),td(3,3),rd(3,3)
complex*16 testm1(6,6),testm2(6,6),evec1(6,6),evec2(6,6)

write(11,*)
write(11,*)'---------------------------------------------'
write(11,*)
write(11,*)
write(11,*)'r/t matrices ... interface ',ilay
write(11,*)'tu'
do 1021 i=1,3
write(11,3100)(tu(i,j),j=1,3)
1021 continue
write(11,*)'ru'
do 1022 i=1,3
write(11,3100)(ru(i,j),j=1,3)
1022 continue
write(11,*)'td'
do 1023 i=1,3
write(11,3100)(td(i,j),j=1,3)
1023 continue
write(11,*)'rd'
do 1024 i=1,3
write(11,3100)(rd(i,j),j=1,3)
1024 continue

3100 format(3(' (',d9.3,',',d9.3,')'))

write(11,*)
write(11,*)'testing b.c.s ... interface',ilay

do 3201 i=1,6
do 3202 j=1,6
testm1(i,j)=(0.d0,0.d0)
3202  continue
3201  continue

do 3200 i=1,3
do 3205 j=1,3
do 3210 k=1,3
testm1(i,j)=testm1(i,j)+evec1(i,k)*ru(k,j)
testm1(i,j+3)=testm1(i,j+3)+evec1(i,k)*td(k,j)
testm1(i+3,j)=testm1(i+3,j)+evec1(i+3,k)*ru(k,j)
testm1(i+3,j+3)=testm1(i+3,j+3)+evec1(i+3,k)*td(k,j)
3210    continue
testm1(i,j)=testm1(i,j) + evec1(i,j+3)
testm1(i+3,j)=testm1(i+3,j) + evec1(i+3,j+3)
3205   continue
3200  continue

write(11,*)testm1(1,1),testm1(2,2),testm1(6,6)

write(11,*)'lower stress-displacement matrix...'
do 3215 i=1,6
write(11,3000)(testm1(i,j),j=1,6)
3215  continue
3000  format(6(' (',d10.4,',',d10.4,')'))

do 3203 i=1,6
do 3204 j=1,6
testm2(i,j)=(0.d0,0.d0)
3204  continue
3203  continue

do 3300 i=1,3
do 3305 j=1,3
do 3310 k=1,3
testm2(i,j)=testm2(i,j)+evec2(i,k+3)*tu(k,j)
testm2(i,j+3)=testm2(i,j+3)+evec2(i,k+3)*rd(k,j)
testm2(i+3,j)=testm2(i+3,j)+evec2(i+3,k+3)*tu(k,j)
testm2(i+3,j+3)=testm2(i+3,j+3)+evec2(i+3,k+3)*rd(k,j)
3310    continue
testm2(i,j+3)=testm2(i,j+3) + evec2(i,j)
testm2(i+3,j+3)=testm2(i+3,j+3) + evec2(i+3,j)
3305   continue
3300  continue

write(11,*)'upper stress-displacement matrix...'

write(11,*)testm2(1,1),testm2(2,2),testm2(6,6)

do 3315 i=1,6
write(11,3000)(testm2(i,j),j=1,6)
3315  continue

bcacc1=0.d0
bcacc2=0.d0
do 4000 i=1,6
do 4001 j=1,6
bcacc1=bcacc1+cdabs(testm1(i,j)-testm2(i,j))
bcacc2=bcacc2+cdabs(testm1(i,j))+cdabs(testm2(i,j))
4001  continue
4000  continue

bcacc=bcacc1/bcacc2
write(11,*)
write(11,*)'overall bc accuracy measure...diff, sum, bcacc ='
write(11,*)bcacc1,bcacc2,bcacc
write(11,*)
write(11,*)'---------------------------------------------'
write(11,*)

return
end
c
c........|.........|.........|.........|.........|.........|.........|..
c
subroutine write4(tus,rus,tds,rds)

c Copyright (c) 1996 C. J. Thomson.

complex*16 tus(3,3),rus(3,3),tds(3,3),rds(3,3)

write(11,*)
write(11,*)'---------------------------------------------'
write(11,*)
write(11,*)'final r/t matrices ...'
write(11,*)'tus'
do 9021 i=1,3
write(11,3100)(tus(i,j),j=1,3)
9021   continue
write(11,*)'rus'
do 9022 i=1,3
write(11,3100)(rus(i,j),j=1,3)
9022   continue
write(11,*)'tds'
do 9023 i=1,3
write(11,3100)(tds(i,j),j=1,3)
9023   continue
write(11,*)'rds'
do 9024 i=1,3
write(11,3100)(rds(i,j),j=1,3)
9024   continue
3100   format(1(' (',d15.9,',',d15.9,')'))

return
end
c
c........|.........|.........|.........|.........|.........|.........|..
c
C                                                                 ```