      SUBROUTINE CANS( T, OMEGA, PULSE, SF, NSD, IT, PULSTITL )

C     Computes the source time series

C     T      is the time
C     OMEGA  is some frequency characterizing the pulse
C     PULSE  is the letter indicating which pulse
C     S      is the time series

      PARAMETER ( PI = 3.141592 )

      REAL S
      COMPLEX SF( NSD, * )
      CHARACTER PULSE*( * ), PULSTITL*( * )

      S = 0.0
      F = OMEGA / ( 2.0 * PI )

      IF ( T .GT. 0.0 ) THEN

C     EVALUATE SELECTED PULSE

      IF ( PULSE(1:1) .EQ. 'P' ) THEN
C        ------ Pseudo gaussian
C        (peak at 0, support [0, 3F]
         IF ( T .LE. 1.0 / F )
     &   S = 0.75 - COS( OMEGA * T ) + 0.25*COS( 2.0 * OMEGA * T )
         PULSTITL = 'Pseudo gaussian'

      ELSE IF ( PULSE(1:1) .EQ. 'R' ) THEN
C        ------ Ricker wavelet
C        (peak at F, support [0, 2F]
         U = OMEGA * T - 5.0
         S = 0.5 * ( 0.25*U*U - 0.5 ) * SQRT( PI ) * EXP( -0.25*U*U )
         PULSTITL = 'Ricker wavelet'

      ELSE IF ( PULSE(1:1) .EQ. 'A' ) THEN
C        ------ Approximate Ricker wavelet
C        (peak at F, support [0, 2.5F]
         TC = 1.55 / F
         IF ( T .LE. TC )
     &   S = 0.48829 *    COS( 2.0 * PI * T / TC )
     &         -0.14128 * 4* COS( 4.0 * PI * T / TC )
     &         +0.01168 * 9* COS( 6.0 * PI * T / TC )
         PULSTITL = 'Approximate Ricker wavelet'

      ELSE IF ( PULSE(1:1) .EQ. 'S' ) THEN
C        ------ Single sine
C        (peak at F, support [0, infinity], nulls at nF
         IF ( T .LE. 1.0 / F ) S = SIN( OMEGA * T )
         PULSTITL = 'Single sine'

      ELSE IF ( PULSE(1:1) .EQ. 'H' ) THEN
C        ------ Hanning weighted four sine
C        (peak at F, support [0, infinity], first null at about 1.5F
         IF ( T .LE. 4.0 / F )
     &   S = 0.5 * SIN( OMEGA * T ) * ( 1 - COS( OMEGA * T / 4.0 ) )
         PULSTITL = 'Hanning weighted four sine'

      ELSE IF ( PULSE(1:1) .EQ. 'N' ) THEN
C        ------ N-wave
C        (peak at 1 F, support [0, 4F], [0,3F] OK
         IF ( T .LE. 1.0 / F )
     &   S = SIN( OMEGA * T ) - 0.5 * SIN( 2.0 * OMEGA * T )
         PULSTITL = 'N-wave'

      ELSE IF ( PULSE(1:1) .EQ. 'M' ) THEN
C        ------ Miracle wave (has a Hilbert transform?)
C        (peak at 0, support [0, infinity]
         A = 1.0 / ( 6.0 * F )
         T0 = 6.0 * A
         TS = ( T - T0 ) / A
         S = 1.0 / ( 1.0 + TS * TS )
C        HS is the Hilbert transform of the time series
C        HS = TS / ( 1.0 + TS * TS )
         PULSTITL = 'Miracle wave'

      ELSE IF ( PULSE(1:1) .EQ. 'G' ) THEN
C        ------ Gaussian
C        (peak at 0, support [0, infinity]
C        T0 is the peak time
C        A is the 3dB down time
C        NOTE S(0) = EXP( -NSIG ** 2 )

         NSIG = 3
         A = 1.0 / F / ( 2.0 * NSIG )
         T0 = NSIG * A
         S = EXP( -( ( T - T0 ) / A ) ** 2 )
         PULSTITL = 'Gaussian'
      ELSE IF ( PULSE(1:1) .EQ. 'T' ) THEN
C        ------ Tone
C        (peak at F, support [0, infinity]

         IF ( T .LE. 0.4 ) S = SIN( OMEGA * T )
         PULSTITL = 'Tone'
      ELSE IF ( PULSE(1:1) .EQ. 'C' ) THEN
C        ------ Sinc
C        ( uniform spectrum from [0, F]

         S = SIN( OMEGA * T ) / ( OMEGA * T )
         PULSTITL = 'Sinc'
      END IF
      END IF

C     *** Copy into source vector ***

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

      RETURN
      END
