'**************************************************************
'* Program to demonstrate integration of a real function F(x) *
'*                    by Simpson's method                     *
'* ---------------------------------------------------------- *
'* REFERENCE: "Mathematiques en Turbo-Pascal (Part 1) By      *
'*             Marc Ducamp and Alain Reverchon, Eyrolles,     *
'*             Paris, 1987" [BIBLI 03].                       *
'* ---------------------------------------------------------- *
'* SAMPLE RUN:                                                *
'* (Integrate sin(x) from x=0 to x=1)                         *
'*                                                            *
'* Integral of a function F(X) by Simpson's method            *
'*                                                            *
'* Input begin and end values for x variable:                 *
'*                                                            *
'*         X0 = 0                                             *
'*         X1 = 1                                             *
'*                                                            *
'* Number of integration steps: 100                           *
'*                                                            *
'*                                                            *
'* Value of integral:  .4596976941334565                      *
'*                                                            *
'*                                                            *
'*                              Basic version by J-P Moreau.  *
'*                                   (www.jpmoreau.fr)        *
'**************************************************************
'PROGRAM Test_Simpson
DEFDBL A-H, O-Z
DEFINT I, N

' double  x0        begin x value
' double  x1        end   x value
' integer nstep     number of integration steps
' double  result    result of integral

  CLS
  PRINT
  PRINT " Integral of a function F(X) by Simpson's method"
  PRINT
  PRINT " Input begin and end values for x variable:"
  PRINT
  INPUT "         x0 = ", x0
  INPUT "         x1 = ", x1
  PRINT
  INPUT " Number of integration steps: ", nstep

  GOSUB 1000 'call Integral_Simpson(x0,x1,nstep,result)

  PRINT
  PRINT
  PRINT " Value of integral: "; result
  PRINT
  PRINT

END

'Given function to integrate
500 'double FUNC(x)
  FUNC = SIN(x)
RETURN

'*******************************************************
'* Integral of a function FUNC(X) by Simpson's method  *
'* --------------------------------------------------- *
'* INPUTS:                                             *
'*          a      begin value of x variable           *
'*          b      end value of x variable             *
'*          n      number of integration steps         *
'*                                                     *
'* OUTPUT:                                             *
'*          res    the integral of FUNC(X) from a to b *
'*                                                     *
'*******************************************************
1000 'Subroutine Integral_Simpson(a, b, n, res)
    n = nstep
    a = x0: b = x1
    xstep = (b - a) / 2 / n
    x = a: GOSUB 500: r = FUNC
    x = b: GOSUB 500
    result = (r + FUNC) / 2#
    FOR i = 1 TO 2 * n - 1
      x = a + i * xstep: GOSUB 500
      r = FUNC
      IF (i MOD 2) <> 0 THEN
	result = result + r + r
      ELSE
	result = result + r
      END IF
    NEXT i
    result = result * xstep * 2# / 3#
RETURN

' End of file tsimpson.bas