      SUBROUTINE HTS( NT, WKMIN, WKDEL, STABIL,
     &                RMIN, NR, RDEL, X, Y, Z, INIFLG, OPTION )

C     Hankel Transform by series approximation - Currently first term only.
C     Formulas from Abramowitz and Stegun, Eqs. 9.2.5-9.2.10

C     X contains input k space spectra, NT points.
C     Z contains output transmission loss, NR points.
C     NT = number of points in transform.
C     WKMIN = minimum wavenumber of input series.
C     WKDEL = spacing of input series.
C     STABIL = amount integration has been moved off of the real axis.
C     NR = number of output transmission loss points.(WITH NR .LE. 2*NT)
C     RDEL = spacing of transmission loss points.
C     RMIN = STARTING POINT FOR TL

C     Adapted from Steve Wales' original program

C     This version uses the FPS library FFT.
C     To revert to original:
C        allocate WSPACE
C        take the commented set of FFT calls

      COMPLEX AI, AIPID4
      PARAMETER ( AI = (0.0, 1.0), PI = 3.141592653589793,
     &            TWOPI = 2. * PI )

      LOGICAL  INIFLG
C     REAL     WSPACE( 4 * MAXNT + 5 )
      COMPLEX  X( * ), Y( * ), Z( * ), CK, AISTAB, AIRMIN, Z2
      CHARACTER OPTION*80

C     ------ Initialization

      AIPID4 = AI * PI / 4.0
      AISTAB = AI * STABIL
      AIRMIN = AI * RMIN
      SQRT2P = 1.0 / SQRT ( TWOPI )
      FACTOR = WKDEL * SQRT2P

C      IF ( INIFLG ) CALL CFFTI ( NT, WSPACE )

C     *** First term, lower-half of cosine transform ***

      IF ( WKMIN .GT. 0.5 * WKDEL )  THEN
         X( 1 ) = FACTOR * SQRT( WKMIN + AISTAB ) * X( 1 )
      ELSE
         X( 1 ) = 0.5 * ( AISTAB + WKDEL ) * SQRT2P
     &                   * SQRT(         AISTAB ) * X( 1 )
      END IF

      Y( 1 ) = X( 1 ) * EXP( AIPID4 )

      DO IS = 2, NT
         WKDIF = ( IS - 1 ) * WKDEL
         X( IS ) = FACTOR * SQRT( WKMIN + WKDIF + AISTAB ) * X( IS )
         Y( IS ) = X( IS ) * EXP( -AIRMIN * WKDIF + AIPID4 )
      END DO

C     ------ (P)ositive or (B)oth positive and negative spectrum

      IF ( OPTION(3:3) .EQ. 'P' .OR. OPTION(3:3) .EQ. 'B' ) THEN
         CALL CFFT( Y, NT, 1 )                 ! EXP (-IKX) transform
C        CALL CFFTF( NT, Y, WSPACE )

         CK = CMPLX ( +STABIL, -WKMIN )
         DO IR = 1, NR
            R = RMIN + ( IR - 1 ) * RDEL
            Z( IR ) = Y( IR ) * EXP ( CK * R )
         END DO
      ENDIF

C     *** Second term, second-half of cosine transform ***

C     ------ (N)egative or (B)oth positive and negative spectrum

      IF ( OPTION(3:3) .EQ. 'N' .OR. OPTION(3:3) .EQ. 'B' ) THEN
         Y( 1 ) = X( 1 ) * EXP( -AIPID4 )
         DO IS = 1, NT
            WKDIF   = ( IS - 1 ) * WKDEL
            Y( IS ) = X( IS ) * EXP( +AIRMIN * WKDIF - AIPID4 )
         END DO

         CALL CFFT( Y, NT, -1 )                  ! EXP (+IKX) transform
C        CALL CFFTB( NT, Y, WSPACE )
      ENDIF

      DO IR = 1, NR
         R  = RMIN + ( IR - 1 ) * RDEL
         Z2 = Y( IR ) * EXP ( -CK * R )

         IF ( OPTION(3:3) .EQ. 'B' ) THEN
            Z( IR ) = Z( IR ) + Z2
         ELSE IF ( OPTION(3:3) .EQ. 'N' ) THEN         
            Z( IR ) = Z2
         ENDIF

         IF ( OPTION(1:1) .EQ. 'R' ) Z( IR ) = Z( IR ) / SQRT( R )
      END DO

      RETURN
      END
C**********************************************************************C
      SUBROUTINE FTS ( NT, WKMIN, WKDEL, STABIL,
     &                RMIN, NR, RDEL, X, Y, Z, INIFLG )

C     Fourier Transform

C     X contains input k space spectra, NT points.
C     Z contains output transmission loss, NR points.
C     NT = number of points in transform.
C     WKMIN = minimum wavenumber of input series.
C     WKDEL = spacing of input series.
C     STABIL = amount integration has been moved off of the real axis.
C     NR = number of output transmission loss points.(WITH NR .LE. 2*NT)
C     RDEL = spacing of transmission loss points.
C     RMIN = STARTING POINT FOR TL

C     Adapted from Steve Wales' original program

      COMPLEX AI
      PARAMETER ( AI = (0.0, 1.0), PI = 3.141592653589793,
     &   TWOPI = 2.*PI )

      LOGICAL INIFLG
      COMPLEX  X( * ), Y( * ), Z( * ), CK, AISTAB, AIRMIN, Z2

C     ------ Initialization
      AISTAB = AI * STABIL
      AIRMIN = AI * RMIN
      SQRT2P = 1. / SQRT ( TWOPI )
      FACTOR = WKDEL * SQRT2P

C     *** First term, lower-half of cosine transform ***

      IF ( WKMIN .GT. 0.5 * WKDEL )  THEN
         X( 1 ) = FACTOR * X( 1 )
      ELSE
         X( 1 ) = ( AISTAB + WKDEL ) / 2.0 * SQRT2P * X( 1 )
      END IF

      Y( 1 ) = X( 1 )

      DO IS = 2, NT
         WKDIF = ( IS - 1 ) * WKDEL
         X( IS ) = FACTOR * X( IS )
         Y( IS ) = X( IS ) * EXP( -AIRMIN * WKDIF )
      END DO

      CALL CFFT (Y, NT, 1)                      ! EXP (-IKX) transform
      CALL CFFTSC ( Y, NT )                     ! Scale
      CK = CMPLX( +STABIL, -WKMIN )

      DO IR = 1, NR
         R = RMIN + ( IR - 1 ) * RDEL
         Z( IR ) = Y( IR ) * EXP ( CK * R )
      END DO

C     *** Second term, second-half of cosine transform ***

      Y( 1 ) = X( 1 )

      DO IS = 1, NT
         WKDIF = ( IS - 1 ) * WKDEL
         Y( IS ) = X( IS ) * EXP( +AIRMIN * WKDIF )
      END DO

      CALL CFFT( Y, NT, -1 )                     ! EXP (+IKX) transform

      DO IR = 1, NR
         R  = RMIN + (IR-1) * RDEL
         Z2 = Y( IR ) * EXP ( -CK * R )
         Z( IR ) = Z( IR ) + Z2
      END DO

      RETURN
      END
