! test_aquad3.f90 adaptive quadrature test 

function f(x) result(y)
  implicit none
  double precision, intent(in) :: x
  double precision :: y
  y = 1.0d0/x
end function f

function f1(x) result(y)
  implicit none
  double precision, intent(in) :: x
  double precision :: y
  y = 1.0d0/((x-0.3d0)*(x-0.3d0)+0.01d0) +  &
      1.0d0/((x-0.9d0)*(x-0.9d0)+0.04d0) -6.0d0
end function f1

program test_aquad3
  implicit none
  double precision :: xmin
  double precision :: xmax
  double precision :: eps  = 0.001d0
  double precision :: area, err, exact

  interface
    function f(x) result(y)
      double precision, intent(in) :: x
      double precision :: y
    end function f
    function f1(x) result(y)
      double precision, intent(in) :: x
      double precision :: y
    end function f1
    function aquad3(f, xmin, xmax, eps) result(integral)
      double precision, intent(in) :: xmin, xmax, eps
      double precision :: integral
      interface
        function f(x) result(y)
          double precision, intent(in) :: x
          double precision :: y
        end function f
      end interface
    end function aquad3
  end interface

  print *, 'test_aquad3.f90 testing aquad3.c 1/x  eps=', eps

  xmin = 0.1d0
  xmax = 2.0d0
  area = aquad3(f, xmin, xmax, eps)
  exact = log(xmax)-log(xmin)
  err = area-exact
  print '("xmin=",f7.5,", xmax=",f3.0,", area=",f12.8,", exact=",f12.8,", err=",e12.4)', &
           xmin, xmax, area, exact, err   

  xmin = 0.01d0
  area = aquad3(f, xmin, xmax, eps)
  exact = log(xmax)-log(xmin)
  err = area-exact
  print '("xmin=",f7.5,", xmax=",f3.0,", area=",f12.8,", exact=",f12.8,", err=",e12.4)', &
           xmin, xmax, area, exact, err   

  xmin = 0.001d0
  area = aquad3(f, xmin, xmax, eps)
  exact = log(xmax)-log(xmin)
  err = area-exact
  print '("xmin=",f7.5,", xmax=",f3.0,", area=",f12.8,", exact=",f12.8,", err=",e12.4)', &
           xmin, xmax, area, exact, err   

  xmin = 0.0001d0
  area = aquad3(f, xmin, xmax, eps)
  exact = log(xmax)-log(xmin)
  err = area-exact
  print '("xmin=",f7.5,", xmax=",f3.0,", area=",f12.8,", exact=",f12.8,", err=",e12.4)', &
           xmin, xmax, area, exact, err   

  xmin = 0.00001d0
  area = aquad3(f, xmin, xmax, eps)
  exact = log(xmax)-log(xmin)
  err = area-exact
  print '("xmin=",f7.5,", xmax=",f3.0,", area=",f12.8,", exact=",f12.8,", err=",e12.4)', &
           xmin, xmax, area, exact, err   

  print *, '1.0/((x-0.3)*(x-0.3)+0.01) + 1.0/((x-0.9)*(x-0.9)+0.04) -6'
  xmin = 0.0d0
  xmax = 1.0d0
  eps = 0.001
  area = aquad3(f1, xmin, xmax, eps)
  exact = 29.85832540
  err = area - exact  
  print '("xmin=",f7.5,", xmax=",f3.0,", area=",f12.8,", exact=",f12.8,", err=",e12.4)', &
           xmin, xmax, area, exact, err   

  print *, 'test_aquad3.f90 finished'
end program test_aquad3