!**************************************************************
!* 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:  0.459697694133457                      *
!*                                                            *
!*                                                            *
!*                                F90 version by J-P Moreau.  *
!*                                    (www.jpmoreau.fr)       *
!**************************************************************
PROGRAM Test_Simpson

  real*8  x0        !begin x value
  real*8  x1        !end   x value
  integer nstep     !number of integration steps
  real*8  result    !result of integral

  print *,' '
  print *,' Integral of a function F(X) by Simpson''s method'
  print *,' '
  print *,' Input begin and end values for x variable:'
  print *,' '
  write(*,10,advance='no'); read *, x0
  write(*,20,advance='no'); read *, x1
  print *,' '
  write(*,30,advance='no'); read *, nstep

  call Integral_Simpson(x0,x1,nstep,result)

  print *,' '
  print *,' '
  print *,' Value of integral: ', result
  print *,' '
  print *,' '

10 format('         X0 = ') 
20 format('         X1 = ')
30 format('  Number of integration steps: ') 

END

! Given function to integrate
real*8 Function FUNC(x)
  real*8 x
  FUNC = dsin(x)  
end


!*******************************************************
!* 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 *
!*                                                     *
!*******************************************************
Subroutine Integral_Simpson(a, b, n, res)
    real*8  a,b,res
    integer n, i
	real*8  step, r
    step=(b-a)/2/n
    r=FUNC(a)
    res=(r+FUNC(b))/2
    do i=1, 2*n-1
      r=FUNC(a+i*step);
      if (MOD(i,2).ne.0) then
	    res = res + r+r
      else 
	    res = res + r
      end if
    end do
    res = res * step*2.d0/3.d0
return
end


! End of file tsimpson.f90