      PROGRAM PLOTSLIC

C     Plot the propagation loss along a specified track

      INTEGER PLPFIL, PRTFIL, SHDFIL
      PARAMETER ( PLPFIL = 1, PRTFIL = 6, 
     &   MAXNSD = 51,  MAXNRD = 4096, MAXNR = 5000 )

      REAL SD( MAXNSD ), RD( MAXNRD ), R( MAXNR ), TL( MAXNR ),
     &   RNEW( MAXNR ), TLNEW( MAXNR )
      COMPLEX P( MAXNR )
      CHARACTER TITLE*80, OPTION*8, TITLEY*80, PLTTYP*10,
     &   FILNAM*80

      DATA RINC, RAXL, TLMIN, TLMAX, TLINC, TLAXL /6* 0.0/

C     *** Read in plot parameters ***

      OPEN( FILE = 'PLPFIL', UNIT = PLPFIL, STATUS = 'OLD' )

      READ( PLPFIL, * ) OPTION
      READ( PLPFIL, * ) NCURVES
      READ( PLPFIL, * ) SRCD, RCVRD
      READ( PLPFIL, * )  RMIN,  RMAX,  RINC,  RAXL
      READ( PLPFIL, * ) TLMIN, TLMAX, TLINC, TLAXL
      READ( PLPFIL, * ) DR3DB

      WRITE( PRTFIL, 100 )
  100 FORMAT( '@   legend on' )

C     *** Loop over curve number ***

      SHDFIL = 0
      DO ICURVE = 1, NCURVES

C        *** Read in header information ***

         READ( PLPFIL, * ) FILNAM
C        --- Append the extension '.shd' to the filename
         JJ = INDEX( FILNAM, ' ' )
         FILNAM( JJ:JJ+3 ) = '.shd'

         CALL RDHEAD( SHDFIL, FILNAM, MAXNR, MAXNSD, MAXNRD,
     &      FREQ, R, NR, DELTAR, PLTTYP, XS, YS, THETA,
     &      SD, NSD, RD, NRD, TITLE )

C        --- Info for legend

         WRITE( PRTFIL, 500 ) ICURVE -1, FILNAM(1:JJ-1)
  500    FORMAT( '@   legend string ', I1, ' "', A, '"' ) 

C        ------ Extract range interval in kilometers

         IFIRST = 1
         ILAST  = NR

         DO IR = 1, NR
            R( IR ) = R( IR ) / 1000.0
            IF ( R( IR ) .LE. RMAX ) ILAST  = IR
            IF ( RMIN .GE. R( IR ) ) IFIRST = IR
         END DO

         NRNEW = ILAST - IFIRST + 1

C        --- Caculate tolerance for accepting source/rcvr depth

         IF ( FREQ .EQ. 0.0 ) THEN
            TOL = 0.001
         ELSE
            TOL = 1500.0 / ( 6.0 * FREQ )  ! one wavelength
         ENDIF

C        *** Calculate source depth index ***

         ISD = NEIGH( SRCD, SD, NSD )

         IF ( ABS( SD( ISD ) - SRCD ) .GT. TOL ) THEN
            WRITE( PRTFIL, * ) 'Sources exist at the following depths:'
            WRITE( PRTFIL, * ) ( SD( I ), I = 1, NSD )
            CALL ERROUT( PRTFIL, 'F', 'PLOTSLICE',
     &                'No source at specified depth' )
         ENDIF

C        *** Calculate receiver depth index ***

         IRD = NEIGH( RCVRD, RD, NRD )

         IF ( ABS( RD( IRD ) - RCVRD ) .GT. TOL ) THEN
            WRITE( PRTFIL, * ) 'Rcvrs exist at the following depths:'
            WRITE( PRTFIL, * ) ( RD( I ), I = 1, NRD )
            CALL ERROUT( PRTFIL, 'F', 'PLOTSLICE',
     &                'No receiver at specified depth' )
         ENDIF

C        *** Read in pressure field and convert to transmission loss ***

         IREC =  6 + NRD * ( ISD - 1 ) + IRD
         READ( SHDFIL, REC = IREC ) ( P( IR ), IR = 1, NR )

         TLMINT =  200.0
         TLMAXT = -200.0

         DO IR = 1, NRNEW
            IRORG = IFIRST + IR - 1
            RNEW( IR ) = R( IRORG )

            IF (OPTION(2:2) .EQ. 'L' .OR. OPTION(2:2) .EQ. 'F') THEN
C              ------ Linear scale
               IF ( OPTION(1:1) .EQ. 'S' ) THEN
                  TL(IR) = 1000.0 * SQRT( R(IRORG) ) * REAL( P(IRORG) )
               ELSE
                  TL(IR) = REAL( P( IRORG ) )
               ENDIF
            ELSE
C              ------ dB scale
               PSQ = P( IRORG ) * CONJG( P( IRORG ) )
               IF ( OPTION(1:1) .EQ. 'S' )
     &            PSQ = 1000.0 * R( IRORG ) * PSQ

               TL(IR) = 200.0
               IF ( PSQ .GT. 1.0E-20 ) TL( IR ) = -10.0 * ALOG10( PSQ )
            ENDIF

            WRITE( PRTFIL, * ) RNEW( IR ), TL( IR )

            TLMAXT = MAX( TL( IR ), TLMAXT )
            TLMINT = MIN( TL( IR ), TLMINT )
         END DO

         WRITE( PRTFIL, * ) '&'

C        Compute limits for transmission loss axis

C         WRITE( PRTFIL, * )
C         WRITE( PRTFIL, * ) 'TLMIN = ', TLMINT, 'TLMAX = ', TLMAXT

         IF ( ICURVE .EQ. 1 ) THEN
            ITLMIN = TLMINT / 10
            ITLMAX = TLMAXT / 10 + 1

            IF ( TLMIN .EQ. TLMAX ) THEN   ! autoscaling?
               TLMIN = 10.0 * ITLMIN
               TLMAX = 10.0 * ITLMAX
            ENDIF

            IF ( OPTION(2:2) .EQ. 'L' ) THEN
               TITLEY = 'Intensity$'
            ELSE
               TITLEY = 'Loss (dB)$'
               TITLEY = 'Power (dB)$'
            ENDIF

         ENDIF

C        Smooth

         NROWS  = 1
         NTHETA = 1
         CALL SMOOTH( TL, TLNEW, NROWS, NTHETA, NRNEW, DELTAR, DR3DB )

         SHDFIL = SHDFIL + 1
      END DO

C     *** Labels ***

      IF ( NCURVES .EQ. 1 ) THEN
      WRITE( PRTFIL, 1000 ) TITLE
 1000 FORMAT( '@ title "', A, '"' )
      ENDIF

      WRITE( PRTFIL, 2000 ) NINT( FREQ ), NINT( SRCD ), NINT( RCVRD )
 2000 FORMAT(  '@ subtitle "',
     &         'F  =', I4, ' Hz, ',
     &         'Sd =', I4, ' m,  ',
     &         'Rd =', I4, ' m"' )
 
      WRITE( PRTFIL, 3000 ) INT( RMIN  ), INT( RMAX  ),
     &                      INT( TLMIN ), INT( TLMAX )
 3000 FORMAT( '@ world xmin ', I5, /,
     &        '@ world xmax ', I5, /,
     &        '@ world ymin ', I5, /,
     &        '@ world ymax ', I5 )

      STOP
      END
