      PROGRAM BOUNCE

C     Program for computing the reflection coefficient, R(theta)

C     Michael B. Porter

      INCLUDE 'comc.f'

C     *** Read in environmental info ***

      TITLE = 'BOUNCE-'

      CALL READIN( TITLE, FREQ, MAXMED, NMEDIA,
     &   TOPOPT, CPT, CST, RHOT, BUMDEN, ETA, XI, NG, SIGMA, DEPTH, 
     &   BOTOPT, CPB, CSB, RHOB,
     &   ENVFIL, PRTFIL )

C     *** Spectral limits ***

      READ(  ENVFIL, *    ) CLOW, CHIGH
      WRITE( PRTFIL, 500  ) CLOW, CHIGH
  500 FORMAT( /, ' CLOW = ', G12.5, '  CHIGH = ', G12.5 )

C     *** Maximum range for calculations ***

      READ(  ENVFIL, * ) RMAX
      WRITE( PRTFIL, * ) 'RMAX = ', RMAX

      OMEGA2 = ( 2.0 * PI * FREQ ) ** 2

C     *** Optionally read in bottom reflection coefficient data ***

      DEGRAD = PI / 180.0

      IF ( BOTOPT(1:1) .EQ. 'F' ) THEN
         WRITE( PRTFIL, * )
     &      'Using tabulated bottom reflection coefficient'
         READ( BRCFIL, * )
     &      NBPTS, ( THETAB( I ), RB( I ), PHIB( I ), I = 1, NBPTS )

         DO I = 1, NBPTS
            PHIB( I ) = DEGRAD * PHIB( I )
         END DO
      ENDIF

C     *** Optionally read in internal reflection coefficient data ***

      IF ( BOTOPT(1:1) .EQ. 'P' ) THEN
         WRITE( PRTFIL, * ) 'Reading precalculated refl. coeff. table'
         OPEN( FILE = 'IRCFIL', UNIT = IRCFIL, STATUS = 'OLD' )

         READ( IRCFIL, * ) NKPTS, ( XTAB( IK ), FTAB( IK ),
     &                 GTAB( IK ), ITAB( IK ), IK = 1, NKPTS )
         IF ( NKPTS .GT. MXREF ) THEN
            WRITE( PRTFIL, * ) 'Too many k-points in IRC file'
            WRITE( PRTFIL, * ) 'NKPTS, MXREF = ', NKPTS, MXREF
            STOP 'FATAL ERROR: Too many k-points in IRC file'
         ENDIF
      ENDIF

      CLOSE( ENVFIL )

C     *** Write internal reflection coefficient data ***

      WRITE( PRTFIL, * ) 'Writing internal refl. coeff. table'
C     --- Compute NKPTS
      IF ( CLOW .EQ. 0.0 ) THEN
         WRITE( PRTFIL, * )
     &         'FATAL ERROR: CLOW must be greater than zero'
         STOP  'FATAL ERROR: CLOW must be greater than zero'
      ENDIF

      RKMIN = SQRT( OMEGA2 ) / CHIGH
      RKMAX = SQRT( OMEGA2 ) / CLOW
      IF ( CHIGH .GE. 1.0E6 ) RKMIN = 0.0
      NKPTS = 1000.0 * RMAX * ( RKMAX - RKMIN ) / ( 2.0 * PI )
      WRITE( PRTFIL, * ) 'NKPTS = ', NKPTS

      IF ( NKPTS .GT. MXREF ) THEN
         WRITE( PRTFIL, * ) 'FATAL ERROR: Too many k-points'
         WRITE( PRTFIL, * ) 'NKPTS, MXREF = ', NKPTS, MXREF
         STOP               'FATAL ERROR: Too many k-points'
      ENDIF

      CALL COMPR( RKMIN, RKMAX )

      STOP
      END
C**********************************************************************C
      SUBROUTINE INIT

C     Initializes arrays defining difference equations

      INCLUDE 'comc.f'

      LOGICAL ELFLAG
      COMPLEX*16 CP( MAXN ), CS( MAXN ), CP2, CS2

      COMMON /SPACE/ CP, CS

      ELFLAG   = .FALSE.
      CMIN     = 1.0E6
      NFACT    = 0
      LOC( 1 ) = 0

C     *** Loop over media ***

      DO MED = 1, NMEDIA
         IF ( MED .NE. 1 ) LOC( MED ) = LOC( MED-1 ) + N( MED-1 ) + 1
         N1 = N(   MED ) + 1
         I  = LOC( MED ) + 1

         IF ( LOC( MED ) + N1 .GT. MAXN ) THEN
            WRITE( PRTFIL, * )
     &           'FATAL ERROR: Insufficient storage for mesh'
            STOP 'FATAL ERROR: Insufficient storage for mesh'
         ENDIF

C        *** PROFIL reads in the data for a medium ***

         CALL PROFIL( DEPTH, CP( I ), CS( I ), RHO( I ),
     &      MED, N1, FREQ, TOPOPT(1:1), TOPOPT(3:3), 'TAB',
     &      ENVFIL, PRTFIL )

C        *** Load diagonals of the finite-difference equations ***

         IF ( CS( I ) .EQ. ( 0.0, 0.0 ) ) THEN

C           --- Case of an acoustic medium ---

            MATER( MED ) = 'ACOUSTIC'
            IF ( NFACT .EQ. 0 ) NFACT = MED
            NLACT = MED

            DO J = I, I + N( MED )
               CMIN = MIN( DBLE( CP( J ) ), CMIN )
               B1( J ) = -2.0 +  H( MED ) ** 2 * OMEGA2 / CP( J ) ** 2
            END DO
         ELSE

C           --- Case of an elastic medium ---

            IF ( SIGMA( MED ) .NE. 0.0 ) THEN
               WRITE( PRTFIL, * ) 'Rough elastic interface not allowed'
               WRITE( PRTFIL, * ) 'PROGRAM ABORTING'
               STOP
     &         'ERROR IN KRAKENC: Rough elastic interface not allowed'
            ENDIF

            MATER( MED ) = 'ELASTIC'
            ELFLAG = .TRUE.
            TWOH   = 2.0 * H( MED )

            DO J = I, I + N( MED )
               CMIN = MIN( DBLE( CS( J ) ), CMIN )

               CP2 = CP( J ) ** 2
               CS2 = CS( J ) ** 2

               B1( J ) = TWOH / ( RHO( J ) * CS2 )
               B2( J ) = TWOH / ( RHO( J ) * CP2 )
               B3( J ) = 4.0 * TWOH * RHO(J) * CS2 * ( CP2 - CS2 ) / CP2
               B4( J ) = TWOH * ( CP2 - 2.0*CS2 ) / CP2

               RHO( J ) = TWOH * OMEGA2 * RHO( J )
            END DO

         ENDIF
      END DO   ! Next MED

C     *** Bottom properties ***

      IF ( BOTOPT(1:1) .EQ. 'A' ) THEN
         IF ( CSB .NE. ( 0.0, 0.0 ) ) THEN
C           --- Elastic bottom:
            ELFLAG = .TRUE.
            CMIN = MIN( CMIN, DBLE( CSB ) )
         ELSE
C           --- Acoustic bottom:
            CMIN = MIN( CMIN, DBLE( CPB ) )
         ENDIF
      ENDIF

C     *** Top properties ***

      IF ( TOPOPT(2:2) .EQ. 'A' ) THEN
         IF ( CST .NE. ( 0.0, 0.0 ) ) THEN
C           --- Elastic top:
            ELFLAG = .TRUE.
            CMIN = MIN( CMIN, DBLE( CST ) )
         ELSE
C           --- Acoustic top:
            CMIN = MIN( CMIN, DBLE( CPT ) )
         ENDIF
      ENDIF

      IF ( ELFLAG ) CMIN = 0.9 * CMIN
      CLOW = MAX( CLOW, 0.99 * CMIN )

      RETURN
      END
C**********************************************************************C
      SUBROUTINE COMPR( RKMIN, RKMAX )

C     Computes the reflection coefficient for k in [RKMIN, RKMAX]

      INCLUDE 'comc.f'

      COMPLEX*16 X, F, G
      CHARACTER BCTYPE*1

      DO MED = 1, NMEDIA
         N( MED ) = NG( MED )
         H( MED ) = ( DEPTH( MED + 1 ) - DEPTH( MED ) ) / N( MED )
      END DO

      HV( 1 ) = H( 1 )
      CALL INIT

      DELTAK = ( RKMAX - RKMIN ) / ( NKPTS - 1 )

      DO IK = 1, NKPTS
         RK = RKMIN + ( IK - 1 ) * DELTAK
         X = RK ** 2

C        --- Bottom impedance
         BCTYPE(1:1) = BOTOPT(1:1)
         CALL BCIMP( X, BCTYPE, 'BOT', CPB, CSB, RHOB, F, G, IPOW )

C        --- Shoot through acoustic layers
         CALL ACOUST( X, F, G, IPOW  )

         XTAB( IK ) = X
         FTAB( IK ) = F
         GTAB( IK ) = G
         ITAB( IK ) = IPOW
      END DO

      OPEN( FILE = 'IRCFIL', UNIT = IRCFIL )
      WRITE( IRCFIL, * ) '''', TITLE, '''', FREQ
      WRITE( IRCFIL, * ) NKPTS, ( XTAB( IK ),
     &           FTAB( IK ), GTAB( IK ), ITAB( IK ), IK = 1, NKPTS )

      RETURN
      END
C**********************************************************************C
      SUBROUTINE ACOUST( X, F, G, IPOW )

C     Shoot through acoustic layers

      INCLUDE 'comc.f'

      COMPLEX*16 X, F, G, P0, P1, P2, H2K2
      PARAMETER ( ROOF = 1.0E5, FLOOR = 1.0E-5, IPOWR = 5, IPOWF = -5 )

      IF ( NFACT .EQ. 0 ) RETURN 

C     *** Loop over successive acoustic media ***

      DO MED = NLACT, NFACT, -1
         H2K2 = H( MED ) ** 2 * X
         I = LOC( MED ) + N( MED ) + 1
         RHOM = RHO( I )
         P1 = -2.0 * G
         P2 = ( B1( I ) - H2K2 ) * G - 2.0 * H( MED ) * F * RHOM

C        *** Shoot through a single medium ***

         DO I = LOC( MED ) + N( MED ), LOC( MED ) + 1, -1

            P0 = P1
            P1 = P2
            P2 = ( H2K2 - B1( I ) ) * P1 - P0

   60       IF ( ABS( DBLE( P2 ) ) .GE. ROOF ) THEN
               P0 = FLOOR * P0
               P1 = FLOOR * P1
               P2 = FLOOR * P2
               IPOW = IPOW - IPOWF
               GOTO 60
            ENDIF
         END DO

C        --- F = P'/RHO and G = -P since FP+GP'/RHO = 0
         F = -( P2 - P0 ) / ( 2.0 * H( MED ) ) / RHOM
         G = -P1
      END DO

      RETURN
      END
