FUNCTION FRCOS ( FC, W, T, ET, HL )

c*********************************************************************72
c
C  THIS ROUTINE COMPUTES THE FOURIER COSINE INTEGRAL FROM
C  ZERO TO INFINITY OF FC(T) * COS(W*T) BY AN ADAPTIVE
C  QUADRATURE METHOD USING A COMBINATION OF FILON AND
C  SIMPSON RULES.
C    PARAMETERS
C  FC -MUST BE DECLARED EXTERNAL IN CALLING PROGRAM.
C  W -VALUE MUST BE NON-NEGATIVE.
C  T -UPPER LIMIT FOR QUADRATURE - SHOULD NORMALLY BE CHOSEN
C  SUCH THAT REST OF INTEGRAL IS NEGLIGIBLE.  THE ACTUAL
C  LIMIT USED BY THE PROGRAM MAY BE SOMEWHAT LARGER THAN
C  THE GIVEN T (SEE INTRODUCTORY COMMENTS).
C  ET -REQUESTED ACCURACY (ABSOLUTE).
C  HL -LIMIT ON STEP SIZE-CONVERGENCE IN ANY SUBINTERVAL IS
C  NOT RECOGNIZED UNLESS SUBINTERVAL IS SMALLER THAN HL.
c
      DIMENSION W1C(9), W2C(9), W3C(9), ER(9)
c
C  ARRAYS ER, W1C, W2C, W3C CONTAIN PRECOMPUTED CONSTANTS
C  NEEDED TO COMPUTE APPROXIMATE VALUES AND ERROR
C  ESTIMATES FOR FILON2 (SEE COMMENTS).
c
      DATA ER(1),ER(2),ER(3),ER(4),ER(5),ER(6),ER(7),ER(8),
     &  ER(9) / 0., .05061,.05969,.06181,.06233,.06246,
     &  .06249,.06249,.0625 /
      DATA W1C(1),W1C(2),W1C(3),W1C(4),W1C(5),W1C(6),W1C(7),
     &  W1C(8),W1C(9) /
     &  4.0528473456934E-01, 3.6761020369133E-01,
     &  3.4316760755741E-01, 3.3587533234962E-01,
     &  3.3397411782348E-01, 3.3349386085934E-01,
     &  3.3337348594489E-01, 3.3334337278212E-01,
     &  3.3333584327653E-01 /
      DATA W2C(1),W2C(2),W2C(3),W2C(4),W2C(5),W2C(6),W2C(7),
     &  W2C(8),W2C(9) /
     &  1.2059522143639E-01, 1.9710810149097E-02,
     &  2.6328277852505E-03, 3.3459141708323E-04,
     &  4.1997086077777E-05, 5.2550600306570E-06,
     &  6.5705211443498E-07, 8.2136815416350E-08,
     &  1.0267267595664E-08 /
      DATA W3C(1),W3C(2),W3C(3),W3C(4),W3C(5),W3C(6),W3C(7),
     &  W3C(8),W3C(9) /
     &  1.0320491018624, 1.2528780015490, 1.3128845799752,
     &  1.3281999871557, 1.3320486708792, 1.3330120847949,
     &  1.3332530160151, 1.3333132536798, 1.3333283133997 /
      DIMENSION FS(61),PVAL(30),AS(30)
C  ARRAYS FS,PVAL,AS ARE STORAGE FOR SAVED VALUES OF F
C  AND BOOK-KEEPING.
      DATA PI2,PI256 / 6.2831853071796, .0122718463 /
C  PI2 = 2*PI, PI256 = PI / 256.
      DATA ALN2, ERC, ROC / .69314718, 1.E30, 1.E-5 /
C  ALN2 = NATURAL LOG OF 2, ERC = ERROR VALUE RETURNED
C  BY FRCOS, ROC = CONSTANT USED TO ELIMINATE ROUNDOFF
C  PROBLEMS IN COMPUTING INTERVAL LIMITS.
      EPS = ET
      VAL = 0.
      N = 1
      AS(1) = 0.
      FS(1) = FC ( 0. )
      PVAL(1) = ERC
C  TEST IF UPPER LIMIT ADJUSTMENT IS NECESSARY.
      WT = W * T
      IF ( WT - PI256 + ROC ) 100, 100, 101
C  NOTE - CONSTANT ROC = 1.E-5 USED THROUGHOUT PROGRAM TO
C  ELIMINATE EFFECT OF FLOATING POINT ROUNDOFF ERROR.
C
C  SET UP FIRST INTERVAL FOR SIMPSON RULE.
100   FS(2) = FC ( .5 * T ) * COS ( .5 * WT )
      FS(3) = FC (      T ) * COS (      WT )
      B = T
      GO TO 105
C  ADJUST UPPER LIMIT.
101   NP = IFIX ( ALOG ( WT / PI256 ) / ALN2 ) + 1
      TA = 2**NP * PI256 / W
      B = TA
C  SET UP FIRST INTERVAL FOR FILON RULE.
      FS(2) = FC ( .5 * TA )
      FS(3) = FC (      TA )
C  TAKE LAST INTERVAL FROM LIST.
105   A = AS(N)
      HI = B - A
      WHI = W * HI
      N2 = 2 * N
      F1 = FS(N2-1)
      F2 = FS(N2)
      F3 = FS(N2+1)
      XQ  = B - .75 * HI
      XQ3 = B - .25 * HI
C  TEST TO DETERMINE WHICH QUADRATURE RULE IS APPLICABLE.
      IF ( WHI - PI256 - ROC ) 110, 110, 111
110   IF ( WHI - PI256 + ROC ) 200, 200, 201
111   IF ( WHI - PI2   - ROC ) 220, 220, 230
C  ESTIMATE BY SIMPSON RULE.
200   FQ = FC ( XQ ) * COS ( W * XQ )
      FQ3 = FC ( XQ3 ) * COS ( W * XQ3 )
      VNEW1 = HI * ( F1 + 4. * FQ  + F2 ) / 12.
      VNEW2 = HI * ( F2 + 4. * FQ3 + F3 ) / 12.
      VNEW = VNEW1 + VNEW2
      ERR = ( PVAL(N) - VNEW ) / 15.
      GO TO 300
C  SWITCH FROM FILON TO SIMPSON RULE.
201   F1 = F1 * COS ( W * A )
      F2 = F2 * COS ( W * ( B - .5 * HI ) )
      F3 = F3 * COS ( W * B )
      PVAL(N) = HI * ( F1 + 4. * F2 + F3 ) / 6.
      GO TO 200
C  ESTIMATE BY FILON2.
220   H = .25 * HI
      FQ  = FC ( XQ )
      FQ3 = FC ( XQ3 )
      NH = IFIX ( ALOG ( PI2 / WHI ) / ALN2 + ROC ) + 1
      W1 =  W1C(NH)
      W2 = -W2C(NH)
      W3 =  W3C(NH)
      WA  = W * A
      WA1 = W * ( B - .5 * HI )
      WB  = W * B
      CO1 = COS ( WA1 )
      SI1 = SIN ( WA1 )
      VNEW1 = H*((W1*COS(WA)+W2*SIN(WA))*F1 + W3*COS(W*XQ)*
     &  FQ+(W1*CO1-W2*SI1)*F2)
      VNEW2 = H*((W1*CO1 +W2*SI1)*F2 + W3*COS(W*XQ3)*FQ3
     &  +(W1*COS(WB) - W2*SIN(WB))*F3)
      VNEW = VNEW1 + VNEW2
      ERT = ER(NH)
      ERR = ERT * ( PVAL(N) - VNEW ) / ( 1. - ERT )
C  SKIP CONVERGENCE TEST IF INTERVAL = ONE PERIOD.
      IF ( WHI - PI2 + ROC ) 300, 300, 400
C  ESTIMATE BY FILON1.
230   FQ = FC ( XQ )
      FQ3 = FC ( XQ3 )
      W2 = W * W
      CONST = 8. / ( W2 * HI )
      VNEW1 = CONST * ( F1 - 2. * FQ  + F2 )
      VNEW2 = CONST * ( F2 - 2. * FQ3 + F3 )
      VNEW = VNEW1 + VNEW2
      W2 = 6. / W2
      W3 = HI * HI
      ERT = ( W3 / 32. - W2 ) / ( W3 / 8. - W2 )
      ERR = ERT * ( PVAL(N) - VNEW ) / ( 1. - ERT )
C  CONVERGENCE TEST.
C  SKIP CONVERGENCE TEST IF HI.GT.HL
300   IF ( HI - HL ) 301, 301, 400
301   IF ( ABS ( ERR ) - EPS * HI / B ) 500, 500, 400
C  CONVERGENCE NOT OBTAINED -SPLIT INTERVAL AND ADD TO LIST.
C  TEST FOR POSSIBLE LIST OVERFLOW.
400   IF ( N - 30 ) 401, 600, 600
401   FS(N2+3) = F3
      FS(N2+2) = FQ3
      FS(N2+1) = F2
      FS(N2)   = FQ
      AS(N+1) = A + .5 * HI
      PVAL(N)   = VNEW1
      PVAL(N+1) = VNEW2
      N = N + 1
      GO TO 105
C  CONVERGENCE OBTAINED -ADD EXTRAPOLATED PARTIAL SUM TO
C  TOTAL--ADJUST ERROR AND INTERVAL.
500   VAL = VAL + VNEW - ERR
      EPS = EPS - ABS ( ERR )
      N = N - 1
      B = A
      IF ( N ) 700, 700, 105
C  CONVERGENCE FAILURE -ROUTINE RETURNS ERC=1.E+30
C  OPTIONAL ERROR MESSAGE MAY BE INSERTED HERE.
600   FRCOS = ERC
      RETURN
C  COMPUTATIONS COMPLETED SUCCESSFULLY.
700   FRCOS = VAL

      RETURN
      END