module integrate
    use constants, only: dp
    use, intrinsic :: ieee_arithmetic, only: ieee_is_nan

    implicit none
    private

    integer, parameter :: quadkind = 8

    abstract interface
        function integrand_i(x)
            use constants, only: dp
            real(dp), intent(in) :: x
            real(dp) :: integrand_i
        end function integrand_i
    end interface

    procedure(integrand_i), private, pointer, save :: integrand => null()

    public :: integrand_i, integrate_1d, sum_trapez_1d

contains

    subroutine integrate_1d(f, a, b, result)
        use quadpack_double, only: qag => dqag
        procedure(integrand_i) :: f
        real(dp), intent(in) :: a, b
        real(dp), intent(out) :: result

        real(quadkind), parameter :: abs_error_tol_quadkind = 0.0_dp
        real(quadkind), parameter :: rel_error_tol_quadkind = 1.0e-6
        integer, parameter :: order_key = 6
        integer, parameter :: limit = 500
        integer, parameter :: lenw = limit*4
        real(quadkind) :: a_quadkind, b_quadkind, result_quadkind
        real(quadkind) :: abs_error, error_limit
        integer :: error_flag

        logical :: error_occurred

        integer :: neval_dummy, last
        integer :: iwork(limit)
        real(quadkind) :: work(lenw)

        if (ieee_is_nan(a) .or. ieee_is_nan(b)) then
            print *, "Integration limits must not be NaN!"
            print *, "a =", a, "b =", b
            error stop
        end if

        ! quadpack operates in and requires real(8)
        a_quadkind = convert_to_quadkind(a)
        b_quadkind = convert_to_quadkind(b)

        integrand => f
        ! integration by globally adaptive interval subdivision (quadpack import)
        call qag(quadkind_integrand, &
                 a_quadkind, &
                 b_quadkind, &
                 abs_error_tol_quadkind, &
                 rel_error_tol_quadkind, &
                 order_key, &
                 result_quadkind, &
                 abs_error, &
                 neval_dummy, &
                 error_flag, &
                 limit, &
                 lenw, &
                 last, &
                 iwork, &
                 work)
        integrand => null()

        error_limit = abs(result_quadkind)*rel_error_tol_quadkind &
                      + abs_error_tol_quadkind
        error_occurred = .false.

        if (abs_error > error_limit) then
            print *, "Integration warning: absolute error =", abs_error
            print *, "bigger than required ", error_limit
            print *, "relative error", abs_error/abs(result_quadkind)
            error_occurred = .true.
        end if

        if (error_flag /= 0) then
            print *, "Integration warning: error =", error_flag
            error_occurred = .true.
        end if

        result = convert_to_dp(result_quadkind)

        if (ieee_is_nan(result)) then
            print *, "Integration result is NaN!"
            error_occurred = .true.
        end if

        if (error_occurred) error stop

    end subroutine integrate_1d

    !> The integrand for quadpack needs to be of the below signature.
    !> Note for developer: due to quadpack interface, functions with parameters
    !> $$
    !> f(x, p_1, p_2, ...)
    !> $$
    !> need to be wrapped
    !>$$
    !>f(x, p_1, p_2, ...) \to f(x)
    !>$$
    !> which at some point in the call chain requires access
    !> to outer scopes where parameters are defined. If one wants to avoid
    !> executable stacks, one has to resort to module variables for that.
    !> Type bound procedures which carry the parameters with them fail latest
    !> at this boundary here. The call chain can not be made thread safe.
    !> We opt therefore for the most straight forward solution i.e. global state.
    function quadkind_integrand(x_quadkind)
        real(quadkind), intent(in) :: x_quadkind
        real(quadkind) :: quadkind_integrand

        real(dp) :: x_dp

        x_dp = real(x_quadkind, kind=dp)
        quadkind_integrand = real(integrand(x_dp), kind=quadkind)
    end function quadkind_integrand

    function convert_to_quadkind(val_dp) result(val_quadkind)
        real(dp), intent(in) :: val_dp
        real(quadkind) :: val_quadkind

        val_quadkind = real(val_dp, kind=quadkind)
    end function convert_to_quadkind

    function convert_to_dp(val_quadkind) result(val_dp)
        real(quadkind), intent(in) :: val_quadkind
        real(dp) :: val_dp

        val_dp = real(val_quadkind, kind=dp)
    end function convert_to_dp

    function sum_trapez_1d(x, y)
        real(dp), dimension(:), intent(in) :: x, y
        real(dp) :: sum_trapez_1d

        integer :: n

        n = size(x)
        if (.not. n == size(y)) then
            print *, "Error in sum_trapez_1d: x and y need to have same length!"
            error stop
        end if
        if (n < 2) then
            print *, "Error in sum_trapez_1d: len(x) must be at least 2!"
            error stop
        end if

        sum_trapez_1d = 0.5_dp*sum((x(2:n) - x(1:n - 1))*(y(2:n) + y(1:n - 1)))
    end function sum_trapez_1d

end module integrate

module integrate_substituted
    use integrate, only: integrand_i, integrate_1d
    use constants, only: dp
    use, intrinsic :: ieee_arithmetic, only: ieee_value, ieee_quiet_nan
    implicit none

    private
    procedure(integrand_i), pointer, save :: integrand => null()
    real(dp), save :: left_limit, right_limit

    public :: integrate_1d_substituted

contains

    subroutine integrate_1d_substituted(f, a, b, result)
        procedure(integrand_i) :: f
        real(dp), intent(in) :: a, b
        real(dp), intent(out) :: result

        real(dp) :: sub_a, sub_b
        real(dp) :: left_result, right_result

        if (b < a) then
            print *, "Error in integrate_1d_substituted: b < a"
            error stop
        end if

        ! integration by globally adaptive interval subdivision (quadpack import)
        ! qags struggels with functions of form sqrt(x-a) and sqrt(b-x)
        ! Terefore we split the integral in the middle and substitute
        ! - left integral: t**2 = x - a
        ! - right integral: t**2 = b - x
        ! resulting in new limits
        ! - left integral: 0 to sqrt((b - a)*0.5_dp)
        ! - right integral: 0 to sqrt((b - a)*0.5_dp) (after absorbing [-] from trafo)
        ! the substitution handled with module states
        left_limit = a
        right_limit = b
        integrand => f
        sub_a = 0.0_dp
        sub_b = sqrt((b - a)*0.5_dp)
        call integrate_1d(left_substituted_integrand, sub_a, sub_b, left_result)
        call integrate_1d(right_substituted_integrand, sub_a, sub_b, right_result)
        left_limit = ieee_value(1.0_dp, ieee_quiet_nan)
        right_limit = ieee_value(1.0_dp, ieee_quiet_nan)
        integrand => null()

        result = left_result + right_result

    end subroutine integrate_1d_substituted

    function left_substituted_integrand(t)
        real(dp), intent(in) :: t
        real(dp) :: left_substituted_integrand

        real(dp) :: x

        x = t**2.0_dp + left_limit
        left_substituted_integrand = integrand(x)*2.0_dp*t
    end function left_substituted_integrand

    function right_substituted_integrand(t)
        real(dp), intent(in) :: t
        real(dp) :: right_substituted_integrand

        real(dp) :: x

        x = right_limit - t**2.0_dp
        right_substituted_integrand = integrand(x)*2.0_dp*t
    end function right_substituted_integrand

end module integrate_substituted
