      PROGRAM PLOTTS

C     Plot the source time series

      INTEGER PLPFIL
      PARAMETER ( PLPFIL = 1, MAXNT = 500000, MAXNSD = 51,
     &            MAXNST = 500000, PI = 3.141592 )

      REAL TV( MAXNT ), SVR( MAXNST ), SVI( MAXNST ), SD( MAXNSD )
      COMPLEX SV( MAXNST )
      CHARACTER PULSE*4, PULSTITL*80

C     *** Read in data ***

      NT = 3126

      OPEN( FILE = 'PLPFIL', UNIT = PLPFIL, STATUS = 'OLD' )
      READ( PLPFIL, * ) FREQ, FMIN, FMAX, PULSE

      IF ( PULSE(4:4) .EQ. 'H' .OR. PULSE(4:4) .EQ. 'N' )
     &   FMIN = 0.0   ! no low-cut filtering

      IF ( PULSE(4:4) .EQ. 'L' .OR. PULSE(4:4) .EQ. 'N' )
     &   FMAX = 1.0E9 ! no high-cut filtering

      READ( PLPFIL, * ) TMIN, TMAX

      OMEGA  = 2.0 * PI * FREQ
      DELTAT = (TMAX - TMIN) / ( NT - 1 )

C     *** Construct the time series ***

      DO IT = 1, NT
         TV( IT ) = TMIN + ( IT - 1 ) * DELTAT
      END DO

      CALL SOURCE( TV, SV, SD, NSD, NT, OMEGA, FMIN, FMAX, PULSE,
     &      PULSTITL )

      CALL SCALETS( TV, SV, SD, NSD, NT, PULSTITL, SVR, SVI,
     &   TMIN, TMAX, TINC, TAXL )

      STOP
      END

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

      SUBROUTINE SCALETS( TV, SV, SD, NSD, NT, PULSTITL, SVR, SVI,
     &   TMIN, TMAX, TINC, TAXL )

C     Scale and plot the time series data

      REAL TV( * ), SVR( NT, * ), SVI( NT, * ), SD( * )
      COMPLEX SV( NSD, * )
      CHARACTER TITLEY*80, PULSTITL*80

C     *** Write data ***

      IF ( NSD .EQ. 1 ) THEN
         TITLEY = 'Amplitude$'   ! Single curve
      ELSE
         TITLEY = 'Depth (m)$'   ! Multiple curves
      ENDIF

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

C     WRITE( *, 2000 ) NINT( FREQ ), NINT( SRCD ), NINT( RCVRD )
 2000 FORMAT(  '@ subtitle "',
     &         'F  =', I4, ' Hz, ',
     &         'Sd =', I4, ' m,  ',
     &         'Rd =', I4, ' m"' )

C     *** compute maximum ***

      SMAX = 1.0E-20

      DO IT = 1, NT
         DO IS = 1, NSD
            SVR( IT, IS ) =  REAL( SV( IS, IT ) )
            SVI( IT, IS ) = AIMAG( SV( IS, IT ) )

C           --- we do this in steps to avoid SGI compiler error:
            SRABS = ABS( SVR( IT, IS ) )
            SIABS = ABS( SVI( IT, IS ) )
            SMAX  = MAX( SRABS, SIABS, SMAX )
         END DO
      END DO

C     --- avoid scaling problems for vanishing time series
      IF ( SMAX .EQ. 0.0 ) SMAX = 1.0

      DO IT = 1, NT
         DO IS = 1, NSD
            SVR( IT, IS ) = SVR( IT, IS )
            SVI( IT, IS ) = SVI( IT, IS )
         END DO
      END DO

C     *** Plot ***

      DO IS = 1, NSD

         WRITE( *, 3000 ) NSD - IS, NSD - IS
 3000    FORMAT( '@WITH G', I1.1, /,
     &           '@G', I1.1, ' ON', /,
     &           '@TYPE xy' )

C        --- real part
         DO IT = 1, NT
            WRITE( *, * ) TV( IT ), SVR( IT, IS )
         END DO
         WRITE( *, * ) '&'

C        --- imaginary part
         DO IT = 1, NT
            WRITE( *, * ) TV( IT ), SVI( IT, IS )
         END DO
         WRITE( *, * ) '&'
      END DO

      RETURN
      END
