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 All rights reserved by the author.

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.
c All rights reserved by the author.

      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.
c All rights reserved by the author.

      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.
c All rights reserved by the author.

      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.
c All rights reserved by the author.

      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 All rights reserved by the author.

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.
c All rights reserved by the author.

      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.
c All rights reserved by the author.

      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 All rights reserved by the author.

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)

      call xadd3(q11,xhelp,td)

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

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

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 All rights reserved by the author.

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 All rights reserved by the author.

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
      subroutine xadd3(a,b,c)

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

c
c Complex matrix addition a+b=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.
c All rights reserved by the author.

      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
      subroutine addit(nlay,ilay,tus,rus,tds,rds,
     .                  tu,ru,td,rd)

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

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,*)
        write(11,*)'inside addit on recursion'

        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 xadd3(ximat,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 addition rules...     with symbolic notation comments

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)
      call xadd3(rd,rdsw,rdsw)

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

c     tdsw = tds * [ ximat + ru*reverbi*rds ] * td
      call xmult3(reverbi,rds,tdsw)
      call xmult3(ru,tdsw,tdsw)
      call xadd3(ximat,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 All rights reserved by the author.

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.
c All rights reserved by the author.

      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 All rights reserved by the author.

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 All rights reserved by the author.

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 All rights reserved by the author.

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 All rights reserved by the author.

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 All rights reserved by the author.

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.
c All rights reserved by the author.

        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.
c All rights reserved by the author.

      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.
c All rights reserved by the author.

      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.
c All rights reserved by the author.

        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