C 
C    Routine COEF52
C    **************
C
C         Routine COEF52 is designed for the computation of the
C    displacement reflection/transmission coefficients (R/T
C    coefficients) of inhomogeneous P, SV and SH plane waves at a
C    stack of homogeneous isotropic dissipative layers between two
C    homogeneous isotropic dissipative halfspaces. The stack of
C    layers may be reduced to a single interface.
C         As special cases, it is also possible to compute the R/T
C    coefficients of pressure waves at a stack of liquid layers
C    between two liquid halfspaces, the R/T coefficients at a
C    surficial stack of layers, and the conversion coefficients
C    (receiver functions) at a free surface of the surficial
C    stack of layers. In all these cases, analogous R/T coefficients
C    for non-dissipative media may be also computed.
C         The routine is written in Fortran 77 programming language.
C
C    Calling:
C    ********
C    CALL COEF52(NC,NH,NQ,NUM,LIN,LOU,ANGLE,GAMMA,FREQ,RMOD,RPHASE,
C    RMOD0,RPH0)
C    Input parameters:
C    NC ..... Type of the 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    NUM .... NUM=0 ... reading the model, no computations
C             NUM=1 ... computations of R/T coefficients
C    LIN .... Number of input file logical unit
C    LOU .... Number of output file logical unit
C    ANGLE .. Angle of incidence (in degrees)
C    GAMMA .. Attenuation angle of the incident inhomogeneous
C             plane wave (in degrees)
C    FREQ ... Frequency (in Hz)
C    Output parameters:
C    RMOD ... Modul of computed R/T coefficient
C    RPHASE.. Phase of computed R/T coefficient
C    RMOD0 .. Modul of analogous R/T coefficient for a single interface
C             between two non-dissipative halfspaces
C    RPH0 ... Phase of analogous R/T coefficient.
C    For NUM=0, all other parameters, with the exception of LIN and
C    LOU, may be arbitrary.
C
C    Specification of the model:
C    ***************************
C         The number of layers in the model (including both
C     halfspaces) is specified by NZ. For NZ=2, the two halfspaces
C     are in contact and the transition layer reduces to a single
C     interface. The layer number 1 corresponds to the incidence
C     (first) halfspace, the layer number NZ to the second halfspace.
C     For each layer (including halfspaces), the following
C     real-valued parameters should be specified:
C          VP ... P velocity (in km/s),
C          VS ... S velocity (in km/s),
C          RHO .. density (in g/cm**3),
C          QP ... quality factor of P waves,
C          QS ... quality factor of S waves,
C          D  ... thickness of the layer (in km).
C         The quality factors QP and QS of P and S waves are assumed
C    to be independent of frequency FREQ (constant-Q model).
C         The real-valued velocities VP and VS of P and S waves are
C    material properties of the model for the reference frequency
C    FREF. They correspond to the real parts of the Lame's elastic
C    moduli for the reference frequency FREF. For known VP(FREF),
C    VS(FREF), QP and QS, the frequency-dependent complex valued
C    velocities CVP(FREQ) and CVS(FREQ) of P and S waves for the
C    constant Q model are given by relations
C    CVP(FREQ) = VP(FREF)*(1.+ln(FREQ/FREF)/(pi*QP) - i/(2.*QP)),
C    CVS(FREQ) = VS(FREF)*(1.+ln(FREQ/FREF)/(pi*QS) - i/(2.*QS)).
C    The frequency-dependent real-valued velocities VP(FREQ) and
C    VS(FREQ) equal the real parts of CVP(FREQ) and CVS(FREQ).
C    They specify the velocity dispersion.
C         In both halfspaces, the thicknesses (D(1) and D(NZ)) may
C    be taken arbitrarily; they are automatically adjusted to zero.
C         For pressure waves in liquid media, VS and QS are not used
C    in the computations. The values of VS and QS may be specified
C    arbitrarily in this case.
C         For SH waves, VP and QP are not used in the computations.
C    The values of VP and QP may be specified arbitrarily in this case.
C         For non-dissipative media (NQ=0), QP and QS are not used
C    in computations. The values of QP and QS may be specified
C    arbitrarily in this case.
C         If the second halfspace represents a vacuum (NH=1),
C    the density  in the second halfspace is automatically adjusted
C    to zero. Similarly, the velocities VP and VS are adjusted to 0.001.
C         For NUM=0, routine COEF52 reads the input data for the
C    model (including the reference frequency FREF), and does not
C    perform any computations. For NUM=1, the computation of the R/T
C    coefficients is performed, for the known model.
C
C    Type of R/T coefficients.
C    *************************
C         The type of the R/T coefficient we wish to compute is
C    specified by NC (input parameter of COEF52). The following
C    table shows the values of NC for individual  R/T coefficients:
C    a) P/SV waves, solids:
C    P1P1...1     P1S1...2     P1P2...3     P1S2...4
C    S1P1...5     S1S1...6     S1P2...7     S1S2...8
C    b) SH waves, solids:
C    S1S1...9     S1S2...10
C    c) Pressure waves, liquids:
C    P1P1...13    P1P2...14
C    For example, P1S1 (NC=2) corresponds to the reflection
C    coefficient, for incident P wave and reflected S wave.
C    Similarly, S1P2 (NC=7) corresponds to the transmission
C    coefficient, for incident S wave and transmitted P wave.
C         The second halfspace may be solid or liquid (NH=0) or
C    a vacuum (NH=1). For NH=1, the interface between the (NZ-1)-th
C    layer and the NZ-th layer (second halfspace) represents a free
C    surface, e.g. the surface of the Earth. For NH=1, the reflection
C    coefficients P1P1, P1S1, S1P1 and S1S1 have a standard meaning.
C    The  physical meaning of the transmission coefficients P1P2,
C    P1S2, S1P2 and S1S2 for NH=1 is, however, different. They
C    give the so-called conversion coefficients, also called receiver
C    functions:
C         For NH=1:   P1P2=PZ, P1S2=PX, S1P2=SZ, S1S2=SX.
C    The conversion coefficients (receiver functions) represent
C    horizontal or vertical displacement components of the free surface
C    (top of the stack of layers) due to a plane wave incident on the
C    bottom of the stack of layers from below. In the notation PZ, PX,
C    SZ and SX for the conversion coefficients, the first letter
C    specifies the incident wave (P or S), and the second letter the
C    Cartesian component of the displacement vector (X or Z). The
C    Cartesian axis X is horizontal, tangential to the free surface,
C    situated in the plane of incidence, oriented "against" the
C    direction of propagation. (For the X-axis oriented "along" the
C    direction of propagation, we must take opposite signs of PX and
C    SX.) The Cartesian axis Z is vertical, perpendicular to the
C    interface, and oriented away from the transition layer. For SH
C    waves and NH=1, the R/T coefficient S1S2 represents the conversion
C    coefficient SY. The Cartesian axis Y is horizontal, perpendicular
C    to the plane of incidence.
C
C    Incident wave:
C    **************
C         In dissipative media, the slowness vector of the incident
C    wave is complex-valued. The real part of the slowness vector
C    is called the propagation vector, and its imaginary part the
C    attenuation vector. The plane of incidence is specified by the
C    normal vector to the interface(s) and by the propagation vector.
C    The attenuation vector of the incident wave is assumed to be
C    situated in the plane of incidence. The direction of the
C    propagation vector is specified by the real-valued angle of
C    incidence ANGLE, and the direction of the attenuation vector is
C    specified by the real-valued attenuation angle GAMMA. The angles
C    ANGLE and GAMMA are input parameters of COEF52. ANGLE represents
C    the acute angle between the propagation vector of the incident
C    wave and normal to the interface. It is specified in degrees, and
C    should be in the range <0.,90.>. The attenuation angle is the
C    angle between the propagation and attenuation vectors of the
C    incident wave, is specified in degrees, and should be in the
C    range (-90.,90.). It is positive if the attenuation vector is
C    pointing to the left from the propagation vector.
C        The complex-valued polarisation vector of the incident P
C    wave is proportional to the slowness vector. The complex-valued
C    polarisation vector of the incident SV wave is perpendicular to
C    the slowness vector (in the complex sense), and its real part is
C    oriented to the left from the propagation vector. The polarisation
C    vector of the incident SH wave is perpendicular to the plane of
C    incidence. Analogous orientation  of polarisation vectors is used
C    even for R/T waves. Consequently, the R/T coefficients are
C    independent of the used coordinate system.
C         For given NC, NH, NQ, ANGLE, GAMMA, and FREQ, and for a
C    given model, the routine COEF52 returns the complex-valued
C    R/T coefficient RCOEF. RCOEF is expressed in terms of its modulus
C    RMOD and phase RPHASE. Both RMOD and RPHASE are real-valued.
C    RPHASE is specified in the range <-180.,180.>.
C         The time factor of the incident wave is exp(-i*omega*t),
C    where t is running time and omega=2.*pi*FREQ. For the time factor
C    exp(i*omega*t), it would be necessary to change the sign of
C    RPHASE.
C         For RMOD=0., RPHASE cannot be computed. Consequently,
C    COEF52 returns 'the warning' RPHASE=999 for RMOD<0.00001. In
C    plotting RPHASE, the points with RPHASE=999 should be ignored.
C         In each computation, also quantities RMOD0 and RPH0 are
C    computed, in addition to RMOD and RPHASE. They correspond to the
C    same R/T coefficient as RMOD and RPHASE, but to the plane
C    interface between the two non-dissipative halfspaces (NZ=2, NQ=0).
C
C    Routines used
C    *************
C         For matrix computations of R/T coefficients for P and SV
C    waves, routine RTMATQ is used in dissipative media, and routine
C    RTMAT in non-dissipative media. They further use routines RTQ,
C    RT, etc.
C         For matrix computation of R/T coefficients of SH waves and
C    for pressure waves in liquids, routine RTMQ2 is used in
C    dissipative media, and routine RTM2 in non-dissipative media.
C         For comparisons, R/T coefficients at a single interface
C    between two non-dissipative media (NZ=2, NQ=0), are computed
C    in routine COEF8, using explicit (non-matrix) expressions for
C    R/T coefficients.
C         Routine CRITAN tests a possible existence of nonelastic
C    discrete critical angles, and computes them if they exist.
C    Nonelastic discrete critical angles correspond to angles of
C    incidence for which both the real and imaginary parts of the
C    vertical component of the slowness vector of any R/T wave in the
C    dissipative model under consideration vanish. The R/T coefficients
C    may display anomalous jumps at the discrete critical angles. The
C    application of CRITAN removes such jumps, see Brokesova and
C    Cerveny (1998). In COEF52, routine CRITAN is applied only if NZ=2,
C    not for NZ.GT.2. Note that the anomalous behaviour of R/T
C    coefficients was also observed in the so-called Rayleigh window.
C    These cases will be subjects of further investigation.
C         Routines RTMAT and RTMATQ are modifications of analogous
C    routines, written by G. Muller, used in his reflectivity packages,
C    and described in the paper Muller (1985). The recursive algorithm
C    is used to compute the R/T coefficients. The potential R/T
C    coefficients, computed by Muller, are adjusted to displacement
C    R/T coefficients, considered here. Although G.Muller
C    considers the time convention exp(i*omega*t), the resulting phase
C    is adjusted to our time convention exp(-i*omega*t). Moreover,
C    inhomogeneous plane incident waves with arbitrary attenuation
C    angles are introduced. The routines are used in COEF52 with the
C    author's permission.
C         Routine COEF8 is a modification of an analogous routine,
C    used in the program package SEIS88 by V. Cerveny and I. Psencik.
C    See explicit analytical expressions for the coefficients in
C    Cerveny (2001). The book also displays pictures of many R/T
C    coefficients at a single interface between two non-dissipative
C    halfspaces.
C         Routine COEF52 itself is a modification of the routine FDEP,
C    used in the program package BEAM87 by V.Cerveny to study the
C    propagation of Gaussian beams in complex, 2-D, laterally varying
C    layered structures, containing thin stacks of layers. See
C    Cerveny (1989).
C         The theoretical treatment of reflection and transmission of
C    inhomogeneous plane waves at a single interface between two
C    dissipative elastic halfspaces can be found in Brokesova and
C    Cerveny (1998). The reference also offers a more detailed
C    description of routine COEF52 and presents numerous examples
C    of R/T coefficients. See also Brokesova (2001).
C
C    Demonstration program for testing
C    *********************************
C         For testing purposes, a brief main program
C    RTCOEF is included.
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
C    in Complex 3-D Structures"
C    Praha, April 2003
C    J.Brokesova, V.Cerveny
C
*********************************************************************
c
      subroutine coef52(nc,nh,nq,num,lin,lou,angle,gamma,freq,rmod,
     1rphase,rmod0,rph0)
      complex rpp,rps,rss,rsp,tpp,tps,tss,tsp,rcoef,rp,tp,rs,ts,cunit,
     1u
      common/model/nz,vp(50),vs(50),rho(50),d(50),qp(50),qs(50),fref,
     1rho2,vp2,vs2,qp2,qs2
c
      cunit=(0.,1.)
      pi=3.141593
      ang=angle*pi/180.
      gam=gamma*pi/180.
      inc=0
      if(nc.gt.4.and.nc.lt.11)inc=1
      ntype=2
      if(nc.ge.13)ntype=0
      if(nc.gt.8.and.nc.lt.13)ntype=1
      if(num.eq.1)go to 20
c
c     reading the model
      read(lin,100)nz
      write(lou,100)nz
      do 1 i=1,nz
      read(lin,101)vp(i),vs(i),rho(i),qp(i),qs(i),d(i)
      write(lou,101)vp(i),vs(i),rho(i),qp(i),qs(i),d(i)
    1 continue
      read(lin,101)fref
      write(lou,101)fref
      rho2=rho(nz)
      vp2=vp(nz)
      vs2=vs(nz)
      d(1)=0.
      d(2)=0.
      return
c
   20 cvr=0.
      if(nq.eq.1)cvr=alog(freq/fref)/pi
      if(nh.eq.1)go to 21
      vp(nz)=vp2
      vs(nz)=vs2
      rho(nz)=rho2
      go to 24
   21 vp(nz)=0.001
      vs(nz)=0.001
      rho(nz)=0.
   24 continue
c
c     computation of the tangential component of the slowness vector
c     of the incident wave
      if(inc.eq.1)go to 22
      v1=vp(1)
      q1=qp(1)
      urr=1./v1
      if(nq.eq.0)go to 23
      v1=vp(1)*(1.+cvr/qp(1))
      go to 23
   22 v1=vs(1)
      q1=qs(1)
      urr=1./v1
      if(nq.eq.0)go to 23
      v1=vs(1)*(1.+cvr/qs(1))
   23 vv1=v1*v1
      par=1/v1
      ur=par*sin(ang)
      urr=urr*sin(ang)
      if(nq.eq.0)go to 26
      gaux=1./(q1*q1*cos(gam)*cos(gam))
      gg=sqrt(gaux+1.)
      g1=sqrt(gg+1.)
      g2=sqrt(gg-1.)
      u=(g1*sin(ang)-cunit*g2*sin(ang-gam))/sqrt(2.*vv1*(1.+
     *1./(q1*q1)))
c
c     computation of R/T coefficients
   26 continue
      call coef8(urr,vp(1),vs(1),rho(1),vp(nz),vs(nz),rho(nz),nc,
     1rmod0,rph0)
      rph0=rph0*180./pi
   31 if(ntype.eq.2)go to 33
      if(nq.eq.1)go to 32
      call rtmat2(nz,vp,vs,rho,d,ur,freq,rp,tp,rs,ts,ntype,nc,rcoef)
      go to 40
   32 call rtmq2(fref,nz,vp,vs,rho,qp,qs,d,u,freq,rp,tp,rs,ts,ntype,
     1nc,rcoef,v1,q1,ang,gam)
      go to 40
   33 if(nq.eq.1)go to 34
      call rtmat(nz,vp,vs,rho,d,ur,freq,rpp,rps,rss,rsp,tpp,tps,tss,
     1tsp,nc,rcoef)
      go to 40
   34 call rtmatq(fref,nz,vp,vs,rho,qp,qs,d,u,freq,rpp,rps,rss,rsp,tpp,
     1tps,tss,tsp,nc,rcoef,v1,q1,ang,gam)
c
   40 b=real(rcoef)
      c=-aimag(rcoef)
      rmod=sqrt(b*b+c*c)
      if(rmod.lt.0.00001)go to 41
      rphase=atan2(c,b)*180./pi
      go to 42
   41 rphase=999.
c
  100 format(8i5)
  101 format(8f10.3)
   42 return
      end
C
C*********************************************************************
c
      SUBROUTINE RTMAT2(N,A,B,RHO,D,U,FREQ,RP,TP,RS,TS,NTYPE,NC,
     1                  RCOEF)
C
C     R/T COEFFICIENTS AT A STACK OF LAYERS FOR  PRESSURE WAVES IN
C     LIQUIDS AND FOR SH WAVES IN NON-DISSIPATIVE MEDIA
C
C     NTYPE=0....PRESSURE WAVES, LIQUID
C     NTYPE=1....SH WAVES
C
C
      DIMENSION A(N),B(N),RHO(N),D(N)
      COMPLEX RPD,RPU,TPD,TPU,RP,TP,RSD,RSU,TSD,TSU,RS,TS,
     1        A1,A2,B1,B2,AA,EP,ES,DD,FF,RCOEF
C
C
      UQ=U*U
      W=6.283185*FREQ
      M=N-1
      X=-W*D(M)
      IF(NTYPE.EQ.1)GO TO 5
C
C     PRESSURE WAVES, LIQUID
C
      AQ1=1./(A(M)*A(M))
      RHO1=RHO(M)
      AQ2=1./(A(N)*A(N))
      RHO2=RHO(N)
      V=SQRT(ABS(AQ1-UQ))
      ARG=X*V
      A1=CMPLX(V,0.)
      if(uq.gt.aq1)a1=cmplx(0.,-v)
      IF(UQ.LE.AQ1)EP=CMPLX(COS(ARG),SIN(ARG))
      IF(UQ.GT.AQ1)EP=CMPLX(EXP(ARG),0.)
      V=SQRT(ABS(AQ2-UQ))
      A2=CMPLX(V,0.)
      IF(UQ.GT.AQ2)A2=CMPLX(0.,-V)
      AA=RHO2*A1+RHO1*A2
      RPD=(RHO2*A1-RHO1*A2)/AA
      RPU=-RPD
      TPD=2.*RHO2*A1/AA
      TPU=2.*RHO1*A2/AA
      RP=RPD*EP*EP
      TP=TPD*EP
      GO TO 10
C
C     SH WAVES
C
   5  BQ1=1./(B(M)*B(M))
      Z1=RHO(M)*B(M)*B(M)
      BQ2=1./(B(N)*B(N))
      Z2=RHO(N)*B(N)*B(N)
      V=SQRT(ABS(BQ1-UQ))
      ARG=X*V
      B1=CMPLX(V,0.)
      IF(UQ.LE.BQ1)ES=CMPLX(COS(ARG),SIN(ARG))
      IF(UQ.GT.BQ1)B1=CMPLX(0.,-V)
      IF(UQ.GT.BQ1)ES=CMPLX(EXP(ARG),0.)
      V=SQRT(ABS(BQ2-UQ))
      B2=CMPLX(V,0.)
      IF(UQ.GT.BQ2)B2=CMPLX(0.,-V)
      AA=Z1*B1+Z2*B2
      RSD=(Z1*B1-Z2*B2)/AA
      RSU=-RSD
      TSD=2.*Z1*B1/AA
      TSU=2.*Z2*B2/AA
      RS=RSD*ES*ES
      TS=TSD*ES
C
C     MATRIX MULTIPLICATION FOR LAYERED MEDIUM
C
  10  IF(N.EQ.2)GO TO 1001
      II=N-2
      DO 1000 I=II,1,-1
      M=I+1
      X=-W*D(I)
      IF(NTYPE.EQ.1)GO TO 15
      AQ1=1./(A(I)*A(I))
      RHO1=RHO(I)
      AQ2=1./(A(M)*A(M))
      RHO2=RHO(M)
      V=SQRT(ABS(AQ1-UQ))
      ARG=X*V
      IF(UQ.LE.AQ1)EP=CMPLX(COS(ARG),SIN(ARG))
      IF(UQ.GT.AQ1)EP=CMPLX(EXP(ARG),0.)
      A1=CMPLX(V,0.)
      IF(UQ.GT.AQ1)A1=CMPLX(0.,-V)
      V=SQRT(ABS(AQ2-UQ))
      A2=CMPLX(V,0.)
      IF(UQ.GT.AQ2)A2=CMPLX(0.,-V)
      AA=RHO2*A1+RHO1*A2
      RPD=(RHO2*A1-RHO1*A2)/AA
      RPU=-RPD
      TPD=2.*RHO2*A1/AA
      TPU=2.*RHO1*A2/AA
      DD=1.-RPU*RP
      RP=(RPD+TPU*TPD*RP/DD)*EP*EP
      IF(I.EQ.1)FF=TPD/DD
      IF(I.GT.1)FF=TPD*EP/DD
      TP=FF*TP
      GO TO 1000
  15  BQ1=1./(B(I)*B(I))
      Z1=RHO(I)*B(I)*B(I)
      BQ2=1./(B(M)*B(M))
      Z2=RHO(M)*B(M)*B(M)
      V=SQRT(ABS(BQ1-UQ))
      ARG=X*V
      IF(UQ.LE.BQ1)ES=CMPLX(COS(ARG),SIN(ARG))
      IF(UQ.GT.BQ1)ES=CMPLX(EXP(ARG),0.)
      B1=CMPLX(V,0.)
      IF(UQ.GT.BQ1)B1=CMPLX(0.,-V)
      V=SQRT(ABS(BQ2-UQ))
      B2=CMPLX(V,0.)
      IF(UQ.GT.BQ2)B2=CMPLX(0.,-V)
      AA=Z1*B1+Z2*B2
      RSD=(Z1*B1-Z2*B2)/AA
      RSU=-RSD
      TSD=2.*Z1*B1/AA
      TSU=2.*Z2*B2/AA
      DD=1.-RSU*RS
      RS=(RSD+TSU*TSD*RS/DD)*ES*ES
      IF(I.EQ.1)FF=TSD/DD
      IF(I.GT.1)FF=TSD*ES/DD
      TS=FF*TS
C
 1000 CONTINUE
 1001 CONTINUE
      if(nc.eq.9)rcoef=rs
      if(nc.eq.10)rcoef=ts
      if(nc.eq.13)rcoef=rp
      if(nc.eq.14)rcoef=tp
C
      RETURN
      END
c
c********************************************************************
c
      SUBROUTINE rtmq2(FR,N,A,B,RHO,QA,QB,D,U,FREQ,RP,TP,RS,TS,NTYPE,
     1                   nc,rcoef,v1,q1,aincr,gamma)
C
C     R/T COEFFICIENTS AT A TRANSITION LAYER FOR PRESSURE WAVES
C     IN LIQUIDS AND FOR SH WAVES IN DISSIPATIVE MEDIA
C
C     NTYPE=0....PRESSURE WAVES, LIQUID
C     NTYPE=1....SH WAVES
C
      DIMENSION A(N),B(N),RHO(N),D(N),QA(N),QB(N)
      COMPLEX RPD,RPU,TPD,TPU,RP,TP,RSD,RSU,TSD,TSU,RS,TS,EIN,
     1        A1,A2,B1,B2,AA,EP,ES,Z1,Z2,CW,CV,aq1,aq2,bq1,bq2,rcoef,
     2        DD,FF,U,UQ
C
      UQ=U*U
      WR=6.283185*FR
      EIN=CMPLX(0.,1.)
      W=6.283185*FREQ
      CW=W
      CV=ALOG(W/WR)/3.141593+0.5*EIN
      M=N-1
      IF(NTYPE.EQ.1)GO TO 5
C
C     PRESSURE WAVES, DISSIPATIVE LIQUIDS
C
      RHO1=RHO(M)
      AQ1=A(M)*(1.+CV/QA(M))
      AQ1=1./(AQ1*AQ1)
      A1=CSQRT(AQ1-UQ)
      RHO2=RHO(N)
      AQ2=A(N)*(1.+CV/QA(N))
      AQ2=1./(AQ2*AQ2)
      A2=CSQRT(AQ2-UQ)
c
      if(n.gt.2)go to 1
      cvr=alog(w/wr)/3.14159
      va1=a(m)*(1.+cvr/qa(m))
      qa1=qa(m)
      va2=a(n)*(1.+cvr/qa(n))
      qa2=qa(n)
        call critan(v1,va1,q1,qa1,gamma,ncrit,acrit1,acrit2)
        if(ncrit.ne.0)then
          if(aincr.gt.acrit1.and.aincr.lt.acrit2)a1=-a1
        endif
        call critan(v1,va2,q1,qa2,gamma,ncrit,acrit1,acrit2)
        if(ncrit.ne.0)then
          if(aincr.gt.acrit1.and.aincr.lt.acrit2)a2=-a2
        endif
    1 continue
c
      AA=RHO2*A1+RHO1*A2
      RPD=(RHO2*A1-RHO1*A2)/AA
      RPU=-RPD
      TPD=2.*RHO2*A1/AA
      TPU=2.*RHO1*A2/AA
      EP=CEXP(-D(M)*CW*EIN*A1)
      RP=RPD*EP*EP
      TP=TPD*EP
      GO TO 10
C
C     SH WAVES, DISSIPATIVE MEDIA
C
   5  continue
      BQ1=B(M)*(1.+CV/QB(M))
      BQ1=1./(BQ1*BQ1)
      B1=CSQRT(BQ1-UQ)
      Z1=RHO(M)/BQ1
      BQ2=B(N)*(1.+CV/QB(N))
      BQ2=1./(BQ2*BQ2)
      B2=CSQRT(BQ2-UQ)
      Z2=RHO(N)/BQ2
c
      if (n.gt.2)go to 2
      cvr=alog(w/wr)/3.141593
      vb1=b(m)*(1.+cvr/qb(m))
      qb1=qb(m)
      vb2=b(n)*(1.+cvr/qb(n))
      qb2=qb(n)
        call critan(v1,vb1,q1,qb1,gamma,ncrit,acrit1,acrit2)
        if(ncrit.ne.0)then
          if(aincr.gt.acrit1.and.aincr.lt.acrit2)b1=-b1
        endif
        call critan(v1,vb2,q1,qb2,gamma,ncrit,acrit1,acrit2)
        if(ncrit.ne.0)then
          if(aincr.gt.acrit1.and.aincr.lt.acrit2)b2=-b2
        endif
    2 continue
c
      AA=Z1*B1+Z2*B2
      RSD=(Z1*B1-Z2*B2)/AA
      RSU=-RSD
      TSD=2.*Z1*B1/AA
      TSU=2.*Z2*B2/AA
      ES=CEXP(-D(M)*CW*EIN*B1)
      RS=RSD*ES*ES
      TS=TSD*ES
C
C     MATRIX MULTIPLICATION FOR LAYERED MEDIUM
C
  10  IF(N.EQ.2)GO TO 1001
      II=N-2
      DO 1000 I=II,1,-1
      M=I+1
      IF(NTYPE.EQ.1)GO TO 15
      AQ1=A(I)*(1.+CV/QA(I))
      AQ1=1./(AQ1*AQ1)
      RHO1=RHO(I)
      AQ2=A(M)*(1.+CV/QA(M))
      AQ2=1./(AQ2*AQ2)
      RHO2=RHO(M)
      A1=CSQRT(AQ1-UQ)
      A2=CSQRT(AQ2-UQ)
      AA=RHO2*A1+RHO1*A2
      RPD=(RHO2*A1-RHO1*A2)/AA
      RPU=-RPD
      TPD=2.*RHO2*A1/AA
      TPU=2.*RHO1*A2/AA
      EP=CEXP(-D(I)*CW*EIN*A1)
      DD=1.-RPU*RP
      RP=(RPD+TPU*TPD*RP/DD)*EP*EP
      IF(I.EQ.1)FF=TPD/DD
      IF(I.GT.1)FF=TPD*EP/DD
      TP=FF*TP
      GO TO 1000
  15  BQ1=B(I)*(1.+CV/QB(I))
      BQ1=1./(BQ1*BQ1)
      Z1=RHO(I)/BQ1
      BQ2=B(M)*(1.+CV/QB(M))
      BQ2=1./(BQ2*BQ2)
      Z2=RHO(M)/BQ2
      B1=CSQRT(BQ1-UQ)
      B2=CSQRT(BQ2-UQ)
      AA=Z1*B1+Z2*B2
      RSD=(Z1*B1-Z2*B2)/AA
      RSU=-RSD
      TSD=2.*Z1*B1/AA
      TSU=2.*Z2*B2/AA
      ES=CEXP(-D(I)*CW*EIN*B1)
      DD=1.-RSU*RS
      RS=(RSD+TSU*TSD*RS/DD)
      IF(I.EQ.1)FF=TSD/DD
      IF(I.GT.1)FF=TSD*ES/DD
      TS=FF*TS
C
 1000 CONTINUE
 1001 CONTINUE
      if(nc.eq.9)rcoef=rs
      if(nc.eq.10)rcoef=ts
      if(nc.eq.13)rcoef=rp
      if(nc.eq.14)rcoef=tp
C
      RETURN
      END
c
c*********************************************************************
c
        SUBROUTINE RTMAT(N,A,B,RHO,D,U,FRQ,RPP,RPS,RSS,RSP,
     *  TPP,TPS,TSS,TSP,NC,RCOEF)
C
C       P/SV REFLECTION AND TRANSMISSION COEFFICIENTS FOR A
C       NON-DISSIPATIVE TRANSITION LAYER, USING RECURSIVE FORMALISM
C
C       N = NUMBER OF LAYERS (HALFSPACES INCLUDED)
C       A,B,RHO,D = MODEL
C       U = TANGENTIAL SLOWNESS
C       FRQ = FREQUENCY
C       RCOEF = COMPLEX-VALUED R/T COEFFICIENT OF THE TYPE NC
C
        DIMENSION A(1),B(1),RHO(1),D(1)
        COMPLEX RPP,RPS,RSS,RSP,TPP,TPS,TSS,TSP,MU11,MU12,MU21,MU22,
     *  RPPD,RSPD,RPSD,RSSD,TPPD,TSPD,TPSD,TSSD,RPPU,RSPU,RPSU,RSSU,
     *  TPPU,TSPU,TPSU,TSSU,F11,F12,F21,F22,G11,G12,G21,G22,H11,H12,
     *  H21,H22,I11,I12,I21,I22,EP,ES,EX,rcoef
        M=N-1
        AQ1=1./(A(M)*A(M))
        BQ1=1./(B(M)*B(M))
        RHO1=RHO(M)
        AQ2=1./(A(N)*A(N))
        BQ2=1./(B(N)*B(N))
        RHO2=RHO(N)
        C=2.*(RHO1/BQ1-RHO2/BQ2)
        CALL RT(AQ1,BQ1,RHO1,AQ2,BQ2,RHO2,C,U,MU11,MU12,MU21,MU22,
     *  TPPD,TSPD,TPSD,TSSD,RPPU,RSPU,RPSU,RSSU,TPPU,TSPU,TPSU,TSSU)
        W=6.283185*FRQ
        UQ=U*U
        CALL PHASE(UQ,W,AQ1,BQ1,D(M),EP,ES)
        EX=EP*ES
        RPP=MU11*EP*EP
        RSP=MU12*EX
        RPS=MU21*EX
        RSS=MU22*ES*ES
        TPP=TPPD*EP
        TSP=TSPD*ES
        TPS=TPSD*EP
        TSS=TSSD*ES
        IF(N.EQ.2)  go to 1001
        II=N-2
        DO  1000  I=II,1,-1
        M=I+1
        AQ1=1./(A(I)*A(I))
        BQ1=1./(B(I)*B(I))
        RHO1=RHO(I)
        AQ2=1./(A(M)*A(M))
        BQ2=1./(B(M)*B(M))
        RHO2=RHO(M)
        C=2.*(RHO1/BQ1-RHO2/BQ2)
        CALL RT(AQ1,BQ1,RHO1,AQ2,BQ2,RHO2,C,U,RPPD,RSPD,RPSD,RSSD,
     *  TPPD,TSPD,TPSD,TSSD,RPPU,RSPU,RPSU,RSSU,TPPU,TSPU,TPSU,TSSU)
        CALL MULMAT(RPP,RSP,RPS,RSS,TPPD,TSPD,TPSD,TSSD,G11,G12,G21,
     *  G22)
        CALL MULMAT(RPP,RSP,RPS,RSS,RPPU,RSPU,RPSU,RSSU,H11,H12,H21,
     *  H22)
        H11=1.-H11
        H12=-H12
        H21=-H21
        H22=1.-H22
        CALL MATINV(H11,H12,H21,H22,I11,I12,I21,I22)
        CALL MULMAT(I11,I12,I21,I22,G11,G12,G21,G22,H11,H12,H21,H22)
        CALL MULMAT(TPPU,TSPU,TPSU,TSSU,H11,H12,H21,H22,G11,G12,G21,
     *  G22)
        MU11=RPPD+G11
        MU12=RSPD+G12
        MU21=RPSD+G21
        MU22=RSSD+G22
        CALL MULMAT(RPPU,RSPU,RPSU,RSSU,RPP,RSP,RPS,RSS,G11,G12,G21,
     *  G22)
        G11=1.-G11
        G12=-G12
        G21=-G21
        G22=1.-G22
        CALL MATINV(G11,G12,G21,G22,H11,H12,H21,H22)
        CALL MULMAT(H11,H12,H21,H22,TPPD,TSPD,TPSD,TSSD,G11,G12,G21,
     *  G22)
        CALL PHASE(UQ,W,AQ1,BQ1,D(I),EP,ES)
        F11=G11*EP
        F12=G12*ES
        F21=G21*EP
        F22=G22*ES
        CALL MULMAT(TPP,TSP,TPS,TSS,F11,F12,F21,F22,G11,G12,G21,G22)
        TPP=G11
        TSP=G12
        TPS=G21
        TSS=G22
        EX=EP*ES
        RPP=MU11*EP*EP
        RSP=MU12*EX
        RPS=MU21*EX
1000    RSS=MU22*ES*ES
1001    CONTINUE
        if(nc.eq.1)rcoef=rpp
        if(nc.eq.2)rcoef=-a(1)*rps/b(1)
        if(nc.eq.3)rcoef=a(1)*tpp/a(n)
        if(nc.eq.4)rcoef=-a(1)*tps/b(n)
        if(nc.eq.5)rcoef=-b(1)*rsp/a(1)
        if(nc.eq.6)rcoef=rss
        if(nc.eq.7)rcoef=-b(1)*tsp/a(n)
        if(nc.eq.8)rcoef=b(1)*tss/b(n)
        return
        END
C
C******************************************************************
C
        SUBROUTINE RT(AQ1,BQ1,RHO1,AQ2,BQ2,RHO2,C,U,RPPD,RSPD,RPSD,
     *  RSSD,TPPD,TSPD,TPSD,TSSD,RPPU,RSPU,RPSU,RSSU,TPPU,TSPU,
     *  TPSU,TSSU)
C
C       A ROUTINE TO RTMAT.
C       MATRICES OF R/T COEFFICIENTS AT A SINGLE INTERFACE BETWEEN
C       TWO HALFSPACES
C
        COMPLEX RPPD,RSPD,RPSD,RSSD,TPPD,TSPD,TPSD,TSSD,RPPU,RSPU,
     *  RPSU,RSSU,TPPU,TSPU,TPSU,TSSU,A1,B1,A2,B2,D1D,D2D,D1U,D2U,
     *  D,A1B1,A2B2,A1B2,A2B1,T1,T2,T3,T4,T5,T6,T7,T8,T9,T10
        UQ=U*U
        V=SQRT(ABS(AQ1-UQ))
        A1=CMPLX(V,0.)
        IF(UQ.GT.AQ1)  A1=CMPLX(0.,-V)
        V=SQRT(ABS(BQ1-UQ))
        B1=CMPLX(V,0.)
        IF(UQ.GT.BQ1)  B1=CMPLX(0.,-V)
        V=SQRT(ABS(AQ2-UQ))
        A2=CMPLX(V,0.)
        IF(UQ.GT.AQ2)  A2=CMPLX(0.,-V)
        V=SQRT(ABS(BQ2-UQ))
        B2=CMPLX(V,0.)
        IF(UQ.GT.BQ2)  B2=CMPLX(0.,-V)
        C0=C*UQ
        C1=C0-RHO1
        C2=C0+RHO2
        C3=C1+RHO2
        A1B1=A1*B1
        A2B2=A2*B2
        A1B2=A1*B2
        A2B1=A2*B1
        RHO12=RHO1*RHO2
        T1=C1*C1*A2B2+RHO12*A2B1
        T2=C2*C2*A1B1+RHO12*A1B2
        T3=C3*C3*UQ
        T4=C*C0*A1B1*A2B2
        D1D=T3+T1
        D2D=T4+T2
        D=D1D+D2D
        D1U=T3+T2
        D2U=T4+T1
        T5=2./D
        T6=A1*T5
        T7=B1*T5
        T8=A2*T5
        T9=B2*T5
        RPPD=(D2D-D1D)/D
        RPPU=(D2U-D1U)/D
        T10=U*(C3*C2+C*C1*A2B2)
        RPSD=-T6*T10
        RSPD=T7*T10
        T10=RHO12*(A1B2-A2B1)*T5
        RSSD=RPPD-T10
        RSSU=RPPU+T10
        T10=U*(C3*C1+C*C2*A1B1)
        RPSU=T8*T10
        RSPU=-T9*T10
        T10=C2*B1-C1*B2
        TPPD=RHO1*T6*T10
        TPPU=RHO2*T8*T10
        T10=U*(C3+C*A2B1)
        TPSD=-RHO1*T6*T10
        TSPU=RHO2*T9*T10
        T10=C2*A1-C1*A2
        TSSD=RHO1*T7*T10
        TSSU=RHO2*T9*T10
        T10=U*(C3+C*A1B2)
        TSPD=RHO1*T7*T10
        TPSU=-RHO2*T8*T10
        RETURN
        END
C
C******************************************************************
C
        SUBROUTINE MULMAT(A11,A12,A21,A22,B11,B12,B21,B22,C11,C12,
     *  C21,C22)
C
C       A ROUTINE TO RTMAT AND RTMATQ.
C       MATRIX MULTIPLICATION
C
        COMPLEX A11,A12,A21,A22,B11,B12,B21,B22,C11,C12,C21,C22
        C11=A11*B11+A12*B21
        C12=A11*B12+A12*B22
        C21=A21*B11+A22*B21
        C22=A21*B12+A22*B22
        RETURN
        END
C
C******************************************************************
C
        SUBROUTINE MATINV(A11,A12,A21,A22,B11,B12,B21,B22)
C
C       A ROUTINE TO RTMAT AND RTMATQ.
C       MATRIX INVERSION
C
        COMPLEX A11,A12,A21,A22,B11,B12,B21,B22,D
        D=1./(A11*A22-A12*A21)
        B11=A22*D
        B12=-A12*D
        B21=-A21*D
        B22=A11*D
        RETURN
        END
C
C******************************************************************
C
       SUBROUTINE PHASE(UQ,W,AQ,BQ,D,EP,ES)
C
C      A ROUTINE TO RTMAT AND RTMATQ
C
        COMPLEX EP,ES
        X=-W*D
        V=SQRT(ABS(AQ-UQ))
        A=X*V
        IF(UQ.GT.AQ)  GO TO 100
        EP=CMPLX(COS(A),SIN(A))
        GO TO 200
100     EP=CMPLX(EXP(A),0.)
200     V=SQRT(ABS(BQ-UQ))
        A=X*V
        IF(UQ.GT.BQ)  GO TO 300
        ES=CMPLX(COS(A),SIN(A))
        GO TO 400
300     ES=CMPLX(EXP(A),0.)
400     RETURN
        END
C
C     ***********************************************************
C
      SUBROUTINE COEF8(P,VP1,VS1,RO1,VP2,VS2,RO2,NCODE,RMOD,RPH)
C
C     ROUTINE COEF8 IS DESIGNED FOR THE COMPUTATION OF REFLECTION
C     AND TRANSMISSION COEFFICIENTS OF P, SV AND SH PLANE WAVES
C     AT A SINGLE PLANE INTERFACE BETWEEN TWO HOMOGENEOUS
C     NON-DISSIPATIVE SOLID HALFSPACES OR AT A FREE SURFACE OF A
C     NON-DISSIPATIVE HOMOGENEOUS SOLID HALFSPACE. EXPLICIT EQUATIONS
C     ARE USED. ALSO PRESSURE WAVES IN LIQUIDS ARE CONSIDERED.
C
      COMPLEX B(4),RR,C1,C2,C3,C4,H1,H2,H3,H4,H5,H6,H,HB,HC
      DIMENSION PRMT(4),D(4),DD(4)
C
      AUX=999.*3.14159/180.
      IF(NCODE.GE.9)GO TO 300
      IF(RO2.LT.0.0001)GO TO 150
      PRMT(1)=VP1
      PRMT(2)=VS1
      PRMT(3)=VP2
      PRMT(4)=VS2
      A1=VP1*VS1
      A2=VP2*VS2
      A3=VP1*RO1
      A4=VP2*RO2
      A5=VS1*RO1
      A6=VS2*RO2
      Q=2.*(A6*VS2-A5*VS1)
      PP=P*P
      QP=Q*PP
      X=RO2-QP
      Y=RO1+QP
      Z=RO2-RO1-QP
      G1=A1*A2*PP*Z*Z
      G2=A2*X*X
      G3=A1*Y*Y
      G4=A4*A5
      G5=A3*A6
      G6=Q*Q*PP
      DO 21 I=1,4
      DD(I)=P*PRMT(I)
   21 D(I)=SQRT(ABS(1.-DD(I)*DD(I)))
      IF(DD(1).LE.1..AND.DD(2).LE.1..AND.DD(3).LE.1..AND.DD(4).LE.
     11.)GO TO 100
C
C     COMPLEX COEFFICIENTS
      DO 22 I=1,4
      IF(DD(I).GT.1.)GO TO 23
      B(I)=CMPLX(D(I),0.)
      GO TO 22
   23 B(I)= CMPLX(0.,D(I))
   22 CONTINUE
      C1=B(1)*B(2)
      C2=B(3)*B(4)
      C3=B(1)*B(4)
      C4=B(2)*B(3)
      H1=G1
      H2=G2*C1
      H3=G3*C2
      H4=G4*C3
      H5=G5*C4
      H6=G6*C1*C2
      H=1./(H1+H2+H3+H4+H5+H6)
      HB=2.*H
      HC=HB*P
      GO TO (1,2,3,4,5,6,7,8),NCODE
    1 RR=H*(H2+H4+H6-H1-H3-H5)
      GO TO 26
    2 RR=VP1*B(1)*HC*(Q*Y*C2+A2*X*Z)
      GO TO 26
    3 RR=A3*B(1)*HB*(VS2*B(2)*X+VS1*B(4)*Y)
      GO TO 26
    4 RR=-A3*B(1)*HC*(Q*C4-VS1*VP2*Z)
      GO TO 26
    5 RR=-VS1*B(2)*HC*(Q*Y*C2+A2*X*Z)
      GO TO 26
    6 RR=H*(H2+H5+H6-H1-H3-H4)
      GO TO 26
    7 RR=A5*B(2)*HC*(Q*C3-VP1*VS2*Z)
      GO TO 26
    8 RR=A5*B(2)*HB*(VP1*B(3)*Y+VP2*B(1)*X)
      GO TO 26
C     REAL COEFFICIENTS
  100 E1=D(1)*D(2)
      E2=D(3)*D(4)
      E3=D(1)*D(4)
      E4=D(2)*D(3)
      S1=G1
      S2=G2*E1
      S3=G3*E2
      S4=G4*E3
      S5=G5*E4
      S6=G6*E1*E2
      S=1./(S1+S2+S3+S4+S5+S6)
      SB=2.*S
      SC=SB*P
      GO TO (101,102,103,104,105,106,107,108),NCODE
  101 R=S*(S2+S4+S6-S1-S3-S5)
      GO TO 250
  102 R=VP1*D(1)*SC*(Q*Y*E2+A2*X*Z)
      GO TO 250
  103 R=A3*D(1)*SB*(VS2*D(2)*X+VS1*D(4)*Y)
      GO TO 250
  104 R=-A3*D(1)*SC*(Q*E4-VS1*VP2*Z)
      GO TO 250
  105 R=-VS1*D(2)*SC*(Q*Y*E2+A2*X*Z)
      GO TO 250
  106 R=S*(S2+S5+S6-S1-S3-S4)
      GO TO 250
  107 R=A5*D(2)*SC*(Q*E3-VP1*VS2*Z)
      GO TO 250
  108 R=A5*D(2)*SB*(VP1*D(3)*Y+VP2*D(1)*X)
      GO TO 250
C
C    EARTHS SURFACE,COMPLEX COEFFICIENTS AND CONVERSION COEFFICIENTS
  150 A1=VS1*P
      A2=A1*A1
      A3=2.*A2
      A4=2.*A1
      A5=A4+A4
      A6=1.-A3
      A7=2.*A6
      A8=2.*A3*VS1/VP1
      A9=A6*A6
      DD(1)=P*VP1
      DD(2)=P*VS1
      DO 151 I=1,2
  151 D(I)=SQRT(ABS(1.-DD(I)*DD(I)))
      IF(DD(1).LE.1..AND.DD(2).LE.1.)GO TO 200
      DO 154 I=1,2
      IF(DD(I).GT.1.)GO TO 155
      B(I)=CMPLX(D(I),0.)
      GO TO 154
  155 B(I)= CMPLX(0.,D(I))
  154 CONTINUE
      H1=B(1)*B(2)
      H2=H1*A8
      H=1./(A9+H2)
      GO TO (161,162,166,165,163,164,168,167),NCODE
  161 RR=(-A9+H2)*H
      GO TO 26
  162 RR=-A5*B(1)*H*A6
      GO TO 26
  163 RR=A5*B(2)*H*A6*VS1/VP1
      GO TO 26
  164 RR=-(A9-H2)*H
      GO TO 26
  165 RR=A5*H1*H
      GO TO 26
  166 RR=-A7*B(1)*H
      GO TO 26
  167 RR=A7*B(2)*H
      GO TO 26
  168 RR=A5*VS1*H1*H/VP1
C
  26  Z2=REAL(RR)
      Z3=AIMAG(RR)
      IF(Z2.EQ.0..AND.Z3.EQ.0.)GO TO 157
      RMOD=SQRT(Z2*Z2+Z3*Z3)
      RPH=ATAN2(Z3,Z2)
      IF(RMOD.LT.0.00001)RPH=AUX
      RETURN
  157 RMOD=0.
      RPH=AUX
      RETURN
C
C     EARTHS SURFACE,REAL COEFFICIENTS AND CONVERSION COEFFICIENTS
  200 S1=D(1)*D(2)
      S2=A8*S1
      S=1./(A9+S2)
      GO TO (201,202,206,205,203,204,208,207),NCODE
  201 R=(-A9+S2)*S
      GO TO 250
  202 R=-A5*D(1)*S*A6
      GO TO 250
  203 R=A5*D(2)*S*A6*VS1/VP1
      GO TO 250
  204 R=(S2-A9)*S
      GO TO 250
  205 R=A5*S1*S
      GO TO 250
  206 R=A7*D(1)*S
      GO TO 250
  207 R=A7*D(2)*S
      GO TO 250
  208 R=-A5*VS1*S1*S/VP1
  250 IF(R.LT.0.)GO TO 251
      RMOD=R
      RPH=0.
      IF(RMOD.LT.0.00001)RPH=AUX
      RR=CMPLX(R,0.)
      RETURN
  251 RMOD=-R
      RPH=-3.14159
      IF(RMOD.LT.0.00001)RPH=AUX
      RR=CMPLX(R,0.)
      RETURN
C
C     SH COEFFICIENTS OF REFLECTION, TRANSMISSION AND CONVERSION
C
  300 IF(NCODE.GE.13)GO TO 400
      IF(RO2.LT.0.0001) GO TO 302
      A1=P*VS1
      A2=P*VS2
      A3=RO1*VS1
      A4=RO2*VS2
      A5=SQRT(ABS(1.-A1*A1))
      A6=SQRT(ABS(1.-A2*A2))
      C1=A5
      IF(A2.LE.1.)C2=A6
      IF(A2.GT.1.)C2=CMPLX(0.,A6)
      C1=C1*A3
      C2=C2*A4
      H=1./(C1+C2)
      IF(NCODE.EQ.10)GO TO 301
      RR=H*(C1-C2)
      GO TO 26
  301 RR=2.*C1*H
      GO TO 26
  302 RMOD=1.
      RPH=0.
      IF(NCODE.EQ.10)GO TO 303
      RETURN
  303 RMOD=2.
      RPH=0.
      RETURN
C
C     PRESSURE R/T COEFFICIENTS, LIQUIDS
C
  400 IF(RO2.LT.0.0001) GO TO 402
      A1=P*VP1
      A2=P*VP2
      A3=RO1*VP1
      A4=RO2*VP2
      A5=SQRT(ABS(1.-A1*A1))
      A6=SQRT(ABS(1.-A2*A2))
      C1=A5
      IF(A2.LT.1.)C2=A6
      IF(A2.GE.1.)C2=CMPLX(0.,A6)
      C1=A4*C1
      C2=A3*C2
      H=1./(C1+C2)
      IF(NCODE.EQ.14)GO TO 401
      RR=H*(C1-C2)
      GO TO 26
  401 RR=2.*C1*H
      GO TO 26
  402 RMOD=1.
      RPH=3.141593
      IF(NCODE.EQ.14)GO TO 403
      RETURN
  403 RMOD=0.
      RPH=AUX
      RETURN
      END
C
C     *************************************************************
C
        SUBROUTINE RTMATQ(FR,N,A,B,RHO,QA,QB,D,U,FRQ,
     *  RPP,RPS,RSS,RSP,TPP,TPS,TSS,TSP,NC,RCOEF,v1,q1,aincr,gamma)
C
C       P/SV REFLECTION AND TRANSMISSION COEFFICIENTS FOR A DISSIPATIVE
C       TRANSITION LAYER, USING RECURSIVE ALGORITHM.
C
C       FR = REFERENCE FREQUENCY
C       N = NUMBER OF LAYERS (HALFSPACES INCLUDED)
C       A,B,RHO,QA,QB,D = MODEL (QA AND QB INDEPENDENT OF FREQUENCY)
C       U = TANGENTIAL SLOWNESS (COMPLEX-VALUED)
C       FRQ = FREQUENCY
C       RCOEF = COMPLEX-VALUED R/T COEFFICIENT OF THE TYPE NC
C
C       LAYERS 1 AND N ARE HALFSPACES
C
C
        DIMENSION A(1),B(1),RHO(1),D(1),QA(1),QB(1),nc(1)
        COMPLEX RPP,RPS,RSS,RSP,TPP,TPS,TSS,TSP,MU11,MU12,MU21,MU22,
     *  RPPD,RSPD,RPSD,RSSD,TPPD,TSPD,TPSD,TSSD,RPPU,RSPU,RPSU,RSSU,
     *  TPPU,TSPU,TPSU,TSSU,F11,F12,F21,F22,G11,G12,G21,G22,H11,H12,
     *  H21,H22,I11,I12,I21,I22,EP,ES,EX,AQ1,BQ1,AQ2,BQ2,C,CW,CV,
     *  EIN,rcoef,apom,bpom,apom2,bpom2,u,uq
C
        WR=6.28319*FR
        EIN=CMPLX(0.,1.)
        W=6.28319*FRQ
        CW=W
        CV=ALOG(W/WR)/3.141593+0.5*EIN
C
        M=N-1
        RHO1=RHO(M)
        AQ1=A(M)*(1.+CV/QA(M))
        APOM=AQ1
        AQ1=1./(AQ1*AQ1)
        BQ1=B(M)*(1.+CV/QB(M))
        BPOM=BQ1
        BQ1=1./(BQ1*BQ1)
        RHO2=RHO(N)
        AQ2=A(N)*(1.+CV/QA(N))
        APOM2=AQ2
        AQ2=1./(AQ2*AQ2)
        BQ2=B(N)*(1.+CV/QB(N))
        BPOM2=BQ2
        BQ2=1./(BQ2*BQ2)
c
c       determination of possible critical angles in critan (for n=2)
        if(n.gt.2) go to 1
        cvr=alog(w/wr)/3.141593
        va1=a(m)*(1.+cvr/qa(m))
        vb1=b(m)*(1.+cvr/qb(m))
        qa1=qa(m)
        qb1=qb(m)
        va2=a(n)*(1.+cvr/qa(n))
        vb2=b(n)*(1.+cvr/qb(n))
        qa2=qa(n)
        qb2=qb(n)
    1   continue
c
        C=2.*(RHO1/BQ1-RHO2/BQ2)
        CALL RTQ(AQ1,BQ1,RHO1,AQ2,BQ2,RHO2,C,U,MU11,MU12,MU21,MU22,
     *  TPPD,TSPD,TPSD,TSSD,RPPU,RSPU,RPSU,RSSU,TPPU,TSPU,TPSU,TSSU,
     *  aincr,gamma,v1,q1,va1,vb1,va2,vb2,qa1,qb1,qa2,qb2,n)
        UQ=U*U
        EP=CEXP(-D(M)*CW*EIN*CSQRT(AQ1-UQ))
        ES=CEXP(-D(M)*CW*EIN*CSQRT(BQ1-UQ))
        EX=EP*ES
        RPP=MU11*EP*EP
        RSP=MU12*EX
        RPS=MU21*EX
        RSS=MU22*ES*ES
        TPP=TPPD*EP
        TSP=TSPD*ES
        TPS=TPSD*EP
        TSS=TSSD*ES
        IF(N.EQ.2)  go to 1001
C
C       LOOP FOR LAYERS
C
        II=N-2
        DO  1000  I=II,1,-1
        M=I+1
       AQ1=A(I)*(1.+CV/QA(I))
       apom=aq1
       AQ1=1./(AQ1*AQ1)
       BQ1=B(I)*(1.+CV/QB(I))
       bpom=bq1
       BQ1=1./(BQ1*BQ1)
       RHO1=RHO(I)
       AQ2=A(M)*(1.+CV/QA(M))
       AQ2=1./(AQ2*AQ2)
       BQ2=B(M)*(1.+CV/QB(M))
       BQ2=1./(BQ2*BQ2)
        RHO2=RHO(M)
        C=2.*(RHO1/BQ1-RHO2/BQ2)
        CALL RTQ(AQ1,BQ1,RHO1,AQ2,BQ2,RHO2,C,U,RPPD,RSPD,RPSD,RSSD,
     *  TPPD,TSPD,TPSD,TSSD,RPPU,RSPU,RPSU,RSSU,TPPU,TSPU,TPSU,TSSU,
     *  aincr,gamma,v1,q1,va1,vb1,va2,vb2,qa1,qb1,qa2,qb2,n)
        CALL MULMAT(RPP,RSP,RPS,RSS,TPPD,TSPD,TPSD,TSSD,G11,G12,G21,
     *  G22)
        CALL MULMAT(RPP,RSP,RPS,RSS,RPPU,RSPU,RPSU,RSSU,H11,H12,H21,
     *  H22)
        H11=1.-H11
        H12=-H12
        H21=-H21
        H22=1.-H22
        CALL MATINV(H11,H12,H21,H22,I11,I12,I21,I22)
        CALL MULMAT(I11,I12,I21,I22,G11,G12,G21,G22,H11,H12,H21,H22)
        CALL MULMAT(TPPU,TSPU,TPSU,TSSU,H11,H12,H21,H22,G11,G12,G21,
     *  G22)
        MU11=RPPD+G11
        MU12=RSPD+G12
        MU21=RPSD+G21
        MU22=RSSD+G22
        CALL MULMAT(RPPU,RSPU,RPSU,RSSU,RPP,RSP,RPS,RSS,G11,G12,G21,
     *  G22)
        G11=1.-G11
        G12=-G12
        G21=-G21
        G22=1.-G22
        CALL MATINV(G11,G12,G21,G22,H11,H12,H21,H22)
        CALL MULMAT(H11,H12,H21,H22,TPPD,TSPD,TPSD,TSSD,G11,G12,G21,
     *  G22)
        EP=CEXP(-D(I)*CW*EIN*CSQRT(AQ1-UQ))
        ES=CEXP(-D(I)*CW*EIN*CSQRT(BQ1-UQ))
        F11=G11*EP
        F12=G12*ES
        F21=G21*EP
        F22=G22*ES
        CALL MULMAT(TPP,TSP,TPS,TSS,F11,F12,F21,F22,G11,G12,G21,G22)
        TPP=G11
        TSP=G12
        TPS=G21
        TSS=G22
        EX=EP*ES
        RPP=MU11*EP*EP
        RSP=MU12*EX
        RPS=MU21*EX
1000    RSS=MU22*ES*ES
1001    continue
        if (nc(1).eq.1)rcoef=rpp
        if(nc(1).eq.2)rcoef=-apom*rps/bpom
        if(nc(1).eq.3)rcoef=apom*tpp/apom2
        if(nc(1).eq.4)rcoef=-apom*tps/bpom2
        if(nc(1).eq.5)rcoef=-bpom*rsp/apom
        if(nc(1).eq.6)rcoef=rss
        if(nc(1).eq.7)rcoef=-bpom*tsp/apom2
        if(nc(1).eq.8)rcoef=bpom*tss/bpom2
        RETURN
        END
C
C      ************************************************************
C
        SUBROUTINE RTQ(AQ1,BQ1,RHO1,AQ2,BQ2,RHO2,C,U,RPPD,RSPD,RPSD,
     *  RSSD,TPPD,TSPD,TPSD,TSSD,RPPU,RSPU,RPSU,RSSU,TPPU,TSPU,
     *  TPSU,TSSU,aincr,gamma,
     *  v1,q1,va1,vb1,va2,vb2,qa1,qb1,qa2,qb2,n)
C
C       A ROUTINE TO RTMATQ
C       R/T COEFFICIENTS AT A SINGLE INTERFACE BETWEEN
C       TWO HOMOGENEOUS DISSIPATIVE HALFSPACES
C
        COMPLEX RPPD,RSPD,RPSD,RSSD,TPPD,TSPD,TPSD,TSSD,RPPU,RSPU,
     *  RPSU,RSSU,TPPU,TSPU,TPSU,TSSU,A1,B1,A2,B2,D1D,D2D,D1U,D2U,
     *  D,A1B1,A2B2,A1B2,A2B1,T1,T2,T3,T4,T5,T6,T7,T8,T9,T10,
     *  AQ1,BQ1,AQ2,BQ2,C,C1,C2,C3,C0,U,UQ
C
        uq=u*u
        A1=CSQRT(AQ1-UQ)
        B1=CSQRT(BQ1-UQ)
        A2=CSQRT(AQ2-UQ)
        B2=CSQRT(BQ2-UQ)
c
c       determination of possible critical angles in critan, for n=2
        if(n.gt.2)go to 1
        call critan(v1,va1,q1,qa1,gamma,ncrit,acrit1,acrit2)
        if(ncrit.ne.0)then
          if(aincr.gt.acrit1.and.aincr.lt.acrit2)a1=-a1
        endif
        call critan(v1,vb1,q1,qb1,gamma,ncrit,acrit1,acrit2)
        if(ncrit.ne.0)then
          if(aincr.gt.acrit1.and.aincr.lt.acrit2)b1=-b1
        endif
        call critan(v1,va2,q1,qa2,gamma,ncrit,acrit1,acrit2)
        if(ncrit.ne.0)then
          if(aincr.gt.acrit1.and.aincr.lt.acrit2)a2=-a2
        endif
        call critan(v1,vb2,q1,qb2,gamma,ncrit,acrit1,acrit2)
        if(ncrit.ne.0)then
          if(aincr.gt.acrit1.and.aincr.lt.acrit2)b2=-b2
        endif
c
    1   C0=C*UQ
        C1=C0-RHO1
        C2=C0+RHO2
        C3=C1+RHO2
        A1B1=A1*B1
        A2B2=A2*B2
        A1B2=A1*B2
        A2B1=A2*B1
        RHO12=RHO1*RHO2
        T1=C1*C1*A2B2+RHO12*A2B1
        T2=C2*C2*A1B1+RHO12*A1B2
        T3=C3*C3*UQ
        T4=C*C0*A1B1*A2B2
        D1D=T3+T1
        D2D=T4+T2
        D=D1D+D2D
        D1U=T3+T2
        D2U=T4+T1
        T5=2./D
        T6=A1*T5
        T7=B1*T5
        T8=A2*T5
        T9=B2*T5
        RPPD=(D2D-D1D)/D
        RPPU=(D2U-D1U)/D
        T10=U*(C3*C2+C*C1*A2B2)
        RPSD=-T6*T10
        RSPD=T7*T10
        T10=RHO12*(A1B2-A2B1)*T5
        RSSD=RPPD-T10
        RSSU=RPPU+T10
        T10=U*(C3*C1+C*C2*A1B1)
        RPSU=T8*T10
        RSPU=-T9*T10
        T10=C2*B1-C1*B2
        TPPD=RHO1*T6*T10
        TPPU=RHO2*T8*T10
        T10=U*(C3+C*A2B1)
        TPSD=-RHO1*T6*T10
        TSPU=RHO2*T9*T10
        T10=C2*A1-C1*A2
        TSSD=RHO1*T7*T10
        TSSU=RHO2*T9*T10
        T10=U*(C3+C*A1B2)
        TSPD=RHO1*T7*T10
        TPSU=-RHO2*T8*T10
        RETURN
        END
C
C     *************************************************************
C
        COMPLEX FUNCTION CSQ(X)
C
C       A ROUTINE TO RTMATQ.
C
        COMPLEX X
        A=CABS(X)
        B=REAL(X)
        RE=0.7071*SQRT(A+B)
        AI=-0.7071*SQRT(A-B)
        CSQ=CMPLX(RE,AI)
        RETURN
        END
c
c**********************************************************************
c
      subroutine critan(v1,v2,q1,q2,gamma,n,acrit1,acrit2)
c
c     A routine to compute the discrete critical angles
c
c     auxw1=sqrt(1.+1./(q1*q1))
c     w1=2./(v1*v1*(1.+auxw1))
c     auxw2=sqrt(1.+1./(q2*q2))
c     w2=2./(v2*v2*(1.+auxw2))
      auxw1=1.+1./(q1*q1)
      auxw2=1.+1./(q2*q2)
      w1=1./(v1*v1*auxw1)
      w2=1./(v2*v2*auxw2)
      c=(2.*w2*q1)/(w1*q2)-1.
      a=(w2*q1*cos(gamma))/(w1*q2)
c
      n = 0
      acrit1=1000.
      acrit2=1000.
c
      caux=1./(c*c)
      cos2g=cos(gamma)*cos(gamma)
      if(cos2g.le.caux)then
        aux=sqrt(1.-c*c*cos2g)
        omega1=(1.+c*cos2g+sin(gamma)*aux)/2.
        omega2=(1.+c*cos2g-sin(gamma)*aux)/2.
        if(0.lt.omega1.and.omega1.lt.1)then
          acrit1=asin(sqrt(omega1))
          if(acrit1.lt.gamma)then
            acrit1=1000.
            go to 2
          endif
          if(abs(acrit1-gamma).lt.0.000001)then
            acrit1=1000.
            go to 2
          endif
          aim1 = a - sin(acrit1)*sin(acrit1-gamma)
          if(abs(aim1).gt.0.0001)then
            acrit1 = 1000.
            go to 2
          endif
c         n=1
          g=sqrt(1.+1./(q1*q1*cos2g))
          rea1=w2-(w1*((1.+g)*sin(acrit1)*sin(acrit1)-
     *    (g-1.)*sin(acrit1-gamma)*sin(acrit1-gamma)))/2.
          if(rea1.gt.0)then
            acrit1=1000.
c           n=0
          else
            n=1
          endif
        endif
c
   2    if(0.lt.omega2.and.omega2.lt.1)then
          acrit2=asin(sqrt(omega2))
          if(acrit2.lt.gamma)then
            acrit2=1000.
            go to 3
          endif
          if(abs(acrit2-gamma).lt.0.000001)then
            acrit2=1000.
            go to 3
          endif
          aim2=a - sin(acrit2)*sin(acrit2-gamma)
          if(abs(aim2).gt.0.0001)then
            acrit2=1000.
            go to 3
          endif
c         n=n+1
          g=sqrt(1.+1./(q1*q1*cos2g))
          rea2=w2-(w1*((1.+g)*sin(acrit2)*sin(acrit2)-
     *    (g-1.)*sin(acrit2-gamma)*sin(acrit2-gamma)))/2.
          if(rea2.gt.0)then
            acrit2=1000.
c           n=n-1
          else
            n=n+1
          endif
        endif
      endif
c
   3  if(acrit2.lt.acrit1)then
        acr11=acrit2
        acrit2=acrit1
        acrit1=acr11
      endif
c
      if(abs(acrit1-acrit2).lt.0.00005)then
      n=1
      acrit2=1000.
      endif
      return
      end
c
c***********************************************************************
c