C 
C    Auxiliary main program RTCOEF
C    (for testing of routine COEF52)
C    *******************************
C         For testing purposes, a brief main program RTCOEF is included
C    to control the computation of R/T coefficients using the routine
C    COEF52. Main program first calls COEF52 to read the data for the
C    model, and, after this, it reads two records specifying the
C    required computations. The computations are performed either in an
C    angle loop, or in a frequency loop. For any required angle of
C    incidence and frequency, routine COEF52 is called to compute the
C    relevant R/T coefficient.
C
C    The complete set of input data is as follows:
C    ********************************************
C    1) One record: NZ
C         NZ ... number of layers, including both halfspaces
C    2) NZ records: VP,VS,RHO,QP,QS,D
C         Parameters of the layers; one record for one layer. The
C         first record corresponds to the halfspace with the incident
C         wave (first halfspace), the last record to the NZ-th layer,
C         that is to the second halfspace.
C    3) One record: FREF
C         FREF ... reference frequency (in Hz).
C    4) One record: NC,NH,NQ,NF,NA
C         NC ... type of R/T coefficient
C         NH ... NH=0...the second halfspace is solid or liquid
C                NH=1...the second halfspace is a vacuum
C         NQ ... NQ=0... model is non-dissipative
C                NQ=1... model is dissipative
C         NF ... number of frequencies in the frequency loop.
C                For NF=1, the frequency loop is not considered
C         NA ... number of angles of incidence in the angle loop.
C                For NA=1, the angle loop is not considered
C     5) One record: FMIN,DF,AMIN,DA,GAMMA
C         FMIN,DF ... specify frequency loop, in Hz
C         AMIN,DA ... specify angle loop, in degrees
C         GAMMA   ... attenuation angle, in degrees
C    The records 4 and 5 can be repeated any number of times. The
C    computation finishes if NC=0 is used in record 4.
C         If we wish to compute the R/T coefficients in an angle loop
C    (for one given frequency FMIN), we use NF=1. If we wish to
C    compute the R/T coefficients in a frequency loop (for one
C    given angle of incidence AMIN), we use NA=1. Use always NF=1
C    or NA=1.
C         The input data are stored in the file input.dat.
C         The results of computations are stored in the file
C    output.dat. For convenience, all input data are first stored
C    in output.dat (even those for model). Then, after an empty
C    line, the results of computations are stored. Each line diplays:
C    ANGLE,FREQ,RMOD,RPHASE,RMOD0,RPH0.
C
C    Test examples
C    *************
C    Four test examples are included. See the description in the file
C    'rtcoef.htm'.
C
C    The test data are located in subdirectory
C    rtcoef
C    of package DATA.
C    The test examples may be executed by command
C      perl go.pl rtcoef.h
C    running the demonstration history file
C    rtcoef.h.
C
C    References
C    **********
C    Brokesova,J. (2000). Reflection/transmission coefficients at a
C         plane interface in dissipative and non-dissipative media:
C         A comparison. J.Comput.Acoustics, 9,623 -641.
C    Brokesova,J., and Cerveny,V. (1998). Inhomogeneous plane waves
C         in dissipative, isotropic and anisotropic media. Reflection/
C         transmission coefficients. In Seismic waves in complex 3-D
C         structures, Report No. 7, pp. 57 - 146. Department of
C         Geophysics, Charles University, Prague.
C    Cerveny,V. (1989). Synthetic body wave seismograms for laterally
C         varying media containing thin transition layers. Geophys. J.
C         Int., 99, 331-349,
C    Cerveny,V. (2001). Seismic ray theory. Cambridge Univ. Press,
C         Cambridge.
C    Muller,G. (1985). The reflectivity method. A tutorial. J.Geophys.,
C         58, 153-174.
C
C
C    Consortium Project "Seismic Waves in Complex 3-D Structures"
C    Praha, April 2003
C    J.Brokesova, V.Cerveny
C
************************************************************************
c
c    program rtcoef
c    **************
c
c    Auxiliary program rtcoef is designed for computations
c    of R/T coefficients of inhomogeneous P, SV and SH plane wave
c    at stack of layers between two isotropic anelastic halfspaces.
c    It uses the routine coef52.
c
      open(7,file='input.dat')
      open(9,file='output.dat')
c
c     reading the model
      call coef52(0,0,0,0,7,9,0.,0.,0.,0.,0.,0.,0.)
c
c     reading the data for coef52
    1 read(7,100)nc,nh,nq,nf,na
      write(9,100)nc,nh,nq,nf,na
      if(nc.eq.0.or.nc.eq.11.or.nc.eq.12.or.nc.gt.14)go to 10
      read(7,101)fmin,df,amin,da,gamma
      write(9,101)fmin,df,amin,da,gamma
      write(9,102)
c
c     angle loop
      if(na.eq.1)go to 5
      freq=fmin
      do 2 i=1,na
      angle=amin+(i-1)*da
      call coef52(nc,nh,nq,1,7,9,angle,gamma,freq,rmod,rphase,
     1rmod0,rph0)
      write(9,101)angle,freq,rmod,rphase,rmod0,rph0
    2 continue
      go to 10
c
c     frequency loop
    5 angle=amin
      do 6 i=1,nf
      freq=fmin+(i-1)*df
      call coef52(nc,nh,nq,1,7,9,angle,gamma,freq,rmod,rphase,
     1rmod0,rph0)
      write(9,101)angle,freq,rmod,rphase,rmod0,rph0
    6 continue
c
    7 go to 1
  100 format(6i5)
  101 format(8f10.3)
  102 format (/)
   10 continue
      stop
      end
c
************************************************************************
c
      INCLUDE 'coef52.for'
c     coef52.for
c
************************************************************************
c