c______________________________________________ProgramaPrincipal_______________
c
c Calculo de una integral mediante la Formula Compuesta del Trapecio(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=20             !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 IntFCTA(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

      read(5,*)               !espera a que se pulse return para finalizar

      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______________________________________________IntFCTA________________________
!             /b
! Cálculo de  | f(x)dx  usando la FORMULA COMPUESTA DEL TRAPECIO (ADAPTATIVA)
!             /a
!
! Otros datos:
!
!     n_ini     = numero inicial de subintervalos <=> n_ini+1 puntos
!     nitr_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 <=> 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 IntFCTA(a,b,f,n_ini,nitr_max,epsl,rltv,n,vint,cnvgc)
      implicit real*8 (a-h,o-z)
      logical cnvgc

c     Formula Compuesta del Trapecio con n_ini subintervalos
        nitr=0
        n=n_ini
        h=(b-a)/dble(n)      !distancia entre puntos
        t1=f(a)+f(b)         !terminos que van multiplicados por 1
        t2=0.d+00            !terminos que van multiplicados por 2
        do i=1,n-1
          x=a+dble(i)*h
          t2=t2+f(x)
        enddo
        vint1=(t1+2.d+00*t2)*h/2.d+00    !Formula del Trapecio Compuesta
        vint=vint1

c     Formula Compuesta del Trapecio Adaptativa
      cnvgc=.false.
      do while(.not.cnvgc.and.nitr.lt.nitr_max)
        nitr=nitr+1
        n=n*2                !duplicamos el numero de subintervalos
        h=h/2.d+00           !reducimos la distancia entre puntos a la mitad
        tn=0.d+00            !nuevos terminos que van multiplicados por 2
        do i=1,n-1,2
          x=a+dble(i)*h
          tn=tn+f(x)
        enddo
        t2=t2+tn
        vint2=(t1+2.d+00*t2)*h/2.d+00    !Formula del Trapecio Compuesta
        e2=(vint2-vint1)/3.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