module numint

  use const_and_precisions, only : wp_, zero, one
  implicit none

contains

  subroutine simpson (n,h,fi,s)     
!   subroutine for integration over f(x) with the simpson rule.  fi:
!   integrand f(x); h: interval; s: integral.  copyright (c) tao pang 1997.
    implicit none
    integer, intent(in) ::  n
    real(wp_), intent(in) :: h
    real(wp_), dimension(n), intent(in) :: fi
    real(wp_), intent(out) ::  s
    integer  :: i
    real(wp_) :: s0,s1,s2

    s  = zero
    s0 = zero
    s1 = zero
    s2 = zero
    do i = 2, n-1, 2
      s1 = s1+fi(i-1)
      s0 = s0+fi(i)
      s2 = s2+fi(i+1)
    end do
    s = h*(s1+4.0_wp_*s0+s2)/3.0_wp_
! if n is even, add the last slice separately     
    if (mod(n,2).eq.0) s = s+h*(5.0_wp_*fi(n)+8.0_wp_*fi(n-1)-fi(n-2))/12.0_wp_
  end subroutine simpson

  subroutine trapezoid(n,xi,fi,s)
!   subroutine for integration with the trapezoidal rule.
!   fi: integrand f(x); xi: abscissa x;
!   s: integral Int_{xi(1)}^{xi(n)} f(x)dx
    implicit none
    integer, intent(in) :: n
    real(wp_), dimension(n), intent(in) :: xi,fi
    real(wp_), intent(out) :: s
    integer :: i

    s  = zero
    do i = 1, n-1
      s = s+(xi(i+1)-xi(i))*(fi(i+1)-fi(i))
    end do
    s = 0.5_wp_*s
  end subroutine trapezoid

  subroutine quanc8(fun,a,b,abserr,relerr,result,errest,nofun,flag)
    implicit none
    real(wp_), intent(in) :: a, b, abserr, relerr
    real(wp_), intent(out) :: result, errest, flag
    integer, intent(out) :: nofun
!
!   estimate the integral of fun(x) from a to b
!   to a user provided tolerance.
!   an automatic adaptive routine based on
!   the 8-panel newton-cotes rule.
!
!   input ..
!
!   fun     the name of the integrand function subprogram fun(x).
!   a       the lower limit of integration.
!   b       the upper limit of integration.(b may be less than a.)
!   relerr  a relative error tolerance. (should be non-negative)
!   abserr  an absolute error tolerance. (should be non-negative)
!
!   output ..
!
!   result  an approximation to the integral hopefully satisfying the
!           least stringent of the two error tolerances.
!   errest  an estimate of the magnitude of the actual error.
!   nofun   the number of function values used in calculation of result.
!   flag    a reliability indicator.  if flag is zero, then result
!           probably satisfies the error tolerance.  if flag is
!           xxx.yyy , then  xxx = the number of intervals which have
!           not converged and  0.yyy = the fraction of the interval
!           left to do when the limit on  nofun  was approached.
!
    real(wp_) :: w0,w1,w2,w3,w4,area,x0,f0,stone,step,cor11,temp
    real(wp_) :: qprev,qnow,qdiff,qleft,esterr,tolerr
    real(wp_), dimension(31) :: qright
    real(wp_), dimension(16) :: f,x
    real(wp_), dimension(8,30) :: fsave,xsave
    integer :: levmin,levmax,levout,nomax,nofin,lev,nim,i,j

    interface
       function fun(x)
         use const_and_precisions, only : wp_
         implicit none
         real(wp_), intent(in) :: x
         real(wp_) :: fun
       end function fun
    end interface
!
!   ***   stage 1 ***   general initialization
!   set constants.
!
    levmin = 1
    levmax = 30
    levout = 6
    nomax = 5000
    nofin = nomax - 8*(levmax-levout+2**(levout+1))
!
!   trouble when nofun reaches nofin
!
    w0 =   3956.0_wp_ / 14175.0_wp_
    w1 =  23552.0_wp_ / 14175.0_wp_
    w2 =  -3712.0_wp_ / 14175.0_wp_
    w3 =  41984.0_wp_ / 14175.0_wp_
    w4 = -18160.0_wp_ / 14175.0_wp_
!
!   initialize running sums to zero.
!
    flag = zero
    result = zero
    cor11  = zero
    errest = zero
    area   = zero
    nofun = 0
    if (a .eq. b) return
!
!   ***   stage 2 ***   initialization for first interval
!
    lev = 0
    nim = 1
    x0 = a
    x(16) = b
    qprev  = zero
    f0 = fun(x0)
    stone = (b - a) / 16.0_wp_
    x(8)  =  (x0  + x(16)) / 2.0_wp_
    x(4)  =  (x0  + x(8))  / 2.0_wp_
    x(12) =  (x(8)  + x(16)) / 2.0_wp_
    x(2)  =  (x0  + x(4))  / 2.0_wp_
    x(6)  =  (x(4)  + x(8))  / 2.0_wp_
    x(10) =  (x(8)  + x(12)) / 2.0_wp_
    x(14) =  (x(12) + x(16)) / 2.0_wp_
    do j = 2, 16, 2
       f(j) = fun(x(j))
    end do
    nofun = 9
!
!   ***   stage 3 ***   central calculation
!   requires qprev,x0,x2,x4,...,x16,f0,f2,f4,...,f16.
!   calculates x1,x3,...x15, f1,f3,...f15,qleft,qright,qnow,qdiff,area.
!
    do
      do
        x(1) = (x0 + x(2)) / 2.0_wp_
        f(1) = fun(x(1))
        do j = 3, 15, 2
          x(j) = (x(j-1) + x(j+1)) / 2.0_wp_
          f(j) = fun(x(j))
        end do
        nofun = nofun + 8
        step = (x(16) - x0) / 16.0_wp_
        qleft = (w0*(f0 + f(8)) + w1*(f(1)+f(7)) + w2*(f(2)+f(6)) &
           + w3*(f(3)+f(5)) + w4*f(4)) * step
        qright(lev+1)=(w0*(f(8)+f(16))+w1*(f(9)+f(15))+w2*(f(10)+f(14)) &
           + w3*(f(11)+f(13)) + w4*f(12)) * step
        qnow = qleft + qright(lev+1)
        qdiff = qnow - qprev
        area = area + qdiff
!
!   ***   stage 4 *** interval convergence test
!
        esterr = abs(qdiff) / 1023.0_wp_
        tolerr = max(abserr,relerr*abs(area)) * (step/stone)
        if (lev .ge. levmin) then
!
!   ***   stage 6   ***   trouble section
!   number of function values is about to exceed limit.
!
          if (lev .ge. levmax) then
!
!   current level is levmax.
!
            flag = flag + one
            exit
          end if
          if (nofun .gt. nofin) then
            nofin = 2*nofin
            levmax = levout
            flag = flag + (b - x0) / (b - a)
            exit
          end if
          if (esterr .le. tolerr) exit
        end if
!
!   ***   stage 5   ***   no convergence
!   locate next interval.
!
        nim = 2*nim
        lev = lev+1
!
!   store right hand elements for future use.
!
        do i = 1, 8
          fsave(i,lev) = f(i+8)
          xsave(i,lev) = x(i+8)
        end do
!
!   assemble left hand elements for immediate use.
!
        qprev = qleft
        do i = 1, 8
          j = -i
          f(2*j+18) = f(j+9)
          x(2*j+18) = x(j+9)
        end do
      end do
!
!   ***   stage 7   ***   interval converged
!   add contributions into running sums.
!
      result = result + qnow
      errest = errest + esterr
      cor11  = cor11  + qdiff / 1023.0_wp_
!
!   locate next interval.
!
      do
        if (nim .eq. 2*(nim/2)) exit
        nim = nim/2
        lev = lev-1
      end do
      nim = nim + 1
      if (lev .le. 0) exit
!
!   assemble elements required for the next interval.
!
      qprev = qright(lev)
      x0 = x(16)
      f0 = f(16)
      do i = 1, 8
        f(2*i) = fsave(i,lev)
        x(2*i) = xsave(i,lev)
      end do
    end do
!
!   ***   stage 8   ***   finalize and return
!
    result = result + cor11
!
!   make sure errest not less than roundoff level.
!
    if (errest .eq. zero) return
    do
      temp = abs(result) + errest
      if (temp .ne. abs(result)) return
      errest = 2.0_wp_*errest
    end do
  end subroutine quanc8

end module numint