d_simps_ab_func Function

function d_simps_ab_func(f,a,b,N) result(int)
  interface
     function f(x)
       real(8) :: x
       real(8) :: f
     end function f
  end interface
  real(8),optional                 :: a,b
  integer,optional                 :: N
  real(8)                          :: dh,a_,b_
  integer                          :: L,M,i
  real(8),dimension(:),allocatable :: xx,wt,dx
  real(8)                          :: int,int1,int2,int3
  L  = 100;if(present(N))L =N
  a_ = 0d0;if(present(a))a_=a
  b_ = 1d0;if(present(b))b_=b
  !
  int=0.d0
  allocate(xx(L),wt(L))
  xx = linspace(a_,b_,L,mesh=dh)
  call get_quadrature_weights(wt)
  do i=1,L
     int = int + f(xx(i))*wt(i)
  enddo
  int = int*dh
  deallocate(xx,wt)
end function d_simps_ab_func