C      ALGORITHM 353, COLLECTED ALGORITHMS FROM ACM.
C      THIS WORK PUBLISHED IN COMMUNICATIONS OF THE ACM
C      VOL. 12, NO. 8, August, 1969, PP.457--458.
#! /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: Wed Jan 18 20:30:14 2006
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 TOMS353_PRB tests FSPL2
c
c  Modified:
c
c    17 January 2006
c
c  Author:
c
c    John Burkardt
c
      implicit none

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'TOMS353_PRB:'
      write ( *, '(a)' ) '  Test ACM algorithm 353, which uses'
      write ( *, '(a)' ) '  Filon quadrature to approximate integrals'
      write ( *, '(a)' ) '  of functions of the form F(X)*COS(M*PI*X)'
      write ( *, '(a)' ) '  or F(X)*SIN(M*PI*X).'
 
      call test01

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

      write ( *, '(a)' ) ' '

      stop
      end
      subroutine test01

c***********************************************************************
c
cc TEST01 tests FSER1 with integrands of the form F(X)*COS(M*PI*X).
c
      implicit none

      real a
      real b
      real c
      real eps
      real exact
      real f1, f2, f3
      external f1
      external f2
      external f3
      integer lc
      integer ls
      integer m
      integer max
      real pi
      real s

      pi = 3.141592653589793E+00
      a = 0.0E+00
      b = 1.0E+00

      eps = 0.00001E+00
      max = 8

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'TEST01'
      write ( *, '(a)' ) '  Use FSER1 to estimate the integrals of'
      write ( *, '(a)' ) '  the form F(X) * COS ( M * PI * X )'
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) '  Use integrand factors:'
      write ( *, '(a)' ) '  F(X) = 1, X, X*X.'
      write ( *, '(a)' ) ' '
      write ( *, '(a,a)' ) 
     &  '               M  Approximate         Exact'
      write ( *, '(a)' ) ' '

      m = 1

      lc = 1
      ls = 0
      call fser1 ( f1, eps, max, m, c, s, lc, ls )

      exact = ( sin ( m * pi * b ) - sin ( m * pi * a ) ) / ( m * pi )
      write ( *, '(a,2x,i6,2x,g14.6,2x,g14.6)' ) '    1:  ', m, c, exact

      lc = 1
      ls = 0
      call fser1 ( f2, eps, max, m, c, s, lc, ls )

      exact = ( ( cos ( m * pi * b ) 
     &   + m * pi * b * sin ( m * pi * b ) ) 
     &   - ( cos ( m * pi * a ) + m * pi * a * sin ( m * pi * a ) ) ) 
     &   / ( m * pi )**2

      write ( *, '(a,2x,i6,2x,g14.6,2x,g14.6)' ) '    X:  ', m, c, exact

      lc = 1
      ls = 0
      call fser1 ( f3, eps, max, m, c, s, lc, ls )

      exact = ( ( 2.0E+00 * m * pi * b * cos ( m * pi * b ) 
     &   + ( ( m * pi )**2 * b**2 - 2.0E+00 ) * sin ( m * pi * b ) ) 
     &     - ( 2.0E+00 * m * pi * a * cos ( m * pi * a ) 
     &   + ( ( m * pi )**2 * a**2 - 2.0E+00 ) * sin ( m * pi * a ) ) ) 
     &   / ( m * pi )**3

      write ( *, '(a,2x,i6,2x,g14.6,2x,g14.6)' ) '  X*X:  ', m, c, exact

      write ( *, '(a)' ) ' '

      m = 2

      lc = 1
      ls = 0
      call fser1 ( f1, eps, max, m, c, s, lc, ls )

      exact = ( sin ( m * pi * b ) - sin ( m * pi * a ) ) / ( m * pi )
      write ( *, '(a,2x,i6,2x,g14.6,2x,g14.6)' ) '    1:  ', m, c, exact

      lc = 1
      ls = 0
      call fser1 ( f2, eps, max, m, c, s, lc, ls )

      exact = ( ( cos ( m * pi * b ) 
     &   + m * pi * b * sin ( m * pi * b ) ) 
     &   - ( cos ( m * pi * a ) + m * pi * a * sin ( m * pi * a ) ) ) 
     &   / ( m * pi )**2

      write ( *, '(a,2x,i6,2x,g14.6,2x,g14.6)' ) '    X:  ', m, c, exact

      lc = 1
      ls = 0
      call fser1 ( f3, eps, max, m, c, s, lc, ls )

      exact = ( ( 2.0E+00 * m * pi * b * cos ( m * pi * b ) 
     &   + ( ( m * pi )**2 * b**2 - 2.0E+00 ) * sin ( m * pi * b ) ) 
     &     - ( 2.0E+00 * m * pi * a * cos ( m * pi * a ) 
     &   + ( ( m * pi )**2 * a**2 - 2.0E+00 ) * sin ( m * pi * a ) ) ) 
     &   / ( m * pi )**3

      write ( *, '(a,2x,i6,2x,g14.6,2x,g14.6)' ) '  X*X:  ', m, c, exact

      write ( *, '(a)' ) ' '

      m = 3

      lc = 1
      ls = 0
      call fser1 ( f1, eps, max, m, c, s, lc, ls )

      exact = ( sin ( m * pi * b ) - sin ( m * pi * a ) ) / ( m * pi )
      write ( *, '(a,2x,i6,2x,g14.6,2x,g14.6)' ) '    1:  ', m, c, exact

      lc = 1
      ls = 0
      call fser1 ( f2, eps, max, m, c, s, lc, ls )

      exact = ( ( cos ( m * pi * b ) 
     &   + m * pi * b * sin ( m * pi * b ) ) 
     &   - ( cos ( m * pi * a ) + m * pi * a * sin ( m * pi * a ) ) ) 
     &   / ( m * pi )**2

      write ( *, '(a,2x,i6,2x,g14.6,2x,g14.6)' ) '    X:  ', m, c, exact

      lc = 1
      ls = 0
      call fser1 ( f3, eps, max, m, c, s, lc, ls )

      exact = ( ( 2.0E+00 * m * pi * b * cos ( m * pi * b ) 
     &   + ( ( m * pi )**2 * b**2 - 2.0E+00 ) * sin ( m * pi * b ) ) 
     &     - ( 2.0E+00 * m * pi * a * cos ( m * pi * a ) 
     &   + ( ( m * pi )**2 * a**2 - 2.0E+00 ) * sin ( m * pi * a ) ) ) 
     &   / ( m * pi )**3

      write ( *, '(a,2x,i6,2x,g14.6,2x,g14.6)' ) '  X*X:  ', m, c, exact

      return
      end
      function f1 ( x )

c***********************************************************************
c
cc F1 evaluates the integrand factor F(X) = 1.
c
      implicit none

      real f1
      real x

      f1 = 1.0E+00

      return
      end
      function f2 ( x )

c***********************************************************************
c
cc F2 evaluates the integrand factor F(X) = X.
c
      implicit none

      real f2
      real x

      f2 = x

      return
      end
      function f3 ( x )

c***********************************************************************
c
cc F3 evaluates the integrand factor F(X) = X*X.
c
      implicit none

      real f3
      real x

      f3 = x*x

      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'
 
TOMS353_PRB:
  Test ACM algorithm 353, which uses
  Filon quadrature to approximate integrals
  of functions of the form F(X)*COS(M*PI*X)
  or F(X)*SIN(M*PI*X).
 
TEST01
  Use FSER1 to estimate the integrals of
  the form F(X) * COS ( M * PI * X )
 
  Use integrand factors:
  F(X) = 1, X, X*X.
 
               M  Approximate         Exact
 
    1:         1   -0.428771E-07   -0.278275E-07
    X:         1   -0.202642       -0.202642    
  X*X:         1   -0.202642       -0.202642    
 
    1:         2    0.255680E-08    0.278275E-07
    X:         2    0.136169E-07    0.271765E-07
  X*X:         2    0.506606E-01    0.506606E-01
 
    1:         3   -0.215751E-08   -0.253054E-08
    X:         3   -0.225158E-01   -0.225158E-01
  X*X:         3   -0.225158E-01   -0.225158E-01
 
TOMS353_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'
      SUBROUTINE FSER1 ( F, EPS, MAX, M, C, S, LC, LS )

      IMPLICIT NONE

      REAL A
      REAL B
      REAL B1
      REAL B2
      REAL B3
      REAL B4
      REAL B5
      REAL C
      REAL C1
      REAL COSIP
      REAL EPS
      REAL F
      REAL F0
      REAL F1
      REAL FRAC
      REAL G
      REAL H
      INTEGER I
      INTEGER IN
      INTEGER LC
      INTEGER LLC
      INTEGER LLS
      INTEGER LN
      INTEGER LS
      INTEGER M
      INTEGER MAX
      INTEGER MSWTCH
      INTEGER N
      INTEGER NST
      INTEGER NSTOP
      REAL ODCOS
      REAL ODSIN
      REAL P
      REAL PI
      REAL PVT1
      REAL PVT2
      REAL S
      REAL S1
      REAL S1SQ
      REAL SUMCOS
      REAL SUMSIN
      REAL T
      REAL T1
      REAL T2
      REAL TEMP1
      REAL TEMP2
      REAL THA
      REAL TMAX
      REAL TP
      REAL TSQ
      REAL XI
      REAL XM

      PI = 3.1415926535898E+00
      XM = M
C
C  F1 = COS ( M * PI ) TEMPORARY.
C
      F1 = 1.0E+00 - 2.0E+00 * ( M - ( M / 2 ) * 2 )
      F0 = F ( 0.0E+00 )
      F1 = F ( 1.0E+00 ) * F1
C
C  'CIR' WILL BE USED THROUGHOUT THESE COMMENTS TO STAND FOR 'SIN' OR
C  'COS' WHENEVER THOSE TWO SYMBOLS MAY OCCUR.
C  NOW DEFINE SUMCIR OF THE ENDPOINTS.
C
      SUMCOS = ( F1 + F0 ) * 0.5E+00
      SUMSIN = 0.0E+00
      B1 = 2.0E+00 / 3.0E+00
C
C  TMAX IS THE SWITCH-OVER POINT IN THE ANGLE T.
C  OUR ANALYSIS INDICATES THAT TMAX = 1/6 IS THE BEST FOR THE ILIAC II,
C  WHICH HAS A 44 BIT FLOATING POINT MANTISSA.
C
      TMAX = 0.166E+00
C
C  N IS THE NUMBER OF ITERATIONS.  
C
      N = LOG ( 2.0E+00 * XM ) / 0.693E+00
C
C  BOTH TMAX AND N MAY BE CHANGED IF THE MACHINE FOR WHICH THIS
C  ROUTINE IS INTENDED HAS GREATER OR LESS ACCURACY THAN ILIAC II.
C  IF N IS CHANGED, THEN THE CORRESPONDING CHANGES MUST BE MADE
C  IN THE ASSIGNMENTS OF H AND NSTOP.
C
      H = 1.0E+00 / REAL( 2**N )
C
C  H = 1 / 2**N.
C
      NSTOP = 2**N - 1
C
C  NSTOP = 2**N - 1.
C
      T = H * XM
      TP = T * PI
      NST = 1
C ***** trh ******
C Initialized T2 to zero following Unassigned trap @ line259
      T2 = 0.0E0
C ***** end trh ******
      ASSIGN 67 TO MSWTCH
C
C  LLC AND LLS ARE USED BY THE ROUTINE IN COMPUTED-GO-TO STATEMENTS.
C  AS SOON AS LLS AND LLC HAVE BEEN DEFINED, WE CAN USE LS AND LC
C  AS RETURN PARAMETERS (SEE ABOVE).
C
      IF ( LS ) 1, 1, 2
1     LLS = 2
      GO TO 3
2     LLS = 1
      LS = MAX
3     IF ( LC ) 4, 4, 5
4     LLC = 2
      GO TO 7
5     LLC = 1
      LC = MAX
7     LN = 1
C
C  ALL OF THE ABOVE IS EXECUTED ONLY ONCE PER CALL.
C  NOW THE ITERATION BEGINS.
C
10    ODCOS = 0.0E+00
      ODSIN = 0.0E+00
C
C  BEGIN SUMMATION FOR ODCOS AND ODSIN.
C
      DO 65 I = 1, NSTOP, NST
        XI = I
        THA = XI * T
C
C  THA * PI IS THE ANGLE USED IN THIS I-TH TERM.
C
C  CIR(I*T*PI) IS CALCULATED HERE USING THE IDENTITY
C  CIR ( INTEGER MULTIPLE OF PI + FRACTIONAL MULTIPLE OF PI )
C  = COS ( INTEGER * PI ) * CIR ( FRAC * PI )
C  = ( + 0R - ) * CIR ( FRAC * PI ).
C
        FRAC = THA
        IN = THA
        THA = IN
        FRAC = ( FRAC - THA ) * PI
C
C  THA IS A FLOATING POINT INTEGER, FRAC IS THE FRACTIONAL PART * PI.
C
        COSIP = 1 - 2 * ( IN - 2 * ( IN / 2 ) )
        TEMP1 = COSIP * F ( XI * H )
C
C  TEMP1 = COS ( INTEGER PART ) * F ( I * H ).
C
        GO TO ( 50, 55 ) LLS
50      ODSIN = TEMP1 * SIN ( FRAC ) + ODSIN
55      GO TO ( 60, 65 ) LLC
60      ODCOS = TEMP1 * COS ( FRAC ) + ODCOS
65      CONTINUE

      GO TO MSWTCH, ( 67, 70 )
67    NST = 2
C
C  NOW HAVE MADE UP FOR THE FIRST 4 ITERATION STEPS, SO RESET THESE
C  TWO NUMBERS TO LOOK LIKE THE GENERAL CASE.
C
      NSTOP = 2**N
C
C  NSTOP = 2**N, IN CASE YOU CHANGE STARTING VALUE OF N.
C
      ASSIGN 70 TO MSWTCH
      GO TO 92
70    TSQ = TP * TP
      IF ( T - TMAX ) 74, 74, 75
C
C  74 IS THE POWER SERIES FOR SMALL T, 75 IS THE CLOSED FORM USED WITH
C  LARGER VALUES OF T.
C  THE POWER SERIES (WITH 'TN' = TP**N)
C  A = (2/45)*T3 - (2/315)*T5 + (2/4725)*T7
C  B = (2/3) + (2/15)*T2 - (4/105)*T4 + (2/567)*T6 - (4/22275)*T8
C  G = (4/3) - (2/15)*T2 - (1/210)*T4 + (1/11340)*T6
C  THE NEXT TERM IN G IS TOO SMALL.  IT IS (1/997920)*T8.
C
74    A = TP * TSQ * ( 1.0E+00 - TSQ * ( 1.0E+00 - TSQ / 15.0E+00 ) 
     &  / 7.0E+00 ) / 22.5E+00
      B2 = B1 * TSQ * 0.2E+00
      B3 = B2 * TSQ * 2.0E+00 / 7.0E+00
      B4 = B3 * TSQ / 10.8E+00
      B5 = B4 * TSQ * 14.0E+00 / 275.0E+00
      B = B1 + B2 - B3 + B4 - B5
      G = 2.0E+00 * B1 - B2 + B3 / 8.0E+00 - B4 / 40.0E+00
C
C  G = 2*B1 - B2 + B3/8 - B4/40 + 5*B5/896.  IF YOU WANT THE T8
C  TERM INCLUDED IN G.
C
      GO TO 80
C
C  CLOSED FORM OF THE COEFFICIENTS, WHERE AGAIN 'TN' MEANS TP**N.
C  A = 1/TP + COS(TP)*SIN(TP)/T2 - 2*(SIN(TP))**2/T3
C  B = 2*((1+(COS(TP))**2)/T2 - 2*SIN(TP)*COS(TP)/T3).
C  G = 4*( SIN(TP)/T3 - COS(TP)/T2 ).
C
75    IN = T
      TEMP1 = 1 - 2 * ( IN - 2 * ( IN / 2 ) )
      TEMP2 = IN
C
C  TEMP1 IS COS ( INTEGER PART OF TP), TEMP2 IS FRACTIONAL PART OF TP.
C
      TEMP2 = ( T - TEMP2 ) * PI
      S1 = TEMP1 * SIN ( TEMP2 )
C
C  S1 = SIN(TP).
C
      C1 = TEMP1 * COS ( TEMP2 )
C
C  C1 = COS(TP).
C
      P = S1 * C1
      S1SQ = S1 * S1
      A = ((( -2.0E+00 * S1SQ / TP ) + P ) / TP + 1.0E+00 ) / TP
      B = 2.0E+00 * ( ( -2.0E+00 * P / TP ) + 2.0E+00 - S1SQ ) / TSQ
      G = 4.0E+00 * ( S1 / TP - C1 ) / TSQ
80    GO TO ( 81, 85 ) LLS
C
C  HAVE CALCULATED THE COEFFICIENTS.  NOW READY FOR THE INTEGRATION
C  FORMULAS.
C
81    T2 = H * ( A * ( F0 - F1 ) + B * SUMSIN + G * ODSIN )
C
C  ENDT1 IS A SUBROUTINE WHICH CHECKS FOR THE CONVERGENCE OF THE
C  ITERATIONS.  ENDT1 REQUIRES THE PRESENT VALUE TO AGREE WITH THE
C  PREVIOUS VALUE TO WITHIN EPS2, WHERE
C  EPS2 = ( 1.0 + ABS ( PRESENT VALUE ) * EPS.
C  EPS IS SUPPLIED BY THE USER.
C
      CALL ENDT1 ( PVT2, T2, EPS, S, LLS, LN )
      GO TO ( 85, 84 ), LLS
84    LS = N
85    GO TO ( 86, 90 ), LLC
C
C  THIS IS THE COSINE INTEGRAL.
C
86     T1 = H * ( B * SUMCOS + G * ODCOS )
       CALL ENDT1 ( PVT1, T1, EPS, C, LLC, LN )
       GO TO ( 90, 89 ), LLC
89     LC = N
90     LN = 2
C
C  NOW TEST TO SEE IF DONE.
C
      IF ( LLC + LLS - 3 ) 92, 92, 100
92    N = N + 1
C
C  THIS IS THE BEGINNING OF THE ITERATION.
C
      IF ( N - MAX ) 95, 95, 100
95    H = 0.5E+00 * H
      T = 0.5E+00 * T
      TP = 0.5E+00 * TP
      NSTOP = 2 * NSTOP
      SUMSIN = SUMSIN + ODSIN
      SUMCOS = SUMCOS + ODCOS
      GO TO 10
100   S = T2
      C = T1
      RETURN
      END
      SUBROUTINE ENDT1 ( PREVQT, QUANT, EPS, VALUE, L1, L2 )

      IMPLICIT NONE

      REAL EPS
      INTEGER L1
      INTEGER L2
      REAL PREVQT
      REAL QUANT
      REAL REPS
      REAL VALUE

      GO TO ( 29, 20 ), L2
20    REPS = EPS * ( 1.0E+00 + ABS ( QUANT ) )
23    IF ( ABS ( PREVQT - QUANT ) - REPS ) 25, 25, 29
25    VALUE = QUANT
      L1 = 2
      GO TO 30
29    PREVQT = QUANT
30    RETURN
      END
SHAR_EOF
fi # end of overwriting check
cd ..
cd ..
cd ..
#       End of shell archive
exit 0