c______________________________________________MainProgram_______________
c
c Calculation of an integral using Simpson's Composite Formula (Adaptive)
c

      implicit real*8 (a-h,o-z)
      logical cnvgc
      external fprueba        !name of the subintegrating function

      a=-4.d+00               !lower integration bound
      b= 4.d+00               !upper integration bound

      n_ini=1                 !We start with 1 subinterval (1+1  points)
      nitr_max=10             !maximum number of iterations

      epsl=1.d-20             !maximum admissible absolute error of the integral
      rltv=1.d-12             !maximum permissible relative error of the integral

      call IntFCSA(a,b,fprueba,n_ini,nitr_max,epsl,rltv,n,vint,cnvgc)

      write(6,*) ' n     = ',n      !number of sub-intervals used,
      write(6,*) ' vint  = ',vint   !numerical approx. of the value of the integral
      write(6,*) ' cnvgc = ',cnvgc  !the errors are (T) or are not (F) admissible

      end

c______________________________________________fprueba(x)=1/(1+x*x)___________
c     Subintegrand function that we use as an example

      real*8 function fprueba(x)
      implicit real*8 (a-h,o-z)

      fprueba=1.d+00/(1+x*x)

      return
      end

c______________________________________________IntFCSA________________________
!             /b
! Calculus of | f(x)dx  using the Composite Simpson Formula (Adaptive)
!             /a
!
! Otros datos:
!
!     n_ini     = initial number of subintervals <=> 2*n_ini+1 points
!     niter_max = maximum number of iterations
!     epsl      = maximum absolute error admissible for the value of the integral
!     rltv      = maximum relative error admissible for the value of the integral
!
! Resultados:
!
!     n         = number of subintervals that have been used <=> 2*n+1 points
!     vint      = numerical approach of the value of the integral
!     cnvgc     = logical variable that indicates if the final errors
!                 are (.TRUE.) or they are not (.FALSE.) admissible
!
! Notas:
!
!   * At each iteration the number of subintervals is doubled.
!   * To estimate the error of the numerical approximation in each iteration
! and to calculate the final value of the integral, a
! RICHARDSON EXTRAPOLATION is performed.
!_____________________________________________________________________________

      Subroutine IntFCSA(a,b,f,n_ini,nitr_max,epsl,rltv,n,vint,cnvgc)
      implicit real*8 (a-h,o-z)
      logical cnvgc

c     Simpson's Composite Formula with n_ini subintervals
        nitr=0
        n=n_ini
        h=(b-a)/dfloat(2*n)  !distance between points
        t1=f(a)+f(b)         !terms to be multiplied times 1
        t4=0.d+00            !terms to be multiplied times 4
        do i=1,n
          x=a+dfloat(2*i-1)*h
          t4=t4+f(x)
        enddo
        t2=0.d+00            !terms to be multiplied times 2
        do i=1,n-1
          x=a+dfloat(2*i)*h
          t2=t2+f(x)
        enddo
        vint1=(t1+4.d+00*t4+2.d+00*t2)*h/3.d+00  !Composite Simpson formula
        vint=vint1

c     Adaptive Composite Simpson formula
      cnvgc=.false.
      do while(.not.cnvgc.and.nitr.lt.nitr_max)
        nitr=nitr+1
        n=n*2                !We duplicate the number of subintervals
        h=h/2.d+00           !reduce the distance between points to the half
        t2=t2+t4             !terms to be multiplied now times 2
        t4=0.d+00
        do i=1,n             !terms to be multiplied now times 4
          x=a+dfloat(2*i-1)*h
          t4=t4+f(x)
        enddo
        vint2=(t1+4.d+00*t4+2.d+00*t2)*h/3.d+00 !Composite Simpson formula
        e2=(vint2-vint1)/15.d+00
        vint=vint2+e2                           !Richardson Extrapolation
        e2mx=abs(vint)*rltv
        if (epsl.gt.e2mx) e2mx=epsl             !Admissible error bound
        cnvgc=(abs(e2).le.e2mx)
        vint1=vint2
      enddo

      return
      end