//#####################################################################
// Copyright 2020, Craig Schroeder.
// This file is part of PhysBAM whose distribution is governed by the license contained in the accompanying file PHYSBAM_COPYRIGHT.txt.
//#####################################################################
// Class Integrate_Function
//#####################################################################
#ifndef __Integrate_Function__
#define __Integrate_Function__

#include <Core/Utilities/TYPE_UTILITIES.h>
#include <limits>
namespace PhysBAM{

template<class T,class F>
auto Integrate_Function_Quadrature(F f,T a,T b)
{
    T dx=(T)0.112701665379258311*(b-a);
    auto i0=f((T).5*(a+b))*(T)0.444444444444444444;
    auto i1=(f(a+dx)+f(b-dx))*(T)0.277777777777777778;
    return (i0+i1)*(b-a);
}

template<class T,class U,class F>
U Integrate_Function_Helper(F f,U fab,T a,T b,T tol,int max_depth)
{
    T c = .5*(a+b);
    U fac=Integrate_Function_Quadrature(f,a,c);
    U fcb=Integrate_Function_Quadrature(f,c,b);
    U facb=fac+fcb;
    if(max_depth==0 || abs(fab-facb)<=tol) return facb;
    U rac=Integrate_Function_Helper(f,fac,a,c,tol/2,max_depth-1);
    U rcb=Integrate_Function_Helper(f,fcb,c,b,tol/2,max_depth-1);
    return rac+rcb;
}

// compute int(f(x),x=a..b);
template<class T,class F>
auto Integrate_Function(F f,T a,T b,T tol=std::numeric_limits<T>::epsilon()*100,
    enable_if_t<IS_FLOATING_POINT<T>::value,int> max_depth=32)
{
    auto fab=Integrate_Function_Quadrature(f,a,b);
    return Integrate_Function_Helper(f,fab,a,b,tol,max_depth);
}
//#####################################################################
}
#endif
