c______________________________________________ProgramaPrincipal_______________
c
c Calculo de una integral mediante la Formula Compuesta de Simpson (Adaptativa)
c

      implicit real*8 (a-h,o-z)
      logical cnvgc
      external fprueba        !nombre de la funcion subintegrando

      a=-4.d+00               !extremo izquierdo de integracion
      b= 4.d+00               !extremo derecho   de integracion

      n_ini=1                 !comenzamos con 1 subintervalo (1+1  puntos)
      nitr_max=10             !numero maximo de iteraciones

      epsl=1.d-20             !error absoluto maximo admisible de la integral
      rltv=1.d-12             !error relativo maximo admisible de la integral

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

      write(6,*) ' n     = ',n      !numero de subintervalos utilizados,
      write(6,*) ' vint  = ',vint   !aprox. numerica del valor de la integral
      write(6,*) ' cnvgc = ',cnvgc  !los errores son (T) o no (F) admisibles

      end

c______________________________________________fprueba(x)=1/(1+x*x)___________
c     Función subintegrando que utilizamos como ejemplo

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

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

      return
      end

c______________________________________________IntFCSA________________________
!             /b
! Cálculo de  | f(x)dx  usando la FORMULA COMPUESTA DE SIMPSON (ADAPTATIVA)
!             /a
!
! Otros datos:
!
!     n_ini     = numero inicial de subintervalos <=> 2*n_ini+1 puntos
!     niter_max = numero maximo de iteraciones
!     epsl      = error absoluto maximo tolerable en el valor de la integral
!     rltv      = error relativo maximo tolerable en el valor de la integral
!
! Resultados:
!
!     n         = numero de subintervalos que se han usado <=> 2*n+1 puntos
!     vint      = aproximacion numerica del valor de la integral
!     cnvgc     = variable lógica que indica si los errores finales
!                 son (.TRUE.) o no son (.FALSE.) admisibles
!
! Notas:
!
!   * En cada iteracion se duplica el numero de subintervalos.
!   * Para estimar el error de la aproximacion numerica en cada iteracion
!     y para calcular el valor final de la integral se realiza una
!     EXTRAPOLACION DE RICHARDSON.
!_____________________________________________________________________________

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

c     Formula Compuesta de Simpson con n_ini subintervalos
        nitr=0
        n=n_ini
        h=(b-a)/dfloat(2*n)  !distancia entre puntos
        t1=f(a)+f(b)         !terminos que van multiplicados por 1
        t4=0.d+00            !terminos que van multiplicados por 4
        do i=1,n
          x=a+dfloat(2*i-1)*h
          t4=t4+f(x)
        enddo
        t2=0.d+00            !terminos que van multiplicados por 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  !Formula de Simpson Compuesta
        vint=vint1

c     Formula Compuesta de Simpson Adaptativa
      cnvgc=.false.
      do while(.not.cnvgc.and.nitr.lt.nitr_max)
        nitr=nitr+1
        n=n*2                !duplicamos el numero de aubintervalos
        h=h/2.d+00           !reducimos la distancia entre puntos a la mitad
        t2=t2+t4             !terminos que ahora van multiplicados por 2
        t4=0.d+00
        do i=1,n             !terminos que ahora van multiplicados por 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 !Formula de Simpson Compuesta
        e2=(vint2-vint1)/15.d+00
        vint=vint2+e2                           !Extrapolacion de Richardson
        e2mx=abs(vint)*rltv
        if (epsl.gt.e2mx) e2mx=epsl             !Error absoluto admisible
        cnvgc=(abs(e2).le.e2mx)
        vint1=vint2
      enddo

      return
      end