!
!  Numerical integration of f(x) on [a,b] 
!  Methods:
!  1. adaptive non-recursive Simpson 
!  2. adaptive non-recursive Gauss-Kronrod 7-15
!  3. adaptive non-recursive Newton-Cotes Quanc8     
!
! Last update: February 2022

! use this modules to count number of function calls
! handy for recursive procedures
module fcalls
  implicit none
  integer count
end module fcalls


    Program main
!===============================================================
! Integration of a function using various integration methods 
!===============================================================
use fcalls
implicit none
double precision f, a, b, eps, sum, sum1, sum2, sum3, sum4, sum5
double precision qsum, abserr, errest, flag
double precision, parameter:: pi = 3.1415926535897932
external f
integer i, c1, c2, c3, c4, c5
integer nfun, nofun
double precision g7, k15, gauss, kronrod

a = 0.0
b = pi

eps  = 1.0e-12

abserr = 0.0
i = 0

  call simpson5(f,a,b,eps,sum,errest,nfun)
  sum2 = sum
  write (*,102) sum2, eps, errest, nfun

!  call gauss715a(f,a,b,eps,sum,errest,nfun)
!  sum3 = sum
!  write (*,103) sum3, eps, errest, nfun

  call quanc8(f,a,b,abserr,eps,qsum,errest,nofun,flag)
  sum1 = qsum
  write (*,101) sum1, eps, errest, nofun

  call GKinteg(f,a,b,eps,sum,errest,nfun,1)
  sum3 = sum
  write (*,103) sum3, eps, errest, nfun

  call GKinteg(f,a,b,eps,sum,errest,nfun,2)
  sum3 = sum
  write (*,104) sum3, eps, errest, nfun

!  call G10K21(f,a,b,gauss,kronrod)
!  write (*,101) gauss, kronrod

100 format('*** Version 20 July 2012 ***',/)
101 format(' Integral Quanc8  = ',f12.8,/,' Tolerance        = ',f12.8,/, &
           ' Tol. estimated   = ',f12.8,/,' Function calls   = ',i12,/)
102 format(' Adaptive Simpson = ',f12.8,/,' Tolerance        = ',f12.8,/, &
           ' Tol. estimated   = ',f12.8,/,' Function calls   = ',i12,/)
103 format(' Adaptive G7K15   = ',f12.8,/,' Tolerance        = ',f12.8,/, &
           ' Tol. estimated   = ',f12.8,/,' Function calls   = ',i12,/)
104 format(' Adaptive G10K21  = ',f12.8,/,' Tolerance        = ',f12.8,/, &
           ' Tol. estimated   = ',f12.8,/,' Function calls   = ',i12,/)
end


  Function f(x)
!----------------------------------------
! Function for integration
!----------------------------------------
implicit none
double precision f, x
double precision, parameter:: pi=3.1415926535897932
  f = sin(x)
!   f = cos(cos(x) + 3.0*sin(x) + 2.0*cos(2.0*x) + 3.0*sin(2.0*x) + 3.0*cos(3.0*x) )
!  f = 2.0/(2.0 + sin(10.0*pi*x))
!  f = (100.0/x**2)*sin(10.0/x)
!  f = exp(x)*sin(pi*x)
!  f = 1.0/((x-0.3)**2+0.01) + 1.0/((x-0.9)**2+0.04) - 6.0
!  f = x*sin(30.0*x)*cos(50.0*x) 
!  f = 1.0/sqrt(x)
return
end

 Subroutine simpson5(f,a,b,eps,sum,errest,nfun)
!==========================================================
! Integration of f(x) on [a,b]
! Method: adaptive non-recursive Simpson rule  
! written by: Alex Godunov (July 2012)
!----------------------------------------------------------
! IN:
! f    - Function to integrate
! a	   - Lower limit of integration
! b	   - Upper limit of integration
! eps  - Error tolerance abs(s1+s2-s0)<=15*eps (using 10*eps recommended)
!
! OUT:
! sum    - the result of integration
! errest - estimated tolerance
! nfun   - the number of function values used in calculation
!
! PARAMETERS:
! im - max number of levels (i.e. the smallest interval (b-a)/2**imax)
! nmax - max number of function calls
!
! Comment:
! reference points for two Simpson's intervals
! numbers      0   1   2   3   4
! in the code  a---1---m---3---b
!
! 14.07.2012 
! so far nmax has not been used
! think about relative error estimation (see the condition)
!==========================================================
implicit none
double precision f, a, b, eps, sum, errest
integer nfun

integer, parameter :: im=32, nmax=5000
double precision tol(im), x(im), h(im), fa(im), fm(im), fb(im), s(im)
double precision step, err, xfail
double precision x0, f0, f1, f2, f3, f4, s0, s1, s2
integer il(im)
integer i, imax, deep

! *** stage 1 ***
! initialization for level 1 (top level)

sum = 0.0
errest = 0.0
xfail = 0.0
i = 1
imax = im

x(1)   = a
h(1)   = (b-a)/2.0
fa(1)  = f(a)
fm(1)  = f(a+h(1))
fb(1)  = f(b)
tol(1) = 15.0*eps    
il(1)  = 1

! Simpson's method for [a,b]
s(1) = h(1)*(fa(1)+4*fm(1)+fb(1))/3.0
nfun = 3

! *** stage 2 ***
! main part: adaptive integration

do while (i > 0)
 
! *** stage 2a ***
! calculate function values at h/2 and 3h/2
  f1 = f(x(i) +     h(i)/2.0)
  f3 = f(x(i) + 3.0*h(i)/2.0)
  nfun = nfun + 2

! Simpson's integrals for the left and right intervals
  s1 = h(i)*(fa(i)+4.0*f1+fm(i))/6.0
  s2 = h(i)*(fm(i)+4.0*f3+fb(i))/6.0

! *** stage 2b ***
! save data at this level
  x0 = x(i)
  f0 = fa(i)
  f2 = fm(i)
  f4 = fb(i)
  step = h(i)
  err  = tol(i)
  s0   = s(i)
  deep = il(i)

! *** stage 2c ***
! the current level has been computed 
  i=i-1

! *** stage 2d ***
! local condition for convergence
  if(abs(s1+s2-s0) <= err ) then
    sum = sum + (s1+s2)
	errest = errest + abs(s1+s2-s0)/15.0
    else
! *** stage 2e *** 
!   check if the code can continue to subdivide intervals
      if( deep >= imax ) then 
      xfail = xfail + step
      else
! *** stage 2f ***
!	  make smaller intervals (i.e. h=h/2)
!     data for the right subinterval
        i = i+1
        x(i)  = x0 + step
        fa(i) = f2
        fm(i) = f3
        fb(i) = f4
        h(i)  = step/2.0
        tol(i)= err/2.0
        s(i)  = s2
        il(i) = deep + 1
!     data for the left subinterval
        i = i+1
        x(i)  = x0
        fa(i) = f0
        fm(i) = f1
        fb(i) = f2
        h(i)  = h(i-1)
        tol(i)= tol(i-1)
        s(i)  = s1
        il(i) = il(i-1)
	end if
  end if
end do

return
end subroutine simpson5


subroutine G7K15(f,a,b,g7,k15)
!==============================================================
! Integration of f(x) on [a,b]
! Method: Gauss7 and Kronrod15 rules applied to whole interval
! written by: Alex Godunov (July 2012)
!--------------------------------------------------------------
! IN:
! f   - Function to integrate
! a	  - Lower limit of integration
! b	  - Upper limit of integration
! !
! OUT:
! g7  - integral Gauss 7 point rule
! k15 - integral Kronrod 15 point rule
!==============================================================
implicit none
double precision f, a, b, g7, k15
double precision h, c

! note: the x points are symmetrical relative to 0
! gauss 7 points are x: -6, -4, -2, 0, 2, 4, 6
double precision x(0:7), wg(0:3), wk(0:7)
double precision f0, fl(7), fr(7)
data  x/0.000000000000000d0, 0.207784955007898d0, 0.405845151377397d0, 0.586087235467691d0, &
        0.741531185599394d0, 0.864864423359769d0, 0.949107912342758d0, 0.991455371120813d0/
data wg/0.417959183673469d0, 0.381830050505119d0, 0.279705391489277d0, 0.129484966168870d0/
data wk/0.209482141084728d0, 0.204432940075298d0, 0.190350578064785d0, 0.169004726639267d0, &
        0.140653259715525d0, 0.104790010322250d0, 0.063092092629979d0, 0.022935322010529d0/
integer i, j
external f

!*** initialization ***
g7  = 0.0
k15 = 0.0

h = (b-a)/2.0
c = (b+a)/2.0

! evaluate function f at points x using the symmetry for -x and +x
! the gauss interval [-1,1] has to be rescaled to [a,b]
! integral {f(x)dx} on [a,b] = (b-a)/2*integral{f([(b-a)/2]*y + (b+a)/2)}dy on [-1,1]
f0 = f(c)
do i = 1,7
  fl(i) = f(h*(-1.0)*x(i)+c)
  fr(i) = f(h*x(i)+c)
end do

! evaluate Gauss 7
  g7 = wg(0)*f0
  do i=1,3
    g7 = g7 + wg(i)*(fl(2*i)+fr(2*i))
  end do
  g7 = g7*h

! evaluate Kronrod 15
  k15 = wk(0)*f0
  do i=1,7
    k15 = k15 + wk(i)*(fl(i)+fr(i))
  end do
  k15 = k15*h

return 
end subroutine G7K15



subroutine G10K21(f,a,b,gauss,kronrod)
!==============================================================
! Integration of f(x) on [a,b]
! Method: Gauss10 and Kronrod21 rules applied to whole interval
! written by: Alex Godunov (July 2012)
!--------------------------------------------------------------
! IN:
! f   - Function to integrate
! a	  - Lower limit of integration
! b	  - Upper limit of integration
! !
! OUT:
! gauss   - integral Gauss 10 point rule
! kronrod - integral Kronrod 21 point rule
!==============================================================
implicit none
double precision f, a, b, gauss, kronrod
double precision h, c, f1

double precision x(11), wk(11), wg(5)
double precision fleft(11), fright(11)

! Kronrod points: 1  2  3  4  5  6  7  8  9 10 11
! Gauss points       1     2     3     4     5
data  x/0.0000000000000000000000000d0, 0.1488743389816312108848260d0, &
        0.2943928627014601981311266d0, 0.4333953941292471907992659d0, &
        0.5627571346686046833390001d0, 0.6794095682990244062343274d0, &
        0.7808177265864168970637176d0, 0.8650633666889845107320967d0, &
	    0.9301574913557082260012072d0, 0.9739065285171717200779640d0, &
        0.9956571630258080807355273d0/
data wk/0.1494455540029169056649365d0, 0.1477391049013384913748415d0, &
        0.1427759385770600807970943d0, 0.1347092173114733259280540d0, &
        0.1234919762620658510779581d0, 0.1093871588022976418992106d0, &
        0.0931254545836976055350655d0, 0.0750396748109199527670431d0, &
        0.0547558965743519960313813d0, 0.0325581623079647274788190d0, &
        0.0116946388673718742780644d0/
data wg/0.2955242247147528701738930d0, 0.2692667193099963550912269d0, &
        0.2190863625159820439955349d0, 0.1494513491505805931457763d0, &
		0.0666713443086881375935688d0/
integer i
external f

!*** initialization ***
gauss   = 0.0
kronrod = 0.0

h = (b-a)/2.0
c = (b+a)/2.0

! evaluate function f at points x using the symmetry for -x and +x
! the gauss interval [-1,1] has to be rescaled to [a,b]
! integral {f(x)dx} on [a,b] = (b-a)/2*integral{f([(b-a)/2]*y + (b+a)/2)}dy on [-1,1]
f1 = f(c)
do i = 2,11
  fleft(i)  = f(h*(-1.0)*x(i)+c)
  fright(i) = f(h*x(i)+c)
end do

! evaluate Gauss 10
  do i=1,5
    gauss = gauss + wg(i)*(fleft(2*i)+fright(2*i))
  end do
  gauss = gauss*h

! evaluate Kronrod 21
  kronrod = wk(1)*f1
  do i=2,11
    kronrod = kronrod + wk(i)*(fleft(i)+fright(i))
  end do
  kronrod = kronrod*h

return 
end subroutine G10K21

 Subroutine GKinteg(f,a,b,eps,sum,errest,nfun,key)
!===============================================================
! Integration of f(x) on [a,b]
! Method: adaptive non-recursive Gauss-Kronrod  
! written by: Alex Godunov (July 2012)
!---------------------------------------------------------------
! IN:
! f    - Function to integrate
! a	   - Lower limit of integration
! b	   - Upper limit of integration
! eps  - Error tolerance (200.0*abs(k15-g7))**1.5 <=err
! key  - select quadrature
!        =1  Gauss  7 - Kronrod 15
!        =2  Gauss 10 - Kronrod 21
!
! OUT:
! sum    - the result of integration
! errest - estimated tolerance
! nfun   - the number of function values used in calculation
!
! PARAMETERS:
! im - max number of levels: the smallest interval (b-a)/2**imax)
!
! Comment:
! the function calls subroutine G7K15
!
! 17.07.2012 
! think about relative error estimation (see the condition)
!===============================================================
implicit none
double precision f, a, b, eps, errest
integer, parameter :: im=32
double precision tol(im), x(im), h(im)
double precision x0, step, err, sum, xfail
double precision gauss, kronrod
integer il(im)
integer i, imax, deep, nfun, key
external f

! *** stage 1 ***
! initialization for level 1 (top level)

sum    = 0.0
errest = 0.0
xfail  = 0.0
nfun = 0
i = 1
imax = im

x(1)   = a
h(1)   = b-a
tol(1) = eps    
il(1)  = 1

! *** stage 2 ***
! main part: adaptive integration

do while (i > 0)
 
! *** stage 2a ***
! calculate integral on [x(i),(x(i)+h(i))]
  select case (key)
    case (1)
      call G7K15(f,x(i),x(i)+h(i),gauss,kronrod)
      nfun = nfun + 15	
	case (2)
      call G10K21(f,x(i),x(i)+h(i),gauss,kronrod)
      nfun = nfun + 21    
  end select 

! *** stage 2b ***
! save data at this level
  x0 = x(i)
  step = h(i)
  err  = tol(i)
  deep = il(i)

! *** stage 2c ***
! the current level has been computed 
  i=i-1

! *** stage 2d ***
! local condition for convergence
   if((200.0*abs(kronrod-gauss))**1.5 <=err) then
    sum = sum + kronrod
    errest = errest + abs(kronrod-gauss)
    else

! *** stage 2e *** 
!   check if the code can continue to subdivide intervals
      if( deep >= imax ) then 
      xfail = xfail + step
      else
! *** stage 2f ***
!	  make smaller intervals (i.e. h=h/2)
!     data for the right subinterval
        i = i+1
        h(i)  = step/2.0
		tol(i)= err/2.0
        x(i)  = x0 + h(i)
        il(i) = deep + 1
!     data for the left subinterval
        i = i+1
        x(i)  = x0
        h(i)  = h(i-1)
        tol(i)= tol(i-1)
        il(i) = il(i-1)
	end if
  end if
end do

xfail = xfail/(b-a)
!write(*,100) xfail
!100 format(/,' xfail = ',e12.6)

return
end subroutine GKinteg



subroutine quanc8(fun,a,b,abserr,relerr,result,errest,nofun,flag)
!
      double precision fun, a, b, abserr, relerr, result, errest, flag
      integer 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.
!
double precision w0,w1,w2,w3,w4,area,x0,f0,stone,step,cor11,temp
double precision qprev,qnow,qdiff,qleft,esterr,tolerr
double precision qright(31),f(16),x(16),fsave(8,30),xsave(8,30)
double precision dabs,dmax1
integer levmin,levmax,levout,nomax,nofin,lev,nim,i,j
!
!   ***   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.0d0 / 14175.0d0
      w1 =  23552.0d0 / 14175.0d0
      w2 =  -3712.0d0 / 14175.0d0
      w3 =  41984.0d0 / 14175.0d0
      w4 = -18160.0d0 / 14175.0d0
!
!   initialize running sums to zero.
!
      flag = 0.0d0
      result = 0.0d0
      cor11  = 0.0d0
      errest = 0.0d0
      area   = 0.0d0
      nofun = 0
      if (a .eq. b) return
!
!   ***   stage 2 ***   initialization for first interval
!
      lev = 0
      nim = 1
      x0 = a
      x(16) = b
      qprev  = 0.0d0
      f0 = fun(x0)
      stone = (b - a) / 16.0d0
      x(8)  =  (x0  + x(16)) / 2.0d0
      x(4)  =  (x0  + x(8))  / 2.0d0
      x(12) =  (x(8)  + x(16)) / 2.0d0
      x(2)  =  (x0  + x(4))  / 2.0d0
      x(6)  =  (x(4)  + x(8))  / 2.0d0
      x(10) =  (x(8)  + x(12)) / 2.0d0
      x(14) =  (x(12) + x(16)) / 2.0d0
      do 25 j = 2, 16, 2
         f(j) = fun(x(j))
   25 continue
      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.
!
   30 x(1) = (x0 + x(2)) / 2.0d0
      f(1) = fun(x(1))
      do 35 j = 3, 15, 2
         x(j) = (x(j-1) + x(j+1)) / 2.0d0
         f(j) = fun(x(j))
   35 continue
      nofun = nofun + 8
      step = (x(16) - x0) / 16.0d0
      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 = dabs(qdiff) / 1023.0d0
      tolerr = dmax1(abserr,relerr*dabs(area)) * (step/stone)
      if (lev .lt. levmin) go to 50
      if (lev .ge. levmax) go to 62
      if (nofun .gt. nofin) go to 60
      if (esterr .le. tolerr) go to 70
!
!   ***   stage 5   ***   no convergence
!   locate next interval.
!
   50 nim = 2*nim
      lev = lev+1
!
!   store right hand elements for future use.
!
      do 52 i = 1, 8
         fsave(i,lev) = f(i+8)
         xsave(i,lev) = x(i+8)
   52 continue
!
!   assemble left hand elements for immediate use.
!
      qprev = qleft
      do 55 i = 1, 8
         j = -i
         f(2*j+18) = f(j+9)
         x(2*j+18) = x(j+9)
   55 continue
      go to 30
!
!   ***   stage 6   ***   trouble section
!   number of function values is about to exceed limit.
!
   60 nofin = 2*nofin
      levmax = levout
      flag = flag + (b - x0) / (b - a)
      go to 70
!
!   current level is levmax.
!
   62 flag = flag + 1.0d0
!
!   ***   stage 7   ***   interval converged
!   add contributions into running sums.
!
   70 result = result + qnow
      errest = errest + esterr
      cor11  = cor11  + qdiff / 1023.0d0
!
!   locate next interval.
!
   72 if (nim .eq. 2*(nim/2)) go to 75
      nim = nim/2
      lev = lev-1
      go to 72
   75 nim = nim + 1
      if (lev .le. 0) go to 80
!
!   assemble elements required for the next interval.
!
      qprev = qright(lev)
      x0 = x(16)
      f0 = f(16)
      do 78 i = 1, 8
         f(2*i) = fsave(i,lev)
         x(2*i) = xsave(i,lev)
   78 continue
      go to 30
!
!   ***   stage 8   ***   finalize and return
!
   80 result = result + cor11
!
!   make sure errest not less than roundoff level.
!
      if (errest .eq. 0.0d0) return
   82 temp = dabs(result) + errest
      if (temp .ne. dabs(result)) return
      errest = 2.0d0*errest
      go to 82
      end