      PROGRAM BELLHOP

C     Beam tracing in cylindrical coordinates
C     Michael B. Porter and Homer P. Bucker

      INCLUDE 'COMBM.FOR'

      PARAMETER ( MXNSD = 51, MXNRZ = 201000, MXNBEAMS = 30001 )

C     The following can be changed to reduce storage:
C
C     MXNSD   max # of Source Depths
C     MXNRD            range-depths coords. for TL calculation
C                      (also used to limit number of arrivals)
C     MXNBEAMS         beams

      INTEGER     NARR( MXNRZ )
      REAL      DELARR( MXNRZ )
      COMPLEX     AARR( MXNRZ )

      REAL      SD( MXNSD ), ALPHA( MXNBEAMS )
      REAL*8    TEND
      COMPLEX   U( MXNRZ ), EPS, PICKEPS
      CHARACTER TITLE*80, BOTOPT*2, RUNTYP*2, BEAMTYP*2
      EQUIVALENCE ( U, AARR )

C     *** Begin by reading in control data ***

      CALL CLTIME   ! clear the time counter

      CALL READIN( TITLE, FREQ, ALPHA, NBEAMS, MXNBEAMS, ISINGL,
     &   NIMAGE, IBWIN, DELTAS, MXN, ZBOX, RBOX, EPMULT, RLOOP,
     &   SD, NSD, MXNSD, RD, NRD, MXNRD, R, NR, MXNR,
     &   RNG, DEP, TBR, TBZ, NBTYPTS, MXBTY,
     &   THETATOP, RTOP, PHITOP, NTOPPTS,
     &   THETABOT, RBOT, PHIBOT, NBOTPTS, MXREF,
     &   TOPOPT, DEPTHT, CPT, RHOT,
     &   BOTOPT, DEPTHB, CPB, RHOB, RUNTYP, BEAMTYP )

      IF ( NR * NRD .GT. MXNRZ ) THEN
         WRITE( PRTFIL, * ) 'NRZ*NRD, MXNRZ = ', NR * NRD, MXNRZ
         CALL ERROUT( PRTFIL, 'F', 'READIN', 'Too many field points'  )
      ENDIF

      MXNARR = MXNRZ / ( NRD * NR )
      WRITE( PRTFIL, * )
      WRITE( PRTFIL, * ) '( Maximum number of arrivals = ', MXNARR, ')'

      OMEGA  = 2.0 * PI * FREQ
      IBWIN2 = IBWIN **2

      IF ( NR .GE. 2 ) THEN
         DELTAR = R( NR ) - R( NR - 1 )
      ELSE
         DELTAR = 0.0
      ENDIF

C     Beam angle stuff

      DALPHA = 0.0
      IF ( NBEAMS .NE. 1 )
     &   DALPHA = ( ALPHA( NBEAMS ) - ALPHA( 1 ) ) / ( NBEAMS - 1 )

C     *** Loop over source depths ***

      DO IS = 1, NSD
         RS = 0.0
         ZS = SD( IS )

C        --- Zero out field vector
         DO I = 1, NR * NRD
            U( I ) = 0.0
         END DO

C        --- Zero out arrival data
         IF( RUNTYP(1:1) .EQ. 'A' ) THEN
           DO I = 1, NRD * NR
             NARR( I ) = 0
           ENDDO
           DO J = 1, MXNRZ
             AARR(   J ) = 0.0
             DELARR( J ) = 0.0
           END DO              
         END IF

         CALL SSP( RS, ZS, C, CR, CZ, CRR, CRZ, CZZ, TOPOPT, 'TAB' )

         RADMAX = 10 * C / FREQ  ! 10 wavelength max radius

C        --- Check that there are enough beams

         DALPHAOPT = SQRT( C / ( 6.0 * FREQ * R( NR ) ) )
         NBEAMSOPT = 2 + ( ALPHA( NBEAMS ) - ALPHA( 1 ) ) / DALPHAOPT

         IF ( RUNTYP(1:1) .EQ. 'C' .AND. NBEAMS .LT. NBEAMSOPT ) THEN
            CALL ERROUT( PRTFIL, 'W', 'BELLHOP', 'Too few beams' )
            WRITE( PRTFIL, * ) 'NBEAMS should be at least = ', NBEAMSOPT
         ENDIF

C        *** Trace successive beams ***

         DO IBEAM = 1, NBEAMS
            WRITE( *, * ) 'Tracing beam ', IBEAM,
     &         ALPHA( IBEAM ) * 180.0 / PI

C           --- Single beam run?
            IF ( ISINGL .EQ. 0 .OR. IBEAM .EQ. ISINGL ) THEN

C           *** Trace a ray ***
            CALL TRACE( DELTAS, RS, ZS, ALPHA( IBEAM ),
     &         BEAMTYP, ZBOX, RBOX, BOTOPT )

            IF ( RUNTYP(1:1) .EQ. 'R' ) THEN
C              *** Write the ray trajectory to RAYFIL ***
               CALL WRTRAY( RV, ZV, NSTEPS, DEPTHT, DEPTHB )
            ELSE
C              *** Compute the contribution to the field ***

C              --- Compute 'optimal' beam constant:
               EPS = PICKEPS( BEAMTYP(1:1), OMEGA, C, CZ,
     &                  ALPHA( IBEAM ), DALPHA, RLOOP, EPMULT )

               IF ( RUNTYP(2:2) .EQ. 'R' ) THEN
                  CALL INFLUR( U, DELTAR, EPS, ALPHA( IBEAM ), NIMAGE,
     &               IBWIN2, RUNTYP, RADMAX, BEAMTYP )
               ELSE IF ( RUNTYP(2:2) .EQ. 'C' ) THEN
                  CALL INFLUC( U, DELTAR, EPS, ALPHA( IBEAM ), NIMAGE,
     &               IBWIN2, RUNTYP, RADMAX, BEAMTYP )
               ELSE
                  CALL INFLUG( U, ZS, ALPHA( IBEAM ), RUNTYP,
     &               DALPHA, NSD, NARR, DELARR, AARR, MXNARR )
               ENDIF
            ENDIF
            ENDIF
         END DO      ! Next beam

C        *** Scale pressure and write to disk ***

         IF ( RUNTYP(1:1) .EQ. 'C' .OR.
     &        RUNTYP(1:1) .EQ. 'I' .OR.
     &        RUNTYP(1:1) .EQ. 'S' ) THEN
            CALL SCALEP( DALPHA, CV, R, U, NRD, NR,
     &                   RUNTYP, TOPOPT, FREQ )
            IREC  = 6 + NRD * ( IS - 1 )
            CALL WRTFLD( U, NRD, NR, IREC )
         ELSE IF ( RUNTYP(1:1) .EQ. 'A' ) THEN
C           --- Write out arrival data
            CALL WRTARR( NARR, AARR, DELARR, R, NRD, NR )
         END IF

      END DO    ! Next source depth

C     *** Display run time ***

      CALL RDTIME( TEND )
      WRITE( PRTFIL, 9500 ) TEND
 9500 FORMAT( /, ' CPU Time = ', G15.3 )

      STOP
      END
C**********************************************************************C
      FUNCTION PICKEPS( BEAMTYP, OMEGA, C, CZ,
     &                  ALPHA, DALPHA, RLOOP, EPMULT )

C     Picks the optimum value for epsilon

      INTEGER PRTFIL
      COMPLEX CI
      PARAMETER ( PRTFIL = 6, CI = ( 0.0, 1.0 ) )

      LOGICAL   INIFLG
      COMPLEX   PICKEPS, EPSOPT
      CHARACTER BEAMTYP*2, TAG*40

      SAVE INIFLG
      DATA INIFLG /.TRUE./

      IF ( BEAMTYP(1:1) .EQ. 'F' ) THEN
         TAG    = 'Space filling beams'
         HWIDTH = 2.0 / ( ( OMEGA / C ) * DALPHA )
         EPSOPT = CI * 0.5 * OMEGA * HWIDTH ** 2

      ELSE IF ( BEAMTYP(1:1) .EQ. 'M' ) THEN
         TAG    = 'Minimum width beams'
         HWIDTH = SQRT( 2.0 * C * 1000.0 * RLOOP / OMEGA )
         EPSOPT = CI * 0.5 * OMEGA * HWIDTH ** 2

      ELSE IF ( BEAMTYP(1:1) .EQ. 'W' ) THEN
         TAG    = 'WKB beams'
         IF ( CZ .EQ. 0.0 ) THEN
            EPSOPT = 1.0E10
         ELSE
            EPSOPT = ( -SIN( ALPHA ) / COS( ALPHA ** 2 ) ) * C * C / CZ
         ENDIF

      ELSE IF ( BEAMTYP(1:1) .EQ. 'C' ) THEN
         TAG    = 'Cerveny style beam'
      END IF

      PICKEPS = EPMULT * EPSOPT

C     *** On first call write info to prt file ***

      IF ( INIFLG ) THEN
         WRITE( PRTFIL, * ) TAG
         WRITE( PRTFIL, * ) 'HWIDTH = ', HWIDTH
         WRITE( PRTFIL, * ) 'EPSOPT = ', EPSOPT
         WRITE( PRTFIL, * ) 'EPMULT = ', EPMULT
         INIFLG = .FALSE.
      END IF

      RETURN
      END
C**********************************************************************C
      SUBROUTINE TRACE( DELTAS, RS, ZS, ALPHA,
     &                  BEAMTYP, ZBOX, RBOX, BOTOPT )

C     Traces the beam corresponding to a particular take-off angle

      INCLUDE 'COMBM.FOR'

      CHARACTER BOTOPT*2, BEAMTYP*2, BC*1

C     *** Initial conditions ***

      CALL SSP( RS, ZS, C, CR, CZ, CRR, CRZ, CZZ, TOPOPT, 'TAB' )

      CV(    1 ) = C

      RV(    1 ) = RS
      ZV(    1 ) = ZS
      RHOV(  1 ) = COS( ALPHA ) / C
      ZETAV( 1 ) = SIN( ALPHA ) / C

      P1V(   1 ) = 1.0
      Q1V(   1 ) = 0.0
      P2V(   1 ) = 0.0
      Q2V(   1 ) = 1.0

      TAUV(  1 ) = 0.0
      RFA(   1 ) = 1.0

C     *** Trace the beam ***
C     (note that REFLECT alters the step index I)

      ISEG = 1
      I = 0
      DO ISTEP = 1, MXN - 1
         I = I + 1

         CALL STEP(
     &      RV(   I   ),  ZV( I   ), RHOV( I   ), ZETAV( I   ),
     &      P1V(  I   ), Q1V( I   ),  P2V( I   ),   Q2V( I   ),
     &      TAUV( I   ), RFA( I   ),   CV( I   ),
     &      RV(   I+1 ),  ZV( I+1 ), RHOV( I+1 ), ZETAV( I+1 ),
     &      P1V(  I+1 ), Q1V( I+1 ),  P2V( I+1 ),   Q2V( I+1 ),
     &      TAUV( I+1 ), RFA( I+1 ),
     &      DELTAS, DEPTHT, DEPTHB, TOPOPT,
     &      RNG( ISEG ), DEP( ISEG ), TBR( ISEG ), TBZ( ISEG ) )

C        *** Reflections? ***

         D1 = RV( I + 1 ) - RNG( ISEG )      ! tangent pointing between bottom and ray term.
         D2 = ZV( I + 1 ) - DEP( ISEG )

         IF (      ZV( I + 1 ) .LE. 0.0    ) THEN   ! Surface reflection
            BC(1:1) = TOPOPT(2:2)
            CALL REFLECT( DEPTHT, I, BEAMTYP, BC,    CPT, RHOT, 'TOP',
     &      RNG( ISEG ), DEP( ISEG), TBR( ISEG ), TBZ( ISEG ) )
         ELSE IF ( TBR( ISEG ) * D2 - TBZ( ISEG ) * D1 .GE. 0.0 ) THEN   ! Bottom reflection
            CALL REFLECT( DEPTHB, I, BEAMTYP, BOTOPT, CPB, RHOB, 'BOT',
     &      RNG( ISEG ), DEP( ISEG), TBR( ISEG ), TBZ( ISEG ) )
         ENDIF

C        *** New bathymetry segment? ***

         IF ( RV( I + 1 ) .GT. RNG( ISEG + 1 ) .AND.
     &        ISEG .LT. NBTYPTS - 1 ) ISEG = ISEG + 1

C        *** Has the ray left the box, loss its energy,
C            or exceeded storage limit? ***
         IF ( RV( I + 1 ) .GT. RBOX .OR.
     &        ZV( I + 1 ) .GT. ZBOX .OR.
     &        RV( I + 1 ) .LT. 0.0  .OR.
     &        ABS( RFA( I + 1 ) ) .LT. 0.001  ) THEN
            NSTEPS = I + 1
            RETURN
         ELSE IF ( I .GE. MXN - 3 ) THEN
            CALL ERROUT( PRTFIL, 'W', 'TRACE',
     &             'Insufficient storage for ray trajectory' )
            NSTEPS = I
            RETURN
         END IF

      END DO   ! Next step

      RETURN
      END

C**********************************************************************C

      SUBROUTINE STEP(
     &   R0,   Z0, RHO0, ZETA0,
     &   P10, Q10,  P20,   Q20, TAU0, RFA0, C0,
     &   R2,   Z2, RHO2, ZETA2,
     &   P12, Q12,  P22,   Q22, TAU2, RFA2,
     &   DELTAS, DEPTHT, DEPTHB, TOPOPT,
     &   RNG1, DEP1, TBR, TBZ )

C     Does a single step along the ray

      REAL*8 R0, Z0, RHO0, ZETA0, P10, Q10, P20, Q20, TAU0,
     &       R1, Z1, RHO1, ZETA1, P11, Q11, P21, Q21,
     &       R2, Z2, RHO2, ZETA2, P12, Q12, P22, Q22, TAU2,
     &       RNG1, DEP1, TBR, TBZ
      COMPLEX*16 RFA0, RFA2
      CHARACTER TOPOPT*4

C     *** Phase 1 of modified polygon method (an Euler step) ***

      CALL SSP( SNGL( R0 ), SNGL( Z0 ), C0, CR0, CZ0, CRR0, CRZ0, CZZ0,
     &          TOPOPT, 'TAB' )

      CSQ0 = C0 * C0
      CNN0 = CSQ0 * ( CRR0 * ZETA0**2 -
     &          2.0 * CRZ0 * RHO0 * ZETA0 + CZZ0 * RHO0**2 )

      H = DELTAS

C     --- If crossing surface or bottom, reduce step to hit interface

      Z = Z0 + H * C0 * ZETA0
      R = R0 + H * C0 * RHO0

      IF ( Z .LT. DEPTHT ) H = 1.001 * ( DEPTHT - Z0 ) / ( ZETA0 * C0 )

      D1 = R - RNG1      ! tangent pointing between bottom and ray term.
      D2 = Z - DEP1
      IF ( TBR * D2 - TBZ * D1 .GT. 0.00001 ) THEN
         A1 = -C0 * RHO0   ! ray tangent
         A2 = -C0 * ZETA0
         E1 = R0 - RNG1    ! tangent between bottom and ray org.
         E2 = Z0 - DEP1
         H = 1.001 * ( E1 * TBZ - E2 * TBR ) / ( A1 * TBZ - A2 * TBR )
      ENDIF

      HH = 0.5 * H

      R1    = R0    + HH * C0 *  RHO0
      Z1    = Z0    + HH * C0 * ZETA0
      RHO1  = RHO0  - HH * CR0 / CSQ0
      ZETA1 = ZETA0 - HH * CZ0 / CSQ0
      P11   = P10   - HH * CNN0 / CSQ0 * Q10
      Q11   = Q10   + HH * C0          * P10
      P21   = P20   - HH * CNN0 / CSQ0 * Q20
      Q21   = Q20   + HH * C0          * P20

C     *** Phase 2 of modified polygon method ***

      CALL SSP( SNGL( R1 ), SNGL( Z1 ), C1, CR1, CZ1, CRR1, CRZ1, CZZ1,
     &          TOPOPT, 'TAB' )

      CSQ1 = C1 * C1
      CNN1 = CSQ1 * ( CRR1 * ZETA1**2 -
     &               2.0 * CRZ1 * RHO1 * ZETA1 + CZZ1 * RHO1**2 )

C     --- If crossing surface or bottom, reduce step to hit interface

      Z = Z0 + H * C1 * ZETA1
      R = R0 + H * C1 * RHO1

      IF ( Z .LE. DEPTHT )  H = 1.001 * ( DEPTHT - Z0 ) / ( ZETA1 * C1 )

      D1 = R - RNG1
      D2 = Z - DEP1
      IF ( TBR * D2 - TBZ * D1 .GT. 0.00001 ) THEN
         A1 = -C1 * RHO0
         A2 = -C1 * ZETA0
         E1 = R0 - RNG1
         E2 = Z0 - DEP1
         H = 1.001 * ( E1 * TBZ - E2 * TBR ) / ( A1 * TBZ - A2 * TBR )
      ENDIF

      R2    = R0    + H * C1 *  RHO1
      Z2    = Z0    + H * C1 * ZETA1
      RHO2  = RHO0  - H * CR1 / CSQ1
      ZETA2 = ZETA0 - H * CZ1 / CSQ1
      P12   = P10   - H * CNN1 / CSQ1 * Q11
      Q12   = Q10   + H * C1          * P11
      P22   = P20   - H * CNN1 / CSQ1 * Q21
      Q22   = Q20   + H * C1          * P21
      TAU2  = TAU0  + H / C1
      RFA2  = RFA0

      RETURN
      END
C**********************************************************************C
      SUBROUTINE REFLECT( ZREFL, I, BEAMTYP, BC, CHS, RHOHS, BOTTOP,
     &   RNG1, DEP1, TBR1, TBZ1 )

      INCLUDE 'COMBM.FOR'
      REAL*8 RHOHS, THEINT, RINT, PHIINT, RNG1, DEP1, TBR1, TBZ1
      COMPLEX GAMMA1, GAMMA2, GAMMA1SQ, GAMMA2SQ, GK, REFL
      COMPLEX*16 CHS
      CHARACTER BEAMTYP*2, BC*1, BOTTOP*3

      I = I + 1

C     *** here's the geometric part, changing the ray direction ***

      RV(    I+1 ) = RV( I )
      ZV(    I+1 ) = ZV( I )

      IF ( BOTTOP .EQ. 'BOT' ) THEN
C        --- component of ray tangent normal to bottom.
C            normal vector is ( -TBZ, TBR ) in ( r, z) coords.
         BETA = -RHOV( I ) * TBZ1 + ZETAV( I ) * TBR1 
         RHOV(  I+1 ) =  RHOV( I ) + 2.0 * BETA * TBZ1
         ZETAV( I+1 ) = ZETAV( I ) - 2.0 * BETA * TBR1
      ELSE
         RHOV(  I+1 ) =   RHOV( I )
         ZETAV( I+1 ) = -ZETAV( I )
      ENDIF

      CALL SSP( SNGL( RV( I + 1 ) ), SNGL( ZV( I + 1 ) ), C, CR, CZ,
     &                  CRR, CRZ, CZZ, TOPOPT, 'TAB' )

      CV( I ) = C

C     --- Following needs updating for c(r,z) problem

      RN = -2.0 * RHOV( I + 1 )**2 / ZETAV( I + 1 ) * CZ / C

      IF ( BEAMTYP(2:2) .EQ. 'D' ) RN = 2.0 * RN
      IF ( BEAMTYP(2:2) .EQ. 'Z' ) RN = 0.0

      P1V( I + 1 ) = P1V( I ) + Q1V( I ) * RN
      Q1V( I + 1 ) = Q1V( I )
      P2V( I + 1 ) = P2V( I ) + Q2V( I ) * RN
      Q2V( I + 1 ) = Q2V( I )

C     --- account for phase change

      IF ( BC .EQ. 'R' ) THEN                  ! rigid
         TAUV( I + 1 ) = TAUV( I )
         RFA(  I + 1 ) = RFA(  I )
      ELSE  IF ( BC .EQ. 'V' ) THEN            ! vacuum
	 TAUV( I + 1 ) =  TAUV( I )
         RFA(  I + 1 ) = -RFA(  I )
      ELSE  IF ( BC .EQ. 'F' ) THEN            ! file
         RADDEG = 180.0 / PI
 	 THEINT = RADDEG * ATAN2( ZETAV( I ), RHOV( I ) )   ! angle of incidence
C        --- Evaluate R(THEINT)
         IF ( BOTTOP(1:3) .EQ. 'TOP' ) THEN
            CALL REFCO( THEINT, RINT, PHIINT,
     &                  THETATOP, RTOP, PHITOP, NTOPPTS )
         ELSE
            CALL REFCO( THEINT, RINT, PHIINT,
     &                  THETABOT, RBOT, PHIBOT, NBOTPTS )
         ENDIF

         TAUV( I + 1 ) = TAUV( I )
         RFA(  I + 1 ) = RFA(  I ) * RINT * EXP( CI * PHIINT )

      ELSE  IF ( BC .EQ. 'A' ) THEN            ! half-space
	 GK     = OMEGA * RHOV( I )   ! hor. wavenumber
	 GAMMA1SQ = ( OMEGA / C   ) ** 2 - GK ** 2
	 GAMMA2SQ = ( OMEGA / CHS ) ** 2 - GK ** 2
	 GAMMA1 = SQRT( GAMMA1SQ )
	 GAMMA2 = SQRT( GAMMA2SQ )

         REFL = ( RHOHS * GAMMA1 - GAMMA2 ) /
     &	        ( RHOHS * GAMMA1 + GAMMA2 )

         IF ( ABS( REFL ) .LE. 1.0E-5 ) THEN
C           --- kill a ray that has lost its energy in reflection
            TAUV( I + 1 ) = TAUV( I )
            RFA(  I + 1 ) = 0.0
         ELSE
            TAUV( I + 1 ) = TAUV( I )
            RFA(  I + 1 ) = REFL * RFA(  I )
C           --- compute beam-displacement Tindle, Eq. (14)
C           --- needs a correction to beam-width as well ...
C            IF ( REAL( GAMMA2SQ ) .LT. 0.0 ) THEN
C               RHOW   = 1.0   ! density of water
C               RHOWSQ  = RHOW  * RHOW
C               RHOHSSQ = RHOHS * RHOHS
C               DELTA = 2 * GK * RHOW * RHOHS * ( GAMMA1SQ - GAMMA2SQ ) /
C     &            ( GAMMA1 * CI * GAMMA2 *
C     &            ( -RHOWSQ * GAMMA2SQ + RHOHSSQ * GAMMA1SQ ) )
C               RV( I + 1 ) = RV( I + 1 ) + DELTA
C            END IF
	 ENDIF
      ENDIF

      RETURN
      END
C**********************************************************************C
      SUBROUTINE INFLUR( U, DELTAR, EPS, ALPHA, NIMAGE,
     &   IBWIN2, RUNTYP, RADMAX, BEAMTYP )

C     Computes the beam influence, i.e. 
C     the contribution of a single beam to the complex pressure

      INCLUDE 'COMBM.FOR'

      LOGICAL IOK
      INTEGER KMAHV( MXN )
      REAL    NA, NB, NSQ
      COMPLEX PV( MXN ), QV( MXN ), P, Q, EPS, CONTRI, U( NRD, NR ),
     &        GAMMAV( MXN ), GAMMA
      CHARACTER RUNTYP*2, BEAMTYP*2

C     Note that during reflection imag(q) is constant and
C     adjacent normals cannot bracket a segment of the TL
C     line, so no special treatment is necessary

C     ------ Lloyd mirror pattern:

      DS = 2.0*( SIN( OMEGA * ZV( 1 ) * ZETAV( 1 ) ) )**2

C     *** BEGIN BY FORMING (P, Q) AND KMAH INDEX ***

C     Note treatment of KMAH index is incorrect for 'Cerveny'
C     style beam width RUNTYP

      IF ( BEAMTYP(1:1) .EQ. 'C' ) EPS = CI * ABS( Q1V( 1 ) / Q2V( 1 ) )

      PV( 1 ) = P1V( 1 ) + EPS * P2V( 1 )
      QV( 1 ) = Q1V( 1 ) + EPS * Q2V( 1 )

C     ------ Following is sloppy

      GAMMAV( 1 ) = 0.0
      IF ( QV( 1 ) .NE. 0 ) GAMMAV( 1 ) = PV( 1 ) / QV( 1 )
      KMAHV( 1 ) = 1

      DO 1000 I = 2, NSTEPS
         IF ( BEAMTYP(1:1) .EQ. 'C' ) EPS = CI * ABS( Q1V(I) / Q2V(I) )

         PV( I ) = P1V( I ) + EPS * P2V( I )
         QV( I ) = Q1V( I ) + EPS * Q2V( I )

         GAMMAV( I ) = PV( I ) / QV( I )
         KMAHV(  I ) = KMAHV( I - 1 )

         CALL BRCUT( QV( I - 1 ), QV( I ), BEAMTYP, KMAHV( I ) )

 1000 CONTINUE

C     Loop over receivers

      DO 9000 IRD = 1, NRD
      ZR = RD( IRD )

C     *** Loop over images ***

      DO 8000 IMAGE = 1, NIMAGE
         IR1 = 9999

C        Loop over steps

         DO 5000 I = 1, NSTEPS

C           *** Compute ray-centered coordinates, (ZNV, RNV) ***

            ZNV = -RHOV( I ) * CV( I )
C           ------ Check for normal parallel to TL-line
            IF ( ABS( ZNV ) .LT. 1.0E-5 ) THEN
               IOK = .FALSE.
               GOTO 5000
            ENDIF

C           Images of beams

            IF ( IMAGE .EQ. 1 ) THEN
C              ------ True beam
               RNV = ZETAV( I ) * CV( I )
               NB  = ( ZR - ZV( I ) ) / ZNV

            ELSE IF ( IMAGE .EQ. 2 ) THEN
C              ------ Surface reflected beam
               RNV = -ZETAV( I ) * CV( I )
               NB  = ( ZR - ( 2.0 * DEPTHT - ZV( I ) ) ) / ZNV

            ELSE IF ( IMAGE .EQ. 3 ) THEN
C              ------ Bottom reflected beam
               RNV = -ZETAV( I ) * CV( I )
               NB  = ( ZR - ( 2.0 * DEPTHB - ZV( I ) ) ) / ZNV
            ENDIF

            RB = RV( I ) + NB * RNV

C           *** Compute influence for each rcvr ***

            IF ( RB .GT. R( 1 ) ) THEN
               IR2 = MIN( INT( ( RB - R( 1 ) ) / DELTAR ) + 1, NR )
            ELSE
               IR2 = MIN( INT( ( RB - R( 1 ) ) / DELTAR ), NR )
            ENDIF

            IF ( ( IR2 .GT. IR1 ) .AND. ( RV( I ) .NE. RV( I - 1 ) )
     &         .AND. IOK ) THEN
               DO 2000 IR = IR1 + 1, IR2
                  W = ( R( IR ) - RA ) / ( RB - RA )

                  P     =     PV(I-1) + W * (     PV(I) -     PV(I-1) )
                  Q     =     QV(I-1) + W * (     QV(I) -     QV(I-1) )
                  GAMMA = GAMMAV(I-1) + W * ( GAMMAV(I) - GAMMAV(I-1) )

                  NSQ = ( NA + W * ( NB - NA ) )**2

                  IF ( AIMAG( GAMMA ) .GT. 0 ) THEN
                     WRITE( PRTFIL, * ) 'Unbounded beam'
                     WRITE( PRTFIL, * ) PV(I-1),QV(I-1),PV(I-1)/QV(I-1)
                     WRITE( PRTFIL, * ) PV(I), QV(I), PV(I) / QV(I)
                     WRITE( PRTFIL, * ) GAMMA, W
                  ENDIF

C                 ------ Within beam window?
                  IF ( -0.5*OMEGA*AIMAG(GAMMA)*NSQ .LT. IBWIN2 ) THEN

                     C   =   CV( I-1 ) + W * (   CV( I ) -   CV( I-1 ) )
                     TAU = TAUV( I-1 ) + W * ( TAUV( I ) - TAUV( I-1 ) )

                     IF ( RUNTYP(2:2) .EQ. 'C' )
     &                  EPS = CI * ABS( Q1V( I ) / Q2V( I ) )

                     CONTRI = SQRT( C * ABS( EPS ) / Q * COS( ALPHA ) )*
     &                  EXP( -CI * OMEGA * ( TAU + 0.5 * GAMMA * NSQ ) )

C                    ------ Get correct branch of SQRT
                     KMAH = KMAHV( I - 1 )
                     CALL BRCUT( QV( I - 1 ), Q, BEAMTYP, KMAH )

                     IF ( KMAH  .LT. 0 ) CONTRI = -CONTRI
                     IF ( IMAGE .EQ. 2 ) CONTRI = -CONTRI

C                    ------ Coherent, semi-coherent, or incoherent TL:

                     IF ( RUNTYP(1:1) .EQ. 'I' ) CONTRI =    ABS(CONTRI)
                     IF ( RUNTYP(1:1) .EQ. 'S' ) CONTRI = DS*ABS(CONTRI)

                     U( IRD, IR ) = U( IRD, IR ) +
     &                   HERMIT( NA, RADMAX, 2.0 * RADMAX ) * CONTRI

                  ENDIF
 2000          CONTINUE
            ELSE
C              ------ Next step maybe ok
               IOK = .TRUE.
            ENDIF
            RA  = RB
            NA  = NB
            IR1 = MAX( IR2, 0 )
 5000    CONTINUE       ! Next step

 8000 CONTINUE      ! Next image

 9000 CONTINUE      ! Next receiver depth

      RETURN
      END
C**********************************************************************C
      SUBROUTINE INFLUC( U, DELTAR, EPS, ALPHA, NIMAGE,
     &   IBWIN2, RUNTYP, RADMAX, BEAMTYP )

C     Computes the beam influence, i.e. 
C     the contribution of a single beam to the complex pressure
C     This version uses a beam representation in Cartesian coordinates

      INCLUDE 'COMBM.FOR'

      INTEGER   KMAHV( MXN )
      COMPLEX   PV( MXN ), QV( MXN ), Q, EPS, CONTRI, U( NRD, NR ),
     &          GAMMAV( MXN ), GAMMA, CONST
      CHARACTER RUNTYP*2, BEAMTYP*2

C     Note that during reflection imag(q) is constant and
C     adjacent normals cannot bracket a segment of the TL
C     line, so no special treatment is necessary

C     ------ Lloyd mirror pattern:
      DS = 2.0 * ( SIN( OMEGA * ZV(1) * ZETAV( 1 ) ) )**2

C     *** Begin by forming (p, q) and KMAH index ***

C     Note treatment of KMAH index is incorrect for 'Cerveny'
C     style beam width BEAMTYP

      DO 1000 I = 1, NSTEPS
         IF ( BEAMTYP(1:1) .EQ. 'C' ) EPS = CI * ABS( Q1V(I) / Q2V(I) )

         PV( I ) = P1V( I ) + EPS * P2V( I )
         QV( I ) = Q1V( I ) + EPS * Q2V( I )

C         RLTEMP = SQRT( -2.0 / ( OMEGA * AIMAG( PV( I ) / QV( I ) ) ) )
C         RKTEMP = -CV( I ) * REAL( PV( I ) / QV( I ) )

C         WRITE( PRTFIL, * ) RLTEMP, RKTEMP

C         IF ( BEAMTYP(2:2) .EQ. 'D' ) THEN
C            PV( I ) = REAL( PV( I ) )
C            P = PV( I )
C            L = 30.0 * 1500.0 / OMEGA
C            QV( I ) = ( CI * OMEGA * L**2 / 2.0) * P
C         ENDIF

         TR =  RHOV( I ) * CV(I)
         TZ = ZETAV( I ) * CV(I)

         CALL SSP( SNGL( RV( I ) ), SNGL( ZV( I ) ), C, CR, CZ,
     &             CRR, CRZ, CZZ, TOPOPT, 'TAB')

         CSQ = C * C
         CS = CR * TR + CZ * TZ
         CN = CR * TZ - CZ * TR 

         GAMMAV( I ) = 0.0
         IF ( QV( I ) .NE. 0.0 )
     &      GAMMAV( I ) = 0.5 * ( PV(I) / QV(I) * TR**2 +
     &                    2.0 * CN / CSQ * TZ * TR - CS / CSQ * TZ**2 )

         IF ( I .EQ. 1 ) THEN
            KMAHV( 1 ) = 1
         ELSE
            KMAHV( I ) = KMAHV( I - 1 )
            CALL BRCUT( QV( I - 1 ), QV( I ), BEAMTYP, KMAHV( I ) )
         ENDIF

 1000 CONTINUE

C     *** Loop over steps ***

      DO 5000 I = 3, NSTEPS
         IF ( RV( I ) .GT. R( NR ) ) RETURN
         RA = RV( I-1 )
         RB = RV( I )
         IF ( RA .EQ. RB ) GOTO 5000

C        ------ Compute upper index on rcvr line
         IR1 = MAX( MIN( INT( ( RA - R(1) ) / DELTAR ) + 1, NR ), 1 )
         IR2 = MAX( MIN( INT( ( RB - R(1) ) / DELTAR ) + 1, NR ), 1 )

C        ------ Following for projecting back reflected beams
C         IF ( RV( I-2 ) .EQ. RV( I-1 ) ) IR1 = 1
C         IF ( RV( I )   .EQ. RV( I+1 ) ) IR2 = NR

         IF ( IR2 .GT. IR1 ) THEN
            DO 2000 IR = IR1 + 1, IR2
               W = ( R( IR ) - RA ) / ( RB - RA )

               Z     =     ZV(I-1) + W * (     ZV(I) -     ZV(I-1) )
               C     =     CV(I-1) + W * (     CV(I) -     CV(I-1) )
               Q     =     QV(I-1) + W * (     QV(I) -     QV(I-1) )
               TAU   =   TAUV(I-1) + W * (   TAUV(I) -   TAUV(I-1) )
               ZETA  =  ZETAV(I-1) + W * (  ZETAV(I) -  ZETAV(I-1) )
               GAMMA = GAMMAV(I-1) + W * ( GAMMAV(I) - GAMMAV(I-1) )

               IF ( AIMAG( GAMMA ) .GT. 0 ) THEN
                  WRITE( PRTFIL, * ) 'Unbounded beam'
                  WRITE( PRTFIL, * ) GAMMAV(I-1), GAMMAV(I), GAMMA
               ENDIF

               IF ( BEAMTYP(1:1) .EQ. 'C' )
     &             EPS = CI * ABS( Q1V( I ) / Q2V( I ) )

               CONST = SQRT( C * ABS(   EPS   ) / Q * COS( ALPHA ) )

C              ------ Get correct branch of SQRT
               KMAH = KMAHV( I - 1 )

               CALL BRCUT( QV( I - 1 ), Q, BEAMTYP, KMAH )

               IF ( KMAH .LT. 0 ) CONST = -CONST

C              *** Loop over receivers ***

               DO 9000 IRD = 1, NRD
                  ZR = RD( IRD )

C                 ------ True beam
                  DELTAZ = ZR - Z
C                  IF ( OMEGA * AIMAG( GAMMA )*DELTAZ**2 .LT. IBWIN2 )
                   CONTRI = HERMIT( DELTAZ, 0.0, RADMAX ) *
     &                EXP( -CI * OMEGA *
     &                ( TAU + ZETA * DELTAZ + GAMMA * DELTAZ**2) )

C                 ------ Surface reflected beam
                  IF ( NIMAGE .GE. 2 ) THEN
                     DELTAZ = -ZR + 2.0 * DEPTHT - Z
                     IF ( OMEGA * AIMAG( GAMMA )*DELTAZ**2 .LT. IBWIN2 )
     &                  CONTRI =  CONTRI -
     &                     HERMIT( DELTAZ, 0.0, RADMAX ) *
     &                     EXP( -CI * OMEGA *
     &                     ( TAU + ZETA * DELTAZ + GAMMA * DELTAZ**2) )
                  ENDIF

C                 ------ Bottom reflected beam
                  IF ( NIMAGE .GE. 3 ) THEN
                     DELTAZ = -ZR + 2.0 * DEPTHB - Z
                     IF ( OMEGA * AIMAG( GAMMA )*DELTAZ**2 .LT. IBWIN2 )
     &                  CONTRI =  CONTRI + 
     &                     HERMIT( DELTAZ, 0.0, RADMAX ) *
     &                     EXP( -CI * OMEGA *
     &                     ( TAU + ZETA * DELTAZ + GAMMA * DELTAZ**2) )
                  ENDIF

C                 ------ Coherent, semi-coherent, or incoherent TL:

                  CONTRI = CONST * CONTRI

                  IF ( RUNTYP(1:1) .EQ. 'I' ) THEN
                     CONTRI =      ABS( CONTRI )
                  ELSE IF ( RUNTYP(1:1) .EQ. 'S' ) THEN
                     CONTRI = DS * ABS( CONTRI )
                  ENDIF

                  U( IRD, IR ) = U( IRD, IR ) + CONTRI

 9000          CONTINUE      ! Next receiver depth
 2000       CONTINUE      ! Next receiver range
         ENDIF
 5000 CONTINUE      ! Next step along the ray

      RETURN
      END
C**********************************************************************C
      SUBROUTINE INFLUG( U, ZS, ALPHA, RUNTYP, DALPHA, NSD,
     &   NARR, DELARR, AARR, MXNARR )

C     Computes the beam influence, i.e. 
C     the contribution of a single beam to the complex pressure
C     This version uses a beam representation in Cartesian coordinates

      INCLUDE 'COMBM.FOR'

      LOGICAL   NEWRAY
      INTEGER     NARR( NRD, NR )
      REAL      DELARR( NRD, NR, * )
      COMPLEX     AARR( NRD, NR, * )

      REAL*8    QV( MXN )
      COMPLEX   U( NRD, NR ), CONST, CONTRI, CAUSTIC
      CHARACTER RUNTYP*2

      EQUIVALENCE ( QV, Q1V )

C     --- Lloyd mirror pattern (pressure form):
      DS = SQRT( 2.0 ) * SIN( OMEGA * ZS * ZETAV( 1 ) )
      Q0 = CV( 1 ) / DALPHA   ! Reference for J = Q0 / Q

      CAUSTIC = 1   ! CAUSTIC is multiplied by i after each caustic
      QOLD    = 1.0
      Q1METER = CV( 1 )
      RAT1    = SQRT( Q1METER / CV ( 1 ) * COS( ALPHA ) )

C     *** Loop over steps ***

      RA = RV( 1 )
      IR = 1

      DO I = 2, NSTEPS

         RB = RV( I )

C        ------ Loop over bracketted receiver ranges

1000     CONTINUE

         IF ( RB .NE. RA .AND. RB .GT. R( IR ) ) THEN

            W = ( R( IR ) - RA ) / ( RB - RA )

            Z     =     ZV(I-1) + W * (     ZV(I) -     ZV(I-1) )
            C     =     CV(I-1) + W * (     CV(I) -     CV(I-1) )
            Q     =     QV(I-1) + W * (     QV(I) -     QV(I-1) )
            RHO   =   RHOV(I-1) + W * (   RHOV(I) -   RHOV(I-1) )
            ZETA  =  ZETAV(I-1) + W * (  ZETAV(I) -  ZETAV(I-1) )
            TAU   =   TAUV(I-1) + W * (   TAUV(I) -   TAUV(I-1) )

            TR = RHO * C
            A  = ABS( TR * Q0 / Q )
            RADMAX = 1.0 / A

            IF ( Q .LT. 0.0 .AND. QOLD .GE. 0.0 .OR.
     &           Q .GT. 0.0 .AND. QOLD .LE. 0.0 )
     &                                   CAUSTIC = CI * CAUSTIC

C           --- Coherent, semi-coherent, or incoherent TL:
            CONST = RAT1 * SQRT( C / ABS( Q ) )
            IF ( RUNTYP(1:1) .EQ. 'S' ) CONST = DS * CONST

C           *** Loop over receiver depths ***

            DO IRD = 1, NRD
               DELTAZ =  RD( IRD ) - Z   ! ray to rcvr distance 
               ADZ    = ABS( DELTAZ )

               IF ( ADZ .LT. RADMAX ) THEN
                 IF ( RUNTYP(1:1) .EQ. 'E' ) THEN
C                  *** Write the eigenray trajectory to RAYFIL ***
                   CALL WRTRAY( RV, ZV, NSTEPS, DEPTHT, DEPTHB )

                 ELSE IF ( RUNTYP(1:1) .EQ. 'A' ) THEN
C                  *** arrivals ***
                   CONTRI = CONST * A * ( RADMAX - ADZ )
                   DELAY  = TAU + ZETA * DELTAZ

C                  --- Is this the second bracketting ray of a pair?
                   NT =  NARR( IRD, IR )    ! # of arrivals
                   NEWRAY = .TRUE.
                   IF ( NT .GE. 1 ) THEN
                     IF( OMEGA * ABS( DELAY - DELARR( IRD, IR, NT ) )
     &                  .LT. 1.0 ) NEWRAY = .FALSE.
                   ENDIF

                   IF ( NEWRAY ) THEN
                     NT =  NT + 1   ! # of arrivals
                     IF ( NT .LE. MXNARR ) THEN
                       NARR( IRD, IR     ) = NT
                       AARR( IRD, IR, NT ) = CONTRI * CAUSTIC * RFA( I )
                       DELARR( IRD, IR, NT ) = DELAY
                     ENDIF
                   ELSE
                     AARR(   IRD, IR, NT ) = AARR( IRD, IR, NT ) +
     &                  CONTRI * CAUSTIC * RFA( I )
                     DELARR( IRD, IR, NT ) =
     &                   0.5 * ( DELARR( IRD, IR, NT ) + DELAY )
                   ENDIF
                 ELSE IF ( RUNTYP(1:1) .EQ. 'C'  ) THEN
                   CONTRI = CONST * A * ( RADMAX - ADZ )
                   DELAY  = TAU + ZETA * DELTAZ
                   CONTRI = CAUSTIC * CONTRI * RFA( I ) *
     &                  EXP( -CI * OMEGA * DELAY )
                    U( IRD, IR ) = U( IRD, IR ) + CONTRI
                 ELSE
C                  *** incoherent field calculation ***
                   CONTRI = CONST ** 2 * A * ( RADMAX - ADZ )
     &                        * ABS( RFA( I ) )
                   U( IRD, IR ) = U( IRD, IR ) + CONTRI
                 ENDIF
               ENDIF
            END DO      ! Next receiver depth

            QOLD = Q
            IR = IR + 1
            IF ( IR .GT. NR ) RETURN
            GOTO 1000     ! Next receiver range
         ENDIF

         RA = RB
      END DO      ! Next step along the ray

      RETURN
      END
C**********************************************************************C
      SUBROUTINE BRCUT( Q1C, Q2C, BEAMTYP, KMAH )

C     Checks for a branch cut crossing

      COMPLEX   Q1C, Q2C
      CHARACTER BEAMTYP*2

      IF ( BEAMTYP(1:1) .EQ. 'W' ) THEN
C        ------ WKBeams:
         Q1 = REAL( Q1C )
         Q2 = REAL( Q2C )
         IF ( ( Q1 .LT. 0.0 .AND. Q2 .GE. 0.0 ) .OR.
     &        ( Q1 .GT. 0.0 .AND. Q2 .LE. 0.0 ) ) KMAH = -KMAH
      ELSE
C        ------ All other beam types
         IF ( REAL( Q2C ) .LT. 0.0 ) THEN
            Q1 = AIMAG( Q1C )
            Q2 = AIMAG( Q2C )
            IF ( ( Q1 .LT. 0.0 .AND. Q2 .GE. 0.0 ) .OR.
     &           ( Q1 .GT. 0.0 .AND. Q2 .LE. 0.0 ) ) KMAH = -KMAH
         ENDIF
      ENDIF

      RETURN
      END
C**********************************************************************C
      FUNCTION HERMIT( X, X1, X2 )

C     Calculates a smoothing function based on the h0 hermite cubic

      AX  = ABS( X  )
      AX1 = ABS( X1 )
      AX2 = ABS( X2 )

      IF ( AX .LE. AX1 ) THEN
         HERMIT = 1.0
      ELSE IF ( AX .GE. AX2 ) THEN
         HERMIT = 0.0
      ELSE
         U = ( AX - AX1 ) / ( AX2 - AX1 )
         HERMIT = ( 1.0 + 2.0 * U ) * ( 1.0 - U )**2
      ENDIF

      hermit = 1.0
C      hermit = hermit / ( 0.5 * ( x1 + x2 ) )

      RETURN
      END
C**********************************************************************C
      SUBROUTINE SCALEP( DALPHA, C, R, U, NRD, NR,
     &                   RUNTYP, TOPOPT, FREQ )

C     Scale the pressure field

      REAL      R( NR )
      COMPLEX   U( NRD, NR )
      CHARACTER TOPOPT*4, RUNTYP*2

C     --- Compute scale factor for field
      IF ( RUNTYP(2:2) .NE. 'G' ) THEN
         CONST = -DALPHA * SQRT( FREQ ) / C
      ELSE
         CONST = -1.0
      ENDIF

C     *** Thorpe attenuation? ***

      IF ( TOPOPT(3:3) .EQ. 'T' ) THEN
         F2 = ( FREQ / 1000.0 ) ** 2
         ALPHA = 40.0 * F2 / ( 4100.0 + F2 ) + 0.1 * F2 / ( 1.0 + F2 )
         ALPHA = ALPHA / 914.4                 ! dB / m
         ALPHA = ALPHA / 8.6858896             ! Nepers / m
      ELSE
         ALPHA = 0.0
      ENDIF

C     *** Loop over range ***

      DO IR = 1, NR

         FACTOR = CONST * EXP( -ALPHA * R( IR ) ) / SQRT( R( IR ) )

C        *** Loop over depth ***

         DO IRD = 1, NRD
C           --- For incoherent RUNTYP, convert intensity to pressure
            IF ( RUNTYP(1:1) .NE. 'C' )
     &         U( IRD, IR ) = SQRT( REAL( U( IRD, IR ) ) )
            
            U( IRD, IR ) = FACTOR * U( IRD, IR )
         END DO

      END DO

      RETURN
      END

C**********************************************************************C

      SUBROUTINE WRTRAY( RV, ZV, NSTEPS, DEPTHT, DEPTHB )

C     Compress the ray data keeping every ISKIP point,
C     points near surface or bottom, and last point.
C     Write to RAYFIL.

      INTEGER RAYFIL
      PARAMETER ( RAYFIL = 21 )

      REAL*8 RV( * ), ZV( * )

C     --- compression
      N2 = 1
      ISKIP = MAX( NSTEPS / 5000, 1 )   ! max #pts written is 500

      DO I = 2, NSTEPS
         IF ( MIN( DEPTHB - ZV( I ),  ZV( I ) - DEPTHT ) .LT. 0.2 .OR.
     &        MOD( I, ISKIP ) .EQ. 0 .OR. I .EQ. NSTEPS ) THEN
             N2 = N2 + 1
             RV( N2 ) = RV( I )
             ZV( N2 ) = ZV( I )
         END IF
      END DO

C     --- write to ray file
      WRITE( RAYFIL, * ) N2
      DO I = 1, N2
         WRITE( RAYFIL, * ) SNGL( RV( I ) ), SNGL( ZV( I ) )
      ENDDO

      RETURN
      END

C**********************************************************************C

      SUBROUTINE WRTARR( NARR, AARR, DELARR, R, NRD, NR )

C     Writes the arrival data (amplitude, delay for each eigenray)

      INTEGER ARRFIL
      PARAMETER ( ARRFIL = 36 )

      INTEGER     NARR( NRD, NR )
      REAL      DELARR( NRD, NR, * ), R( NR )
      COMPLEX     AARR( NRD, NR, * ), AARRT

      DO IRD = 1, NRD
        DO IR = 1, NR
          WRITE( ARRFIL, * ) NARR( IRD, IR )
          DO JJ = 1, NARR( IRD, IR )
            AARRT = AARR( IRD, IR, JJ ) / SQRT( R( IR ) )   ! cylindrical spreading
            WRITE( ARRFIL, * )
     &         REAL( AARRT ), AIMAG( AARRT ),
     &         DELARR( IRD, IR, JJ )
          END DO   ! next arrival
        END DO   ! next receiver depth
      END DO   ! next range

      RETURN
      END
