!****************************************************************
!* Integration of a discrete function F(x) by Simpson's method  *
!* ------------------------------------------------------------ *
!* SAMPLE RUN                                                   *
!*                                                              *
!* (Find integral of digitalized function sin(x) from x=0 to    *
!*  x=1).                                                       *
!*                                                              *
!* Integration of a discrete real function F(x)                 *
!*             by Simpson's Method                              *
!*                                                              *
!* Number of points: 25                                         *
!* Begin x         : 0                                          *
!* End x           : 1                                          *
!*                                                              *
!* Value of integral: 0.4596978                                 *
!*                                                              *
!*                                                              *
!*                               F90 version by J-P Moreau.     *
!*                                   (www.jpmoreau.fr)          *
!****************************************************************
PROGRAM Test_Discreet_Simpson

parameter(NSIZE=100)


real    Xi(NSIZE), Yi(NSIZE)
integer i, npoints
real    dx,result,x,x1,x2

  print *,'    Integration of a discrete real function F(x)'
  print *,'                by Simpson''s Method'
  print *,' '
  write(*,10,advance='no'); read *, npoints
  write(*,20,advance='no'); read *, x1
  write(*,21,advance='no'); read *, x2
  ! build up test discreet curve F(x)=sin(x)
  dx=(x2-x1)/(npoints-1)
  x=x1-dx
  do i=0, npoints-1
    x = x + dx
    Xi(i)=x; Yi(i)=sin(x)
  end do

  ! call integration function
  result = DiscreetIntegral(npoints,Xi,Yi)

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

10 format('  Number of points: ')
20 format('  Begin x         : ')
21 format('  End x           : ')

END


!***********************************************************
!* This function returns the integral value of a discrete  *
!* real function F(x) defined by n points from x1 to xn.   *
!* ------------------------------------------------------- *
!* INPUTS:                                                 *
!*          n      Number of points (xi, yi)               *
!*          X      Table of n xi values in ascending order *
!*          Y      Table of n yi values                    *
!*                                                         *
!* OUTPUT          The function returns the integral of    *
!*                 F(x) from x1 to xn.                     *
!* ------------------------------------------------------- *
!* Ref.: "Mathematiques en Turbo Pascal by Alain Reverchon *
!*        and Marc Ducamp, Armand Colin Editor, Paris,     *
!*        1990".                                           *
!***********************************************************
real Function DiscreetIntegral(n, X, Y)
  parameter(NSIZE=100)
  real X(NSIZE), Y(NSIZE)
  real h,k,j,r
  r=0.
  i=0
  do while (i<=n-3)
    i1=i+1; i2=i+2; h=X(i1)-X(i)
    k=X(i2)-X(i1); j=h+k
    if (h==0) then
	  DiscreetIntegral=0.
	  return
    end if
    r = r + ((h+h-k)*Y(i)+(j*j*Y(i1)+h*(k+k-h)*Y(i2))/k)*j/6/h
    i=i2
  end do
  if (MOD(n,2)==0)  &
    r = r + (X(n-1)-X(n-2))*(Y(n-1)+Y(n-2))/2.
  DiscreetIntegral = r
end

! end of file disinteg.f90