C      ALGORITHM 427, COLLECTED ALGORITHMS FROM ACM.
C      THIS WORK PUBLISHED IN COMMUNICATIONS OF THE ACM
C      VOL. 15, NO. 5, May, 1972, PP.358--360.
#! /bin/sh
# This is a shell archive, meaning:
# 1. Remove everything above the #! /bin/sh line.
# 2. Save the resulting text in a file.
# 3. Execute the file with /bin/sh (not csh) to create the files:
#	Fortran/
#	Fortran/Sp/
#	Fortran/Sp/Drivers/
#	Fortran/Sp/Drivers/Makefile
#	Fortran/Sp/Drivers/driver.f
#	Fortran/Sp/Drivers/res
#	Fortran/Sp/Src/
#	Fortran/Sp/Src/src.f
# This archive created: Thu Dec 15 13:28:02 2005
export PATH; PATH=/bin:$PATH
if test ! -d 'Fortran'
then
	mkdir 'Fortran'
fi
cd 'Fortran'
if test ! -d 'Sp'
then
	mkdir 'Sp'
fi
cd 'Sp'
if test ! -d 'Drivers'
then
	mkdir 'Drivers'
fi
cd 'Drivers'
if test -f 'Makefile'
then
	echo shar: will not over-write existing file "'Makefile'"
else
cat << "SHAR_EOF" > 'Makefile'
all: Res

src.o: src.f
	$(F77) $(F77OPTS) -c src.f

driver.o: driver.f
	$(F77) $(F77OPTS) -c driver.f

DRIVERS= driver
RESULTS= Res

Objs1= driver.o src.o
driver: $(Objs1)
	$(F77) $(F77OPTS) -o driver $(Objs1) $(SRCLIBS)
Res: driver 
	./driver >Res

diffres:Res res
	echo "Differences in results from driver"
	$(DIFF) Res res

clean: 
	rm -rf *.o $(DRIVERS) $(CLEANUP) $(RESULTS)
SHAR_EOF
fi # end of overwriting check
if test -f 'driver.f'
then
	echo shar: will not over-write existing file "'driver.f'"
else
cat << "SHAR_EOF" > 'driver.f'
      program main

c***********************************************************************
c
cc TOMS427_PRB tests FRCOS.
c
      implicit none

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'TOMS427_PRB'
      write ( *, '(a)' ) '  Test TOMS algorithm 427,'
      write ( *, '(a)' ) '  Fourier Cosine Integrals.'
      write ( *, '(a)' ) ' '

      call test01
      call test02

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'TOMS427_PRB'
      write ( *, '(a)' ) '  Normal end of execution.'

      stop
      end
      subroutine test01

c***********************************************************************
c
cc TEST01 handles F(X) = 1/(X*X+1).
c
      implicit none

      real et
      real exact
      external f1
      real f1
      real frcos
      real g1
      real hl
      integer i
      real t
      real value
      real w

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'TEST01:'
      write ( *, '(a)' ) '  f(x) = 1/(X*X+1),'
      write ( *, '(a)' ) '  Look at effect of upper limit T.'
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 
     &  '       W                T             VALUE           EXACT'
      write ( *, '(a)' ) ' '

      t = 5.0E+00
      w = 1.0E+00
      et = 0.0001E+00
      hl = 1.0E+00

      do i = 1, 3

        value = frcos ( f1, w, t, et, hl )

        exact = g1 ( w, t )

        write ( *, '(4(2x,g14.6))' ) w, t, value, exact

        t = t * 4.0

      end do

      return
      end
      function f1 ( x )

c***********************************************************************
c
cc F1 evaluates F(X) = 1/(X*X+1).
c
      implicit none

      real f1
      real x

      f1 = 1.0E+00 / ( x * x + 1.0E+00 )

      return
      end
      function g1 ( w, t )

c***********************************************************************
c
cc G1 evaluates G(W,T) = Integral ( 0 <= X <= T ) F1(X) * cos ( W * X ) dX.
c
      implicit none

      real g1
      real pi
      real t
      real w

      pi = 3.141592653589793

      g1 = pi * exp ( -w ) / 2.0E+00

      return
      end
      subroutine test02

c***********************************************************************
c
cc TEST02 handles F(X) = EXP(-X).
c
      implicit none

      real et
      real exact
      external f2
      real f2
      real frcos
      real g2
      real hl
      integer i
      real t
      real value
      real w

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'TEST02:'
      write ( *, '(a)' ) '  f(x) = exp(-x),'
      write ( *, '(a)' ) '  Look at effect of upper limit T.'
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 
     &  '       W                T             VALUE           EXACT'
      write ( *, '(a)' ) ' '

      t = 5.0E+00
      w = 1.0E+00
      et = 0.0001E+00
      hl = 1.0E+00

      do i = 1, 3

        value = frcos ( f2, w, t, et, hl )

        exact = g2 ( w, t )

        write ( *, '(4(2x,g14.6))' ) w, t, value, exact

        t = t * 4.0

      end do

      return
      end
      function f2 ( x )

c***********************************************************************
c
cc F2 evaluates F(X) = 1/(X*X+1).
c
      implicit none

      real f2
      real x

      f2 = exp ( -x )

      return
      end
      function g2 ( w, t )

c***********************************************************************
c
cc G2 evaluates G(W,T) = Integral ( 0 <= X <= T ) F2(X) * cos ( W * X ) dX.
c
      implicit none

      real g2
      real pi
      real t
      real w

      pi = 3.141592653589793

      g2 = 1.0E+00 / ( w * w + 1.0E+00 )

      return
      end

    

    
SHAR_EOF
fi # end of overwriting check
if test -f 'res'
then
	echo shar: will not over-write existing file "'res'"
else
cat << "SHAR_EOF" > 'res'
 
TOMS427_PRB
  Test TOMS algorithm 427,
  Fourier Cosine Integrals.
 
 
TEST01:
  f(x) = 1/(X*X+1),
  Look at effect of upper limit T.
 
       W                T             VALUE           EXACT
 
     1.00000         5.00000        0.571638        0.577864    
     1.00000         20.0000        0.577773        0.577864    
     1.00000         80.0000        0.577895        0.577864    
 
TEST02:
  f(x) = exp(-x),
  Look at effect of upper limit T.
 
       W                T             VALUE           EXACT
 
     1.00000         5.00000        0.499066        0.500000    
     1.00000         20.0000        0.500000        0.500000    
     1.00000         80.0000        0.500000        0.500000    
 
TOMS427_PRB
  Normal end of execution.
SHAR_EOF
fi # end of overwriting check
cd ..
if test ! -d 'Src'
then
	mkdir 'Src'
fi
cd 'Src'
if test -f 'src.f'
then
	echo shar: will not over-write existing file "'src.f'"
else
cat << "SHAR_EOF" > 'src.f'
      REAL FUNCTION FRCOS ( FC, W, T, ET, HL )
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  ARRAYS ER, W1C, W2C, W3C CONTAIN PRECOMPUTED CONSTANTS
C  NEEDED TO COMPUTE APPROXIMATE VALUES AND ERROR
C  ESTIMATES FOR FILON2 (SEE COMMENTS).
C     .. Scalar Arguments ..
      REAL ET,HL,T,W
C     ..
C     .. Function Arguments ..
      REAL FC
      EXTERNAL FC
C     ..
C     .. Local Scalars ..
      REAL A,ALN2,B,CO1,CONST,EPS,ERC,ERR,ERT,F1,F2,F3,FQ,FQ3,H,HI,PI2,
     +     PI256,ROC,SI1,TA,VAL,VNEW,VNEW1,VNEW2,W1,W2,W3,WA,WA1,WB,WHI,
     +     WT,XQ,XQ3
      INTEGER N,N2,NH,NP
C     ..
C     .. Local Arrays ..
      REAL AS(30),ER(9),FS(61),PVAL(30),W1C(9),W2C(9),W3C(9)
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC ABS,ALOG,COS,IFIX,SIN
      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 /
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
SHAR_EOF
fi # end of overwriting check
cd ..
cd ..
cd ..
#       End of shell archive
exit 0