Program main
!------------------------------------------------------------------
! Integration of a function using adaptive procedure 
! with Gauss 8 and 16 points 
!------------------------------------------------------------------
implicit none
double precision gauss_a, a, b, f, integral, eps
external f
integer nint
double precision, parameter :: pi = 3.1415926
 a = 0.0
 b = 2.0*pi
 eps = 1.0e-7
 integral = gauss_a(f,a,b,eps)
 write(*,103) integral
103 format(' integral  = ',f12.8)
stop
End

  Function f(x)
!----------------------------------------
! Function for integration
!----------------------------------------
implicit none
double precision f, x
! f = sin(x)
! f = x*cos(10.0*x**2)/(x**2 + 1.0)
 f = x*sin(30.0*x)*cos(50.0*x)
return
end

 recursive function gauss_a(f,a,b,eps)
!==========================================================
! Integration of f(x) on [a,b]
! Method: Gauss quadratures with adaptive integration  
! for left/right intervals till  error = |I_16 - I_8| < eps
! written by: Alex Godunov (October 2009)
!----------------------------------------------------------
! IN:
! f   - Function to integrate (supplied by a user)
! a	  - Lower limit of integration
! b	  - Upper limit of integration
! eps - tolerance (should not be less than 1.0e-8)
! OUT:
! gauss_a - Result of integration
!==========================================================
implicit none
double precision gauss_a, f, a, b, eps, gauss8, gauss16
double precision s1, s2, h, ax, bx, sum
integer i
external f

if(eps <= 1.0e-8) eps = 1.0e-8

h = (b-a)/2.0
sum = 0.0

do i=1,2
   ax = a + h*dfloat(i-1)
   bx = ax + h
   s1 =  gauss8(f,ax,bx)
   s2 = gauss16(f,ax,bx)
   if(abs(s2-s1)<= eps .and. abs(s2-s1)/abs(s2+s1)<= eps) then
      sum = sum + s2
   else
      sum = sum + gauss_a(f,ax,bx,eps)
   end if
end do
gauss_a = sum
return
end function gauss_a

  Function gauss8(f,a,b)
!==========================================================
! Integration of f(x) on [a,b]
! Method: Gauss 8 points  
! written by: Alex Godunov (October 2009)
!----------------------------------------------------------
! IN:
! f   - Function to integrate (supplied by a user)
! a	  - Lower limit of integration
! b	  - Upper limit of integration
! OUT:
! gauss8 - Result of integration
!==========================================================
implicit none
integer, parameter :: n=4
double precision gauss8, f, a, b
double precision ti(n), ci(n)
data ti/0.1834346424, 0.5255324099, 0.7966664774, 0.9602898564/
data ci/0.3626837833, 0.3137066458, 0.2223810344, 0.1012285362/ 
double precision r, m, c
integer i

r = 0.0;
m = (b-a)/2.0;
c = (b+a)/2.0;

do i = 1,n 
   r = r + ci(i)*(f(m*(-1.0)*ti(i) + c) + f(m*ti(i) + c))
end do
gauss8 = r*m
return
end function gauss8

  Function gauss16(f,a,b)
!==========================================================
! Integration of f(x) on [a,b]
! Method: Gauss 16 points  
! written by: Alex Godunov (October 2009)
!----------------------------------------------------------
! IN:
! f   - Function to integrate (supplied by a user)
! a	  - Lower limit of integration
! b	  - Upper limit of integration
! OUT:
! gauss16 - Result of integration
!==========================================================
implicit none
integer, parameter :: n=8
double precision gauss16, f, a, b
double precision ti(n), ci(n)
data ti/0.0950125098, 0.2816035507, 0.4580167776, 0.6178762444, &  
        0.7554044083, 0.8656312023, 0.9445750230, 0.9894009349/ 
data ci/0.1894506104, 0.1826034150, 0.1691565193, 0.1495959888, &
        0.1246289712, 0.0951585116, 0.0622535239, 0.0271524594/ 
double precision r, m, c
integer i

r = 0.0;
m = (b-a)/2.0;
c = (b+a)/2.0;

do i = 1,n 
   r = r + ci(i)*(f(m*(-1.0)*ti(i) + c) + f(m*ti(i) + c))
end do
gauss16 = r*m
return
end function gauss16