      PROGRAM BEAM3D

C     Gaussian beam tracing in three dimensions
C     Michael B. Porter

      IMPLICIT REAL ( A-H, O-Z )
      INTEGER BTPFIL, PRTFIL, RAYFIL
      COMPLEX CI
      PARAMETER ( BTPFIL = 5, PRTFIL = 6, RAYFIL = 21,
     &   PI = 3.14159265, CI = (0.0, 1.0), DEGRAD = PI / 180.0, 
     &   MAXNR = 1000, MAXNT = 50, MAXNRT = 50000 )

      REAL*4 R( MAXNR ), THETA( MAXNT )
      COMPLEX EPS1, EPS2
      COMPLEX*8 P( MAXNRT )
      CHARACTER TITLE*80, OPTION*8, OPT2*4

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

C      CALL SYS$GETTIME( TSTART, WTIME1 )
      READ( BTPFIL, * ) TITLE
      READ( BTPFIL, * ) OPTION
      READ( BTPFIL, * ) FREQ
      READ( BTPFIL, * ) ALPHA0, ALPHA1, NALPHA
      READ( BTPFIL, * ) BETA0, BETA1, NBETA
      READ( BTPFIL, * ) IBWIN
      READ( BTPFIL, * ) NSTEPS, DELTAS
      READ( BTPFIL, * ) EPMULT
      READ( BTPFIL, * ) ZS, ZR
      READ( BTPFIL, * ) THETA0, THETA1, NTHETA
      READ( BTPFIL, * ) R0, R1, NR
      READ( BTPFIL, * ) OPT2

C     ------ For Nx2D run beams must lie on rcvr radial
      IF ( OPT2 .EQ. 'NX2D' ) THEN
         BETA0 = THETA0
         BETA1 = THETA1
         NBETA = NTHETA
      ENDIF

      TITLE = 'BEAM'//OPT2//'- '//TITLE

C     *** Echo to print file ***

      WRITE( PRTFIL, 2000 ) TITLE, FREQ, ALPHA0, ALPHA1, NALPHA, 
     &   BETA0, BETA1, NBETA, IBWIN, 
     &   NSTEPS, DELTAS, EPMULT, ZS, ZR, 
     &   THETA0, THETA1, NTHETA, R0, R1, NR

 2000 FORMAT( ' ', A80, //, ' FREQ = ', G14.6, //, 
     &   ' BEAM FAN PARAMETERS:', /, 
     &   ' ALPHA0 = ', G14.6, ' ALPHA1 = ', G14.6, ' NALPHA = ', I5, /,
     & '  BETA0 = ', G14.6, '  BETA1 = ', G14.6, '  NBETA = ', I5, ///,
     &   ' IBWIN = ', I5, ///,
     &   ' NSTEPS = ', I5, ' DELTAS = ', G14.6, //,
     &  ' EPMULT = ', G14.6, //, ' ZS = ', G14.6, ' ZR = ', G14.6, ///,
     &   ' RECEIVER LINES:', /,
     &   ' THETA0 = ', G14.6, ' THETA1 = ', G14.6, ' NTHETA = ', I5, /
     &   '     R0 = ', G14.6, '     R1 = ', G14.6, ' NR = ', I5, /// )
      WRITE( PRTFIL, * ) 'SOLUTION TECHNIQUE:', OPT2

      IF ( OPTION( 1:1 ) .EQ. 'R' ) THEN
         WRITE( PRTFIL, * ) 'RAY FILE OPTION INVOKED'
         WRITE( RAYFIL ) TITLE, FREQ, NALPHA, NSTEPS
      ENDIF

      IF ( OPTION( 3:3 ) .EQ. 'D' ) THEN
         WRITE( PRTFIL, * ) '****** CURVATURE DOUBLING INVOKED'
      ENDIF
      WRITE( PRTFIL, * )

      R0 = 1000.0 * R0
      R1 = 1000.0 * R1
      IF ( R0 .EQ. 0.0 ) R0 = 1.0

      ALPHA0 = DEGRAD * ALPHA0
      ALPHA1 = DEGRAD * ALPHA1
      DALPHA = ( ALPHA1 - ALPHA0 ) / ( NALPHA - 1 )

      BETA0 = DEGRAD * BETA0
      BETA1 = DEGRAD * BETA1
      DBETA = 0.0
      IF ( NBETA .NE. 1 ) DBETA = ( BETA1 - BETA0 ) / ( NBETA - 1 )

      XS = 0.0
      YS = 0.0

C     *** Compute 'optimal' beam constant ***

      CALL SSP(XS, YS, ZS, C0, CX, CY, CZ, CXX, CYY, CZZ, CXY, CXZ, CYZ)
      OMEGA = 2.0 * PI * FREQ

C     ------ Space filling in far field:
      IF ( OPTION( 2:2 ) .EQ. 'F' ) THEN
         WRITE( PRTFIL, * ) 'SPACE FILLING BEAMS'
         HWIDTH1 = 2.0 / ( ( OMEGA / C0 ) * DALPHA )
         EPSOPT1 = PI * FREQ * HWIDTH1 **2
         IF ( NBETA .GT. 1 ) HWIDTH2 = 2.0/( ( OMEGA/C0 ) * DBETA )
         EPSOPT2 = PI*FREQ*HWIDTH2**2
      ENDIF

C     ------ Minimum width at R1
      IF ( OPTION( 2:2 ) .EQ. 'M' ) THEN
         WRITE( PRTFIL, * ) 'MINIMUM WIDTH BEAMS'
         HWIDTH = SQRT( 2.0 * C0 * R1 / OMEGA )
         EPSOPT1 =  PI * FREQ * HWIDTH **2
         EPSOPT2 =  EPSOPT1
      ENDIF

      EPS1 = EPMULT * CI * EPSOPT1
      EPS2 = EPMULT * CI * EPSOPT2

      WRITE( PRTFIL, * ) 'HWIDTH1 = ', HWIDTH1, 'HWIDTH2 = ', HWIDTH2
      WRITE( PRTFIL, * ) 'EPSOPT1 = ', EPSOPT1, 'EPSOPT2 = ', EPSOPT2
      WRITE( PRTFIL, * ) 'EPMULT = ', EPMULT

      DELTAR = ( R1 - R0 ) / ( NR - 1 )
      DO 5000 IR = 1, NR
         R( IR ) = R0 + ( IR - 1 ) * DELTAR
 5000 CONTINUE

      DTHETA = 0.0
      IF ( NTHETA .NE. 1 ) DTHETA = ( THETA1-THETA0 ) / ( NTHETA-1 )
      DO 6000 ITHETA = 1, NTHETA
         THETA( ITHETA ) = THETA0 + ( ITHETA - 1 ) * DTHETA
 6000 CONTINUE

      NSD = 1
      CALL HEADER( TITLE, ZS, NSD, THETA, NTHETA,
     &   R, NR, FREQ, ATTEN, 'POLAR', XS, YS, ZR )

C     *** Construct a fan of beams and their influence ***

      CALL FAN( ALPHA0, DALPHA, NALPHA, BETA0, DBETA, NBETA, IBWIN,
     &   THETA, NTHETA, R, NR, DELTAR, DELTAS, NSTEPS, XS, YS, ZS, 
     &   ZR, EPS1, EPS2, OMEGA, C0, OPTION, OPT2, P )

C      CALL SYS$GETTIME( TEND, WTIME2 )
      WRITE( PRTFIL, 9500 ) TEND - TSTART
 9500 FORMAT( /, ' CPU Time = ', G15.3 )

      STOP
      END
C**********************************************************************C
      SUBROUTINE FAN( ALPHA0, DALPHA, NALPHA, BETA0, DBETA, NBETA,
     &  IBWIN, THETA, NTHETA, R, NR, DELTAR, DELTAS, NSTEPS, XS, YS, ZS,
     &   ZR, EPS1, EPS2, OMEGA, C0, OPTION, OPT2, P )

C     Computes a fan of beams from the source and writes the shade file

      IMPLICIT REAL ( A-H, O-Z )
      PARAMETER ( PI = 3.14159265 )

      REAL THETA( * ), R( * )
      COMPLEX P( NTHETA, * ), CONST, EPS1, EPS2
      CHARACTER TITLE*80, OPTION*8, OPT2*4

C     ------ Zero out field matrix
      DO ITHETA = 1, NTHETA
         DO IR = 1, NR
            P( ITHETA, IR ) = 0.0
         END DO
      END DO

C     *** Trace successive beams ***

      DO 8000 IBETA = 1, NBETA

         DO 7000 IALPHA = 1, NALPHA
            ALPHA = ALPHA0 + ( IALPHA - 1 ) * DALPHA
            BETA  = BETA0  + ( IBETA  - 1 ) * DBETA

C           ------ 3D or Nx2D computation
            IF ( OPT2 .EQ. '3D' ) THEN
               CALL TRACE( XS, YS, ZS, ALPHA, BETA, EPS1, EPS2, 
     &            OMEGA, DELTAS, NSTEPS, OPTION )
            ELSE
               CALL TRACE2D( XS, YS, ZS, ALPHA, BETA, EPS1, EPS2, 
     &            OMEGA, DELTAS, NSTEPS, OPTION )
            ENDIF

C           ------ Optionally dump rays to a disk file
            IF ( OPTION( 1:1 ) .EQ. 'R' .AND. IBETA .EQ. 1 ) THEN
C               WRITE( RAYFIL ) ( XV( I ), ZV( I ), XIV( I ), ZETAV( I ), 
C     &            FV( I ), GV( I ), I = 1, NSTEPS )
            ENDIF

C           *** Compute contribution to field ***

            IF ( OPT2 .EQ. '3D' ) THEN
               CALL INFLU( ALPHA, IBWIN, ZR, THETA, NTHETA, R, NR, 
     &            DELTAR, P, NSTEPS, OMEGA )
            ELSE
               CALL INFLU2D( ALPHA, IBWIN, ZR, BETA, IBETA, NTHETA,
     &            R, NR, DELTAR, P, NSTEPS, OMEGA, EPS1, OPTION )
            ENDIF

 7000    CONTINUE
 8000 CONTINUE

C     *** Scale the complex pressure field ***

      IF ( OPT2 .EQ. '3D' ) THEN
         CONST = SQRT( EPS1 * EPS2 ) * OMEGA * DBETA * DALPHA/
     &      ( 2.0 * PI* ( SQRT( C0 ) ) **3 )
      ELSE
         CONST = SQRT( ABS( EPS1 ) * OMEGA * DALPHA**2/
     &      ( 2.0 * PI * C0 **2 ) )
      ENDIF

      DO ITHETA = 1, NTHETA
         DO IR = 1, NR
            IF ( OPT2 .EQ. '3D' ) THEN
               P( ITHETA, IR ) = CONST * P( ITHETA, IR )
            ELSE
               P( ITHETA, IR ) = CONST * P( ITHETA, IR )/SQRT( R( IR ) )
            ENDIF
         END DO
      END DO

C     *** Write the field to the shade file ***

      IREC = 6
      CALL WRTFLD( P, NTHETA, NR, IREC )

      RETURN
      END
C**********************************************************************C
      SUBROUTINE TRACE( XS, YS, ZS, ALPHA, BETA, EPS1, EPS2,
     &   OMEGA, DELTAS, NSTEPS, OPTION )

C     Traces the beam corresponding to a particular take off angle

      IMPLICIT REAL ( A-H, O-Z )
      PARAMETER ( MAXN = 10003 , PI = 3.14159265 )

      COMMON /BEAM/XV, YV, ZV, XIV, ETAV, ZETAV, FV, GV, HV, DETQV,
     &   PHIV, TAUV, CV

      REAL    XV( MAXN ),    YV( MAXN ),    ZV( MAXN ),
     &       XIV( MAXN ),  ETAV( MAXN ), ZETAV( MAXN ),
     &      PHIV( MAXN ),  TAUV( MAXN ),    CV( MAXN )
      COMPLEX FV( MAXN ),    GV( MAXN ),    HV( MAXN ),
     &     DETPV( MAXN ), DETQV( MAXN ),
     &   F, G, H, DETP, DETQ, EPS1, EPS2
      CHARACTER OPTION*80

C     *** Initial conditions ***

      XV( 1 ) = XS
      YV( 1 ) = YS
      ZV( 1 ) = ZS

      CALL SSP( XV( 1 ), YV( 1 ), ZV( 1 ), C, CX, CY, CZ,
     &   CXX, CYY, CZZ, CXY, CXZ, CYZ )

      XIV(   1 ) = COS( ALPHA ) * COS( BETA ) / C
      ETAV(  1 ) = COS( ALPHA ) * SIN( BETA ) / C
      ZETAV( 1 ) = SIN( ALPHA ) / C

      FV( 1 ) = EPS2
      GV( 1 ) = EPS1
      HV( 1 ) = 0.0

      DETPV( 1 ) = 1.0
      DETQV( 1 ) = EPS1 * EPS2

      PHIV( 1 ) = 0.0
      TAUV( 1 ) = 0.0

C     *** Trace the beam ***

      I = 0
      DO 3000 ISTEP = 1, NSTEPS
         I = I + 1

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

         CALL SSP( XV( I ), YV( I ), ZV( I ), C, CX, CY, CZ,
     &      CXX, CYY, CZZ, CXY, CXZ, CYZ )
         CV( I ) = C
         C2 = C * C

C        ------ Compute ray normals E1 and E2
         CALL RAYNOR( XIV( I ), ETAV( I ), ZETAV( I ), PHIV( I ), C, 
     &      E1X, E1Y, E1Z, E2X, E2Y, E2Z )

C        ------ Compute second partials of c along ray normals
         CALL VCALC( CXX, CXY, CXZ, CYY, CYZ, CZZ, E1X, E1Y, E1Z,
     &      E2X, E2Y, E2Z, CNN, CNM, CMM )

         STEP = DELTAS

C        ------ If crossing surface reduce step to land at surface
         Z = ZV( I ) + STEP * C * ZETAV( I )
         IF ( Z .LT. 0.0 ) STEP = -ZV( I ) / ( ZETAV( I ) * C )
         HSTEP = 0.5 * STEP

         X = XV( I ) + HSTEP * C *   XIV( I )
         Y = YV( I ) + HSTEP * C *  ETAV( I )
         Z = ZV( I ) + HSTEP * C * ZETAV( I )

         XI   = XIV( I )   - HSTEP * CX / C2
         ETA  = ETAV( I )  - HSTEP * CY / C2
         ZETA = ZETAV( I ) - HSTEP * CZ / C2

         F = FV( I ) + HSTEP * ( C*DETPV( I ) - ( CNN/C2 )*DETQV( I ) )
         G = GV( I ) + HSTEP * ( C*DETPV( I ) - ( CMM/C2 )*DETQV( I ) )
         H = HV( I ) + HSTEP * ( - ( CNM/C2 )*DETQV( I ) )

         DETP = DETPV( I ) + HSTEP / C2 * 
     &      ( -CMM*FV( I ) - CNN*GV( I ) + 2.0*CNM*HV( I ) )
         DETQ = DETQV( I ) + HSTEP * C * ( FV( I ) + GV( I ) )

         PHI = PHIV( I ) + HSTEP *( 1.0 / C ) * ZETAV( I )*
     &      ( ETAV( I )*CX-XIV( I )*CY )/( XIV( I )**2+ETAV( I )**2 )

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

         CALL SSP( X, Y, Z, C, CX, CY, CZ, CXX, CYY, CZZ, CXY, CXZ, CYZ)

         C2 = C * C

C        ------ Compute ray normals E1 and E2
         CALL RAYNOR( XI, ETA, ZETA, PHI, C, 
     &      E1X, E1Y, E1Z, E2X, E2Y, E2Z )

C        ------ Compute second partials of c along ray normals
         CALL VCALC( CXX, CXY, CXZ, CYY, CYZ, CZZ, E1X, E1Y, E1Z,
     &      E2X, E2Y, E2Z, CNN, CNM, CMM )

C        ------ If crossing surface reduce step to land at surface
         Z = ZV( I ) + STEP*C*ZETAV( I )
         IF ( Z .LE. 0.0 ) STEP = -1.00000001*ZV( I )/( ZETAV( I )*C )

         XV( I+1 ) = XV( I ) + STEP * C * XI
         YV( I+1 ) = YV( I ) + STEP * C * ETA
         ZV( I+1 ) = ZV( I ) + STEP * C * ZETA

         XIV( I+1 )   = XIV( I )   - STEP * CX / C2
         ETAV( I+1 )  = ETAV( I )  - STEP * CY / C2
         ZETAV( I+1 ) = ZETAV( I ) - STEP * CZ / C2

         FV( I+1 ) = FV( I ) + STEP * ( C*DETP - ( CNN/C2 )*DETQ )
         GV( I+1 ) = GV( I ) + STEP * ( C*DETP - ( CMM/C2 )*DETQ )
         HV( I+1 ) = HV( I ) + STEP * ( - ( CNM/C2 )*DETQ )

         DETPV( I+1 ) = DETPV( I ) + STEP / C2 * 
     &           ( -CMM*F - CNN*G + 2.0*CNM*H )
         DETQV( I+1 ) = DETQV( I ) + STEP * C * ( F + G )

         PHIV( I+1 ) = PHIV( I ) + STEP * ( 1.0 / C )*
     &      ZETA * ( ETA * CX - XI * CY ) / ( XI**2 + ETA**2 )

         TAUV( I+1 ) = TAUV( I ) + STEP / C

C        *** If reflecting off surface need jump conditions ***

         IF ( ZV( I+1 ) .LE. 0.0 ) THEN
            I = I + 1

               XV( I+1 ) =  XV( I )
               YV( I+1 ) =  YV( I )
               ZV( I+1 ) = -ZV( I )

              XIV( I+1 ) =    XIV( I )
             ETAV( I+1 ) =   ETAV( I )
            ZETAV( I+1 ) = -ZETAV( I )

            CALL SSP( XV( I+1 ), YV( I+1 ), ZV( I+1 ), C, CX, CY, CZ, 
     &         CXX, CYY, CZZ, CXY, CXZ, CYZ )
            CV( I ) = C

C           ------ Following is suspect
            RN1 = 0.0
            RN2 = 0.0
C            RN1 = -2.0*XIV(  I+1 )**2 / ZETAV( I+1 ) * CZ / C
C            RN2 = -2.0*ETAV( I+1 )**2 / ZETAV( I+1 ) * CZ / C

            IF ( OPTION( 3:3 ) .EQ. 'D' ) THEN
               RN1 = 2.0 * RN1
               RN2 = 2.0 * RN2
            ENDIF

            FV( I+1 ) = FV( I ) + DETQV( I ) * RN1
            GV( I+1 ) = GV( I ) + DETQV( I ) * RN1
            HV( I+1 ) = HV( I ) + DETQV( I ) * RN2

            DETPV( I+1 ) = DETPV( I )
            DETQV( I+1 ) = DETQV( I )

C           ------ What happens to torsion?
            PHIV( I+1 ) = PHIV( I )
            TAUV( I+1 ) = TAUV( I ) - PI / OMEGA
         ENDIF

 3000 CONTINUE

      RETURN
      END
C**********************************************************************C
      SUBROUTINE INFLU( ALPHA, IBWIN, ZR, THETA, NTHETA, R, NR, DELTAR,
     &   P, NSTEPS, OMEGA )

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

      IMPLICIT REAL ( A-H, O-Z )
      COMPLEX CI
      PARAMETER ( MAXN = 10003 , PI = 3.14159265, CI = ( 0.0, 1.0 ), 
     &   DEGRAD = PI / 180.0 )
      COMMON /BEAM/XV, YV, ZV, XIV, ETAV, ZETAV, FV, GV, HV, DETQV,
     &   PHIV, TAUV, CV

      INTEGER KMAHV( MAXN )

      REAL XV( MAXN ),   YV( MAXN ),    ZV( MAXN ),
     &    XIV( MAXN ), ETAV( MAXN ), ZETAV( MAXN ),
     &   PHIV( MAXN ), TAUV( MAXN ),    CV( MAXN ), R( * ),
     &   THETA( * ), NA, NB, NSQ, LSQ1,
     &               MA, MB, MSQ, LSQ2, NM

      COMPLEX FV( MAXN ), GV( MAXN ), HV( MAXN ), DETQV( MAXN ),
     &   F, G, H, DETQ, CONTRI, P( NTHETA, * )

      IBWIN2 = IBWIN

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     *** Begin by forming KMAH index ***

      KMAHV( 1 ) = 1
      DO 1000 I = 2, NSTEPS
         KMAHV( I ) = KMAHV( I - 1 )

         IF ( REAL( DETQV( I ) ) .LT. 0.0 ) THEN
            QL = AIMAG( DETQV( I - 1 ) )
            QR = AIMAG( DETQV( I ) )
            IF ( ( QL .LT. 0.0 .AND. QR .GE. 0.0 ) .OR.
     &           ( QL .GT. 0.0 .AND. QR .LE. 0.0 ) )
     &         KMAHV( I ) = -KMAHV( I )
         ENDIF

 1000 CONTINUE

C     *** Loop over radials of receiver line ***

      DO 9500 ITHETA = 1, NTHETA
      TH = DEGRAD * THETA( ITHETA )
C     DO 9000 IMAGE = 1, 1
         IMAGE = 1
         IR1   = 9999
         DO 5000 I = 1, NSTEPS

C           ------ Compute ray normals E1 and E2
            CALL RAYNOR( XIV( I ), ETAV( I ), ZETAV( I ),
     &         PHIV( I ), CV( I ), E1X, E1Y, E1Z, E2X, E2Y, E2Z )

C           ------ Compute coordinates of intercept: NB, MB, RB
            ZT = ZV( I ) - ZR
            CT = COS( TH )
            ST = SIN( TH )

            DELTA = -CT * ( E1Y * E2Z - E1Z * E2Y ) +
     &               ST * ( E1X * E2Z - E1Z * E2X )

C           ------ Check for ray normal || radial of rcvr line
            IF ( ABS( DELTA ) .LT. 1.0E-10 ) GOTO  9500
            NB =   ( -CT * ( YV( I ) * E2Z - ZT * E2Y ) +
     &                ST * ( XV( I ) * E2Z - ZT * E2X ) ) / DELTA
            MB =   ( -CT * ( E1Y * ZT - E1Z * YV( I ) ) +
     &                ST * ( E1X * ZT - E1Z * XV( I ) ) ) / DELTA
            RB = ( -XV( I )* ( E1Y * E2Z - E1Z * E2Y )
     &             +YV( I )* ( E1X * E2Z - E1Z * E2X )
     &               -ZT   * ( E1X * E2Y - E1Y * E2X ) ) / DELTA

C           ------ Compute nearest rcvr before normal ( avoiding intovfl )
            IF ( RB .GT. R( NR ) + DELTAR ) THEN
               IR2 = NR
            ELSE IF ( RB .LT. R( 1 ) ) THEN
               IR2 = 0
            ELSE
               IR2 = IFIX( REAL( ( RB - R( 1 ) ) / DELTAR ) ) + 1
            ENDIF

C           *** Compute contributions to bracketted receivers ***

            IF ( (IR2 .GT. IR1) .AND. (XV( I ) .NE. XV( I-1 )) ) THEN
               DO 2000 IR = IR1 + 1, IR2
                  W = ( R( IR )-RA ) / ( RB-RA )
                  NSQ = ( NA + W * ( NB-NA ) )**2
                  MSQ = ( MA + W * ( MB-MA ) )**2

C                 ------ Within beam window?
                  LSQ1 = -2.0/( OMEGA * AIMAG( FV( I )/DETQV( I ) ) )
                  LSQ2 = -2.0/( OMEGA * AIMAG( GV( I )/DETQV( I ) ) )

                  IF ( NSQ .LT. IBWIN2 * LSQ1 .AND.
     &                 MSQ .LT. IBWIN2 * LSQ2  ) THEN

                     F = FV( I-1 ) + W * ( FV( I ) - FV( I-1 ) )
                     G = GV( I-1 ) + W * ( GV( I ) - GV( I-1 ) )
                     H = HV( I-1 ) + W * ( HV( I ) - HV( I-1 ) )

                     DETQ = DETQV( I-1 ) + W*( DETQV( I )-DETQV( I-1 ) )
                     C    = CV(    I-1 ) + W*(    CV( I )-CV(    I-1 ) )
                     TAU  = TAUV(  I-1 ) + W*(  TAUV( I )-TAUV(  I-1 ) )
                     KMAH = KMAHV( I-1 )

                     IF ( REAL( DETQ ) .LT. 0.0 ) THEN
                        QL = AIMAG( DETQV( I - 1 ) )
                        QR = AIMAG( DETQ )
                        IF ( ( QL .LT. 0.0 .AND. QR .GE. 0.0 ) .OR.
     &                       ( QL .GT. 0.0 .AND. QR .LE. 0.0 ) )
     &                     KMAH = -KMAH
                     ENDIF

                     NM = SQRT( NSQ * MSQ )
                     CONTRI = COS( ALPHA )*SQRT( C / DETQ ) *
     &                  EXP( -CI * OMEGA *( TAU +
     &                     0.5*( F*NSQ + 2.0*H*NM + G*MSQ ) / DETQ ) )

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

                     P( ITHETA, IR ) = P( ITHETA, IR ) + CONTRI
                  ENDIF
 2000          CONTINUE
            ENDIF

            NA = NB
            MA = MB
            RA = RB

            IR1 = IR2
C           ------ Next step
 5000    CONTINUE
C      ------ Next image
C 9000 CONTINUE
C     ------ Next radial
 9500 CONTINUE

      RETURN
      END
C**********************************************************************C
      SUBROUTINE RAYNOR( XI, ETA, ZETA, PHI, C, 
     &   E1X, E1Y, E1Z, E2X, E2Y, E2Z )

C     Computes the ray normals, E1 and E2
C     following error involving XIV existed until April 95:
C     E1Y = ( C * ETA * ZETA * COS( PHI ) -
C    &                        XIV * SIN( PHI ) ) / RL

      IMPLICIT REAL ( A-H, O-Z )

      RL = SQRT( XI**2 + ETA**2 )

C     ------ E1
      E1X = ( C * XI  * ZETA * COS( PHI ) +
     &                        ETA * SIN( PHI ) ) / RL
      E1Y = ( C * ETA * ZETA * COS( PHI ) -
     &                        XI  * SIN( PHI ) ) / RL
      E1Z = -C * RL * COS( PHI )

C     ------ E2
      E2X = ( C * XI  * ZETA * SIN( PHI ) -
     &                        ETA * COS( PHI ) ) / RL
      E2Y = ( C * ETA * ZETA * SIN( PHI ) +
     &                        XI  * COS( PHI ) ) / RL
      E2Z = -C * RL * SIN( PHI )

      RETURN
      END
C**********************************************************************C
      SUBROUTINE VCALC( CXX, CXY, CXZ, CYY, CYZ, CZZ, E1X, E1Y, E1Z,
     &   E2X, E2Y, E2Z, CNN, CNM, CMM )

C     Computes the second partials of c along ray normals

      IMPLICIT REAL ( A-H, O-Z )

      CNN = CXX * E1X**2 + CYY * E1Y**2 + CZZ * E1Z**2 +
     &   2.0*CXY * E1X*E1Y + 2.0*CXZ * E1X*E1Z + 2.0*CYZ * E1Y*E1Z

      CNM = CXX * E1X*E2X + CYY * E1Y*E2Y + CZZ * E1Z*E2Z +
     &   CXY * ( E1X*E2Y + E2X*E1Y ) + CXZ * ( E1X*E2Z + E2X*E1Z ) + 
     &   CYZ * ( E1Y*E2Z+E2Y*E1Z )

      CMM = CXX * E2X**2 + CYY * E2Y**2 + CZZ * E2Z**2 +
     &   2.0*CXY * E2X*E2Y + 2.0*CXZ * E2X*E2Z + 2.0*CYZ * E2Y*E2Z

      RETURN
      END
C**********************************************************************C
      SUBROUTINE TRACE2D( XS, YS, ZS, ALPHA, BETA, EPS1, EPS2, 
     &   OMEGA, DELTAS, NSTEPS, OPTION )

C     Traces the beam corresponding to a particular take off angle

      IMPLICIT REAL ( A-H, O-Z )

      PARAMETER ( MAXN = 10003 , PI = 3.14159265 )

      COMMON /BEAM/ XV, YV, ZV, XIV, ETAV, ZETAV, FV, GV, HV, DETQV,
     &   PHIV, TAUV, CV
      REAL XV( MAXN ),   YV( MAXN ),    ZV( MAXN ),
     &    XIV( MAXN ), ETAV( MAXN ), ZETAV( MAXN ),
     &   PHIV( MAXN ), TAUV( MAXN ),    CV( MAXN )
      COMPLEX FV( MAXN ), GV( MAXN ), HV( MAXN ),
     &   DETQV( MAXN ), EPS1, EPS2
      CHARACTER OPTION*80

      COMMON /BEAM2D/ RV, RHOV, P1V, Q1V, P2V, Q2V
      REAL*4 RV( MAXN ), RHOV( MAXN ), 
     &   P1V( MAXN ), Q1V( MAXN ), P2V( MAXN ), Q2V( MAXN )

C     *** Initial conditions ***

      DEPTH = 1.0E6
      RS = SQRT( XS**2 + YS**2 )
      RV( 1 ) = RS
      ZV( 1 ) = ZS

      CALL SSP2D( RV( 1 ), ZV( 1 ), C, CR, CZ, CRR, CRZ, CZZ, BETA )

      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

C     *** Trace the beam ***

      I = 0
      DO 3000 ISTEP = 1, NSTEPS
         I = I + 1

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

         CALL SSP2D( RV( I ), ZV( I ), C, CR, CZ, CRR, CRZ, CZZ, BETA )

         CV( I ) = C
         C2 = C * C
         CNN = C2*( CRR*ZETAV( I )**2+2.0*CRZ*RHOV( I )*ZETAV( I )
     &             +CZZ*RHOV( I )**2 )
         STEP = DELTAS

C        ------ If crossing surface or bottom, reduce step to hit interface
         Z = ZV( I ) + STEP * C * ZETAV( I )
         IF ( Z .LT. 0.0   ) STEP =        -ZV( I )   / ( ZETAV( I )*C )
         IF ( Z .GT. DEPTH ) STEP = ( DEPTH-ZV( I ) ) / ( ZETAV( I )*C )

         R = RV( I ) + 0.5*STEP * C * RHOV(  I )
         Z = ZV( I ) + 0.5*STEP * C * ZETAV( I )

         RHO  = RHOV( I )  - 0.5 * STEP * CR/C2
         ZETA = ZETAV( I ) - 0.5 * STEP * CZ/C2

         P1 = P1V( I ) - 0.5*STEP * CNN / C2 * Q1V( I )
         Q1 = Q1V( I ) + 0.5*STEP * C * P1V( I )
         P2 = P2V( I ) - 0.5*STEP * CNN / C2 * Q2V( I )
         Q2 = Q2V( I ) + 0.5*STEP * C * P2V( I )

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

         CALL SSP2D( R, Z, C, CR, CZ, CRR, CRZ, CZZ, BETA )

         C2 = C * C
         CNN = CRR * ( ZETA * C )**2 +
     &      2.0 * CRZ * ( RHO * C * ZETA * C ) + CZZ * ( RHO * C )**2

C        ------ If crossing surface or bottom, reduce step to hit interface
         Z = ZV( I ) + STEP * C * ZETAV( I )
         IF ( Z .LE. 0.0 ) STEP = -1.00000001*ZV( I )/( ZETAV( I )*C )
         IF ( Z .GE. DEPTH )
     &      STEP = 1.0000001 * ( DEPTH - ZV( I ) ) / ( ZETAV( I ) * C )

         RV( I+1 ) = RV( I ) + STEP * C * RHO
         ZV( I+1 ) = ZV( I ) + STEP * C * ZETA

         RHOV(  I+1 ) = RHOV(  I ) - STEP * CR / C2
         ZETAV( I+1 ) = ZETAV( I ) - STEP * CZ / C2

         P1V( I+1 ) = P1V( I ) - STEP * CNN / C2 * Q1
         Q1V( I+1 ) = Q1V( I ) + STEP * C * P1
         P2V( I+1 ) = P2V( I ) - STEP * CNN / C2 * Q2
         Q2V( I+1 ) = Q2V( I ) + STEP * C * P2

         TAUV( I+1 ) = TAUV( I ) + STEP / C

C        *** If reflecting off surface, need jump conditions ***

         IF ( ZV( I+1 ) .LE. 0.0 ) THEN
            I = I + 1

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

             RHOV( I+1 ) =   RHOV( I )
            ZETAV( I+1 ) = -ZETAV( I )

            CALL SSP2D( RV( I+1 ), ZV( I+1 ), C, CR, CZ,
     &         CRR, CRZ, CZZ, BETA )

            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 ( OPTION( 3:3 ) .EQ. 'D' ) RN = 2.0 * RN

            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 )
            TAUV( I+1 ) = TAUV( I ) - PI / OMEGA
         ENDIF

C        *** If reflecting off bottom, need jump conditions ***

         IF ( ZV( I+1 ) .GT. DEPTH ) THEN
            I = I + 1

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

            ZETAV( I+1 ) = -ZETAV( I )
            RHOV(  I+1 ) =  RHOV(  I )

            CALL SSP2D( RV( I+1 ), ZV( I+1 ), C, CR, CZ,
     &                                   CRR, CRZ, CZZ, BETA )
            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 ( OPTION( 3:3 ) .EQ. 'D' ) RN = 2.0 * RN

            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 )

            TAUV( I+1 ) = TAUV( I )
         ENDIF
 3000 CONTINUE

      RETURN
      END
C**********************************************************************C
      SUBROUTINE INFLU2D( ALPHA, IBWIN, ZR, BETA, IBETA, NTHETA,
     &   R, NR, DELTAR, PRESS, NSTEPS, OMEGA, EPS, OPTION )

C     Computes the beam influence, i.e. 
C     the contribution of a single beam to the COMPLEX*8 pressure

      IMPLICIT REAL ( A-H, O-Z )
      COMPLEX CI
      PARAMETER ( MAXN = 10003, CI = ( 0.0, 1.0 ), NRD = 1 )
      COMMON /BEAM/XV, YV, ZV, XIV, ETAV, ZETAV, FV, GV, HV, DETQV,
     &   PHIV, TAUV, CV

      REAL XV( MAXN ),   YV( MAXN ),    ZV( MAXN ),
     &    XIV( MAXN ), ETAV( MAXN ), ZETAV( MAXN ),
     &   PHIV( MAXN ), TAUV( MAXN ), CV( MAXN ), R( * )

      COMPLEX FV( MAXN ), GV( MAXN ), HV( MAXN ), DETQV( MAXN ),
     &   PRESS( NTHETA, * )

      COMMON /BEAM2D/ RV, RHOV, P1V, Q1V, P2V, Q2V
      INTEGER KMAHV( MAXN )
      REAL*4 RV( MAXN ), RHOV( MAXN ), 
     &   P1V( MAXN ), Q1V( MAXN ), P2V( MAXN ), Q2V( MAXN ),
     &   NA, NB, NSQ, RD( NRD )
      COMPLEX*8 PV( MAXN ), QV( MAXN ), P, Q, EPS, CONTRI
      CHARACTER OPTION*80

      IBWIN2  = IBWIN**2
      DEPTH   = 1.0E6
      RD( 1 ) = ZR

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     *** Begin by forming ( P, Q ) and KMAH index ***

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

      KMAHV( 1 ) = 1

      DO 1000 I = 2, NSTEPS
         PV( I ) = P1V( I ) + EPS * P2V( I )
         QV( I ) = Q1V( I ) + EPS * Q2V( I )
         KMAHV( I ) = KMAHV( I-1 )

         IF ( REAL( QV( I ) ) .LT. 0.0 ) THEN
            Q1 = AIMAG( QV( I-1 ) )
            Q2 = AIMAG( QV( I ) )
            IF ( ( Q1 .LT. 0.0 .AND. Q2 .GE. 0.0 ) .OR.
     &           ( Q1 .GT. 0.0 .AND. Q2 .LE. 0.0 ) )
     &                       KMAHV( I ) = -KMAHV( I )
         ENDIF
 1000 CONTINUE

C     *** Loop over receivers ***

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

C     Loop over images

      IF ( OPTION( 4:4 ) .EQ. 'R' ) NIMAGE = 1
      IF ( OPTION( 4:4 ) .EQ. 'S' ) NIMAGE = 2
      IF ( OPTION( 4:4 ) .EQ. 'B' ) NIMAGE = 3

      DO 8000 IMAGE = 1, NIMAGE
         IR1 = 9999

C        *** Loop over steps ***

         DO 5000 I = 1, NSTEPS

C           Compute ray-centered coordinates

            ZNV = -RHOV( I ) * CV( I )

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

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

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

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

            IR2 = MIN( IFIX( REAL( ( RB-R( 1 ) )/DELTAR ) ) + 1, NR )
            IF ( (IR2 .GT. IR1) .AND. (RV( I ) .NE. RV( I-1 )) ) 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 ) )
                  NSQ = ( NA + W * ( NB-NA ) )**2

C                 ------ Within beam window?
                  IF ( -0.5*OMEGA*AIMAG( P/Q )*NSQ .LT. IBWIN2 ) THEN
                     C   =   CV( I-1 ) + W * (   CV( I ) -   CV( I-1 ) )
                     TAU = TAUV( I-1 ) + W * ( TAUV( I ) - TAUV( I-1 ) )

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

C                    ------ Get correct branch of sqrt
                     KMAH = KMAHV( I - 1 )
                     IF ( REAL( Q ) .LT. 0.0 ) THEN
                        Q1 = AIMAG( QV( I - 1 ) )
                        Q2 = AIMAG( Q )
                        IF ( ( Q1 .LT. 0.0 .AND. Q2 .GE. 0.0 ) .OR.
     &                       ( Q1 .GT. 0.0 .AND. Q2 .LE. 0.0 ) )
     &                                KMAH = -KMAH
                     ENDIF

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

                     PRESS( IBETA, IR ) = PRESS( IBETA, IR ) + CONTRI

                  ENDIF
 2000          CONTINUE
            ENDIF

            RA  = RB
            NA  = NB
            IR1 = MAX( IR2, 0 )
C           ------ Next step
 5000    CONTINUE
C      ------ Next image
 8000 CONTINUE
C     ------ Next receiver
 9000 CONTINUE

      RETURN
      END
C**********************************************************************C
      SUBROUTINE SSP2D( R, Z, C, CR, CZ, CRR, CRZ, CZZ, BETA )

C     Converts cartesian gradients to polar

      X = R * COS( BETA )
      Y = R * SIN( BETA )

      CALL SSP( X, Y, Z, C, CX, CY, CZ, CXX, CYY, CZZ, CXY, CXZ, CYZ )

      CR  = COS( BETA ) * CX  + SIN( BETA ) * CY
      CRZ = COS( BETA ) * CXZ + SIN( BETA ) * CYZ

      CRR =       CXX * ( COS( BETA ) )**2
     &    + 2.0 * CXY *   COS( BETA ) * SIN( BETA )
     &    +       CYY * ( SIN( BETA ) )**2

      RETURN
      END
