      SUBROUTINE SOURCE( T, S, SD, NSD, NTOUT, OMEGA, FMIN, FMAX,
     &       PULSE, PULSTITL )

C     Evaluate the source time series

C     STEMP is a temporary workspace
C     Assumes that the filter parameters will not change on
C     subsequent calls. (The time series is only filtered
C     when it's first read.)

      PARAMETER ( MAXNT = 1000000, MAXNSD = 11,
     &            MAXNST = 1000000 )

      LOGICAL INIFLG
      REAL T( * ), TF( MAXNT ), SD( MAXNSD ), STEMPR( MAXNT )
      COMPLEX S( NSD, * ), SF( MAXNST ), STEMP( MAXNT )
      CHARACTER PULSE*( * ), PULSTITL*( * )
      EQUIVALENCE ( STEMPR, STEMP )

      SAVE INIFLG, TF, SF, NT
      DATA INIFLG /.TRUE./

C     *** If first call, tabulate the time series ***

      IF ( INIFLG ) THEN
         IF ( PULSE(1:1) .EQ. 'F' .OR. PULSE(1:1) .EQ. 'B' ) THEN
C           ------ From a file
            CALL STSHDR( PULSTITL, SD, NSD, MAXNSD )
            CALL SFILE( PULSE, PULSTITL, NSD, TF,
     &                  STEMPR, SF, NT, MAXNT, MAXNST )
            DELTAT = TF( 2 ) - TF( 1 )
         ELSE
C           ------ Use one of the canned wavelets
            NSD = 1
            DELTAT = 0.2 / OMEGA
            TSTART = -200.0 * DELTAT
            NT = 1024
            DO IT = 1, NT
               TF( IT ) = TSTART + ( IT - 1 ) * DELTAT
               CALL CANS( TF( IT ), OMEGA, PULSE,
     &                    SF, NSD, IT, PULSTITL )
            END DO
         ENDIF

C        *** Filter the time series ***
         CALL FILTER( PULSE, DELTAT, STEMP, SF, NSD, NT, FMIN, FMAX )

         INIFLG = .FALSE.
      ENDIF

      CALL EVALU( T, PULSE, TF, SF, NSD, NT, S, NTOUT )

      RETURN
      END
C**********************************************************************C
      SUBROUTINE STSHDR( PULSTITL, SD, NSD, MAXNSD )

C     Read header from time series file


      INTEGER PRTFIL, STSFIL
      PARAMETER ( PRTFIL = 6, STSFIL = 10 )

      REAL SD( * )
      CHARACTER PULSTITL*( * )

C     *** Read in time series data ***

      OPEN( FILE = 'STSFIL', UNIT = STSFIL, STATUS = 'OLD',
     &       FORM = 'FORMATTED' )

      READ( STSFIL, * ) PULSTITL

C     --- SD is used by PLOTTS
      READ( STSFIL, * ) NSD, ( SD( IS ), IS = 1, NSD )

      IF ( NSD .GT. MAXNSD )
     &   CALL ERROUT( PRTFIL, 'F', 'SOURCE:SFILE',
     &             'Max. number of sources exceeded' )

      RETURN
      END
C**********************************************************************C
      SUBROUTINE SFILE( PULSE, PULSTITL, NSD, TF,
     &    STEMP, SF, NT, MAXNT, MAXNST )

C     Time series from file

C     STEMP is used to hold a temporary real vector

      INTEGER PRTFIL, STSFIL
      PARAMETER ( PRTFIL = 6, STSFIL = 10 )

      REAL    STEMP( * ), TF( * )
      COMPLEX SF( NSD, * )
      CHARACTER PULSE*( * ), PULSTITL*( * )

C     *** Read in time series data ***

      NT = 0
      DO IT = 1, MAXNT

         READ( STSFIL, *, END = 2000 )
     &      TF(IT), ( STEMP( IS ), IS = 1, NSD )

         IF ( IT .GT. 1 ) THEN
            IF ( TF( IT ) .LT. TF( IT - 1 ) )
     &         CALL ERROUT( PRTFIL, 'F', 'SOURCE:SFILE',
     &             'Time series not ordered in time' )
         ENDIF

         DO IS = 1, NSD
            SF( IS, IT ) = STEMP( IS )
         END DO

         NT = NT + 1

         IF ( NT * NSD .GT. MAXNST ) 
     &        CALL ERROUT( PRTFIL, 'F', 'SOURCE:SFILE',
     &                     'Too many time series points' )
      END DO

C     --- normal exit is by EOF in file, not by falling through loop
      CALL ERROUT( PRTFIL, 'F', 'SOURCE:SFILE',
     &             'Too many time series points' )

C     *** Time reversal ***
C     note if NT is odd, middle point doesn't move

 2000 IF ( PULSE(1:1) .EQ. 'B' ) THEN
         TMAX = TF( NT )
         DO IT = 1, NT / 2

            DO IS = 1, NSD
               TEMP                  = SF( IS, IT )
               SF( IS, IT )          = SF( IS, NT - IT + 1 )
               SF( IS, NT - IT + 1 ) = TEMP
            END DO

            TEMP                  = TMAX - TF( IT )
            TF( IT )              = TMAX - TF( NT - IT + 1 )
            TF( NT - IT + 1 )     = TEMP

         END DO
      ENDIF

      CLOSE( STSFIL )

      RETURN
      END
C**********************************************************************C

      SUBROUTINE EVALU( T, PULSE, TF, SF, NSD, NT, S, NTOUT )

C     *** Evaluate the source function ***

      REAL T( * ), TF( * )
      COMPLEX S( NSD, * ), SF( NSD, * )
      CHARACTER PULSE*( * )

      SAVE IT
      DATA IT /1/

      DO ITOUT = 1, NTOUT

C        Necessary to advance time window?

C        ------ Case of T reset to origin for new march
         IF ( T( ITOUT ) .LT. TF( IT ) ) IT = 1
 5000    IF ( T( ITOUT ) .GT. TF( IT + 1 ) .AND. IT .LT. NT - 1 ) THEN
            IT = IT + 1
            GO TO 5000
         ENDIF

C        Linear interpolation in time

         DO IS = 1, NSD
            S( IS, ITOUT ) = 0.0   ! in case T outside time window
         END DO

         IF ( T( ITOUT ) .GE. TF( 1  ) .AND.
     &        T( ITOUT ) .LE. TF( NT ) ) THEN
            W = ( T( ITOUT ) - TF( IT ) ) / ( TF( IT + 1 ) - TF( IT ) )
            DO IS = 1, NSD
               S( IS, ITOUT ) = SF( IS, IT )
     &            + W * ( SF( IS, IT+1 ) - SF(IS, IT) )
               IF ( PULSE(3:3) .EQ. '-' )
     &                 S( IS, ITOUT ) = -S( IS, ITOUT )
            END DO
         ENDIF

      END DO   ! next ITOUT

      RETURN
      END
C**********************************************************************C

      SUBROUTINE FILTER( PULSE, DELTAT, STEMP, SF, NSD, NT, FMIN, FMAX )

C     *** Evaluate the source function ***

      COMPLEX STEMP( * ), SF( NSD, * )
      CHARACTER PULSE*( * )

C     NOTE: program will abort if NT is not a power of 2

      DO IS = 1, NSD
         DO I = 1, NT   ! Copy time series into STEMP
            STEMP( I ) = SF( IS, I ) 
         END DO

C        --- Band-pass filter
         IF ( PULSE(4:4) .NE. 'N' )
     &      CALL BNDPAS( STEMP, NT, DELTAT, FMIN, FMAX )

C        --- Form the pre-envelope
         IF ( PULSE(2:2) .EQ. 'H' ) CALL PREENV(  STEMP, NT )

C        --- Hilbert transform
         IF ( PULSE(2:2) .EQ. 'Q' ) CALL HILBERT( STEMP, NT )

         DO I = 1, NT   ! Copy time series back
            SF( IS, I ) = STEMP( I )
         END DO
      END DO

      RETURN
      END
