      PROGRAM PLOTIRC

C     Plots the internal reflection coefficient

      INTEGER PLPFIL
      PARAMETER( PLPFIL = 1, PI = 3.141592 )

      INTEGER ITAB( 5000 )
      REAL C( 5000 ), Z( 5000 ), PHIR( 5000 )
      REAL*8 XTAB( 5000 )
      COMPLEX*8 CK( 5000 )
      COMPLEX*16 FTAB( 5000 ), GTAB( 5000 )
      CHARACTER TITLE*80

C     User parameters

      CMIN = 1495.0
      CMAX = 1515.0
      CINC = 0.0

      OPEN( FILE = 'PLPFIL', UNIT = PLPFIL, STATUS = 'OLD' )
      READ( PLPFIL, * ) CMIN, CMAX

C     &   'Phase Speed (m/s)$', 100, 'Z$', 100, CAXL,  ZAXL )

C     *** Plot positions of eigenvalues ***

      CALL GETMOD( 1, ' ', MAXM, RD, 1, 'N',
     &   CK, PHIR, M, FREQ, TITLE )

      OMEGA = 2.0 * PI * FREQ

      DO MODE = 1, M
         C( MODE ) = OMEGA / REAL( CK( MODE ) )
      END DO

C     *** Plot the two IRC files ***

      DO IRCFIL = 10, 11
      READ( IRCFIL ) NKPTS, (XTAB( IK ), FTAB( IK ), GTAB( IK ), 
     &                                      ITAB(IK), IK = 1, NKPTS)

C     Search for maximum ITAB

      ITMAX = -99999
      ITMIN =  99999

      DO IK = 1, NKPTS
         ITAB( IK ) = 0
         ITMAX = MAX ( ITAB( IK ),  ITMAX )
         ITMIN = MIN ( ITAB( IK ),  ITMIN )
      END DO

C     Convert to impedance and plot

      ISTART = 1

      DO IK = 1, NKPTS

         C(IK) = OMEGA / SQRT( XTAB( IK ) )
         IF ( ABS( FTAB( IK ) ) .LT. 1.0E-25 ) THEN
            Z( IK ) = 0.0
         ELSE
            Z( IK ) = REAL( GTAB( IK ) / FTAB( IK ) )
         ENDIF

         IF ( IRCFIL .EQ. 11 ) Z( IK ) = -Z( IK )

         IF ( REAL( FTAB( IK ) ) * REAL( FTAB( IK + 1 ) ) .LT. 0.0 .OR.
     &                                   IK .GE. NKPTS ) THEN

            IF ( C( ISTART ) .GT. CMIN .AND. C( IK ) .LT. CMAX ) THEN
               DO I = ISTART, IK - ISTART + 1
                  WRITE( *, * ) C( I ), Z( I )
               END DO
            ENDIF

            ISTART = IK + 1     ! RESET LEFT SIDE OF WINDOW
         ENDIF

      END DO   ! next k-point

      WRITE( *, * ) '&'

      END DO   ! next IRCFIL

C     *** Write data ***

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

      WRITE( *, 2000 ) NINT( FREQ )
 2000 FORMAT(  '@ subtitle "',
     &         'F  =', I4, ' Hz"' )
 

      STOP
      END
