      PROGRAM PLOTRTH

C     Plots the reflection coefficient

      INTEGER IRCFIL, PLPFIL, PRTFIL
      COMPLEX CI
      PARAMETER( PLPFIL = 1, PRTFIL = 6, IRCFIL = 10,
     &   PI = 3.141592, CI = ( 0.0, 1.0 ),
     &   RADDEG = 180.0 / PI, MXNKPTS = 30000 )

      INTEGER ITAB( MXNKPTS )
      REAL R( MXNKPTS ), TH( MXNKPTS ), PHASE( MXNKPTS ), K0, KX, KZ
      REAL*8 XTAB( MXNKPTS )
      COMPLEX Z, RCMPLX
      COMPLEX*16 FTAB( MXNKPTS ), GTAB( MXNKPTS )
      CHARACTER TITLE*80, TITLEY*80, OPT*5

C     Note: we assume unit density in the upper medium

C     Read user parameters

      OPEN( FILE = 'PLPFIL', UNIT = PLPFIL, STATUS = 'OLD' )
      READ( PLPFIL, * ) OPT
      READ( PLPFIL, * ) C0

      THMIN =  0.0
      THMAX = 90.0

      READ( PLPFIL, * ) THMIN, THMAX

C     Initialize plotting

      TITLEY = 'Refl. coef.$'
      IF ( OPT(1:1) .EQ. 'D' ) TITLEY = 'Refl. loss (dB)$'

C     *** Read IRCFIL and convert to R( theta ) ***

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

      OMEGA = 2.0 * PI * FREQ
      K0 = OMEGA / C0

C     *** Convert to impedance and plot ***

      DO IK = 1, NKPTS

C        --- compute the impedance Z = G / F
         IF ( ABS( FTAB( IK ) ) .LT. 1.0E-25 ) THEN
            Z = 0.0
         ELSE
            Z = GTAB( IK ) / FTAB( IK )
         ENDIF

         C  = OMEGA / SQRT( XTAB( IK ) )
         KX = SQRT( XTAB( IK ) )

C        --- compute the vertical wavenumber kz
         IF ( K0 ** 2 .LT. KX ** 2 ) THEN
            KZ = 0.0
         ELSE
            KZ = SQRT( K0 ** 2 - KX ** 2 )
         ENDIF

C        --- convert wavenumbers (kx, kz) to angles
C         IF ( KZ .GT. 1.0E6 * KX ) THEN
C            TH( IK ) = 0.0
C         ELSE
            TH( IK ) = RADDEG * ATAN( KZ / KX )
C         ENDIF

C        --- convert impedance Z to reflection coef. R

         RCMPLX =  - ( 1.0 - CI * KZ * Z )
     &             / ( 1.0 + CI * KZ * Z )

         R(     IK ) = ABS( RCMPLX )
         PHASE( IK ) = ATAN2( AIMAG( RCMPLX ), REAL( RCMPLX ) )

C        --- dB, degrees for units?
         IF ( OPT(1:1) .EQ. 'D' ) R( IK ) = -20.0 * ALOG10( R( IK ) )
         IF ( OPT(2:2) .EQ. 'D' ) PHASE( IK ) = RADDEG * PHASE( IK )

      END DO   ! next k-point

C     *** Write data ***

      LAST = INDEX( TITLE, '    ' ) - 1
      WRITE( PRTFIL, 1000 ) TITLE(1:LAST)
 1000 FORMAT( '@ title "', A, '"' )

      WRITE( PRTFIL, 2000 ) NINT( FREQ )
 2000 FORMAT(  '@ subtitle "',
     &         'F  =', I4, ' Hz"' )
 
      WRITE( *, 5000 )
 5000 FORMAT( '@WITH G1', /,
     &        '@G1 ON', /,
     &        '@TYPE xy' )

      DO IK = 1, NKPTS
         WRITE( PRTFIL, * ) TH( IK ), R( IK )
      END DO
      WRITE( PRTFIL, * ) '&'

      WRITE( *, 6000 )
 6000 FORMAT( '@WITH G0', /,
     &        '@G0 ON', /,
     &        '@TYPE xy' )
 
      DO IK = 1, NKPTS
         WRITE( PRTFIL, * ) TH( IK ), PHASE( IK )
      END DO

C     *** produce a reflection coefficient file in the BELLHOP format

      IF( OPT(3:3) .EQ. 'F' ) THEN
         WRITE( PRTFIL, * ) NKPTS
         DO IK = NKPTS, 1, -1
            IF ( PHASE( IK ) .LT. 0.0 ) PHASE( IK ) = PHASE( IK ) + 360.0
            WRITE( PRTFIL, * ) TH( IK ), R( IK ), PHASE( IK )
         END DO
      END IF

      STOP
      END
