c
c........|.........|.........|.........|.........|.........|.........|..
c
c Generic program to run routine Rmatrix, which finds the R/T
c coefficient matrices for a stack of anisotropic layers.
c
c Except where otherwise acknowldeged, rmatrix was written by
c C. J. Thomson.
c
c This version allows only real elasticity.
c
c C. J. Thomson, started August 1993.
c Modified for complex horizontal slowness and frequency, February 1995.
c
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 Copyright (c) 1996 C. J. Thomson.
c All rights reserved by the author.

      program rmatrix1

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

      complex*16 px,py,omega

c overal r/t matrices for stack

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

c overal r/t matrices for stack with reversed slowness (primed)
      complex*16 tuspr(3,3),ruspr(3,3),tdspr(3,3),rdspr(3,3)
c and the e-values and vectors for symmetry testing
      complex*16 evalsf(6,2),evecsf(6,6,2)

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

c open i/o files

      open(unit=10,file='rmatrix.inp',form='FORMATTED')
      open(unit=11,file='rmatrix.out',form='FORMATTED')

c-LK: Originally the second line of the input data is now read here
c-LK: as the first line of the input data:
      read(10,*)omega,px,py
c-LK
c-LK: Input logical unit number 10 added as the first argument,
c-LK: arguments PX,PY,OMEGA discarded:
      call input(10,nlay,ifree)


c call rmatrix for px, py pair just selected... and free surface
c not considered here

       iter=1
       print *,'iter=',iter
       call rmatrix(px,py,omega,nlay,ifree,tus,rus,tds,rds,iter)

       iter=2
       print *,'iter=',iter
       call rmatrix(px,py,omega,nlay,ifree,tus,rus,tds,rds,iter)

c test reversal symmetry ...
c first store eigenvector matrices at top and bottom

       call trnsf6(evalsf(1,1),evals(1,1),evecsf(1,1,1),evecs(1,1,1))
       call trnsf6(evalsf(1,2),evals(1,nlay),
     .                   evecsf(1,1,2),evecs(1,1,nlay))

       iter=1
       px=-px
       py=-py
       call rmatrix(px,py,omega,nlay,ifree,tuspr,ruspr,tdspr,rdspr,iter)

       call jprtst(evalsf,evecsf,nlay,
     .               tus,rus,tds,rds,tuspr,ruspr,tdspr,rdspr)

       iter=2
       call rmatrix(px,py,omega,nlay,ifree,tuspr,ruspr,tdspr,rdspr,iter)

       call jprtst(evalsf,evecsf,nlay,
     .               tus,rus,tds,rds,tuspr,ruspr,tdspr,rdspr)


      close(10)
      close(11)


      stop
      end
c
c........|.........|.........|.........|.........|.........|.........|..
c