dqc25c Subroutine

    subroutine dqc25c(f, a, b, c, Result, Abserr, Krul, Neval)
        implicit none

        procedure(func) :: f !! function subprogram defining the integrand function `f(x)`.
        real(wp), intent(in) :: a !! left end point of the integration interval
        real(wp), intent(in) :: b !! right end point of the integration interval, `b>a`
        real(wp), intent(in) :: c !! parameter in the weight function
        real(wp), intent(out) :: Result !! approximation to the integral.
                                        !! `result` is computed by using a generalized
                                        !! clenshaw-curtis method if `c` lies within ten percent
                                        !! of the integration interval. in the other case the
                                        !! 15-point kronrod rule obtained by optimal addition
                                        !! of abscissae to the 7-point gauss rule, is applied.
        real(wp), intent(out) :: Abserr !! estimate of the modulus of the absolute error,
                                        !! which should equal or exceed `abs(i-result)`
        integer, intent(inout) :: Krul !! key which is decreased by 1 if the 15-point
                                       !! gauss-kronrod scheme has been used
        integer, intent(out) :: Neval !! number of integrand evaluations

        real(wp) :: ak22, amom0, amom1, amom2, cc, &
                    p2, p3, p4, resabs, resasc, u
        integer :: i, isym, k
        real(wp) :: fval(25) !! value of the function `f` at the points
                             !! `cos(k*pi/24)`, `k = 0, ..., 24`
        real(wp) :: cheb12(13) !! chebyshev series expansion coefficients,
                               !! for the function `f`, of degree 12
        real(wp) :: cheb24(25) !! chebyshev series expansion coefficients,
                               !! for the function `f`, of degree 24
        real(wp) :: res12 !! approximation to the integral corresponding
                          !! to the use of cheb12
        real(wp) :: res24 !! approximation to the integral corresponding
                          !! to the use of cheb24
        real(wp) :: hlgth !! half-length of the interval
        real(wp) :: centr !! mid point of the interval

        integer,parameter :: kp = 0 !! unused variable for [[dqwgtc]]

        real(wp), dimension(11), parameter :: x = [(cos(k*pi/24.0_wp), k=1, 11)]
            !! the vector x contains the values `cos(k*pi/24)`,
            !! `k = 1, ..., 11`, to be used for the chebyshev series
            !! expansion of `f`

        ! check the position of c.

        cc = (2.0_wp*c - b - a)/(b - a)
        if (abs(cc) < 1.1_wp) then

            ! use the generalized clenshaw-curtis method.

            hlgth = 0.5_wp*(b - a)
            centr = 0.5_wp*(b + a)
            Neval = 25
            fval(1) = 0.5_wp*f(hlgth + centr)
            fval(13) = f(centr)
            fval(25) = 0.5_wp*f(centr - hlgth)
            do i = 2, 12
                u = hlgth*x(i - 1)
                isym = 26 - i
                fval(i) = f(u + centr)
                fval(isym) = f(centr - u)
            end do

            ! compute the chebyshev series expansion.

            call dqcheb(x, fval, cheb12, cheb24)

            ! the modified chebyshev moments are computed by forward
            ! recursion, using amom0 and amom1 as starting values.

            amom0 = log(abs((1.0_wp - cc)/(1.0_wp + cc)))
            amom1 = 2.0_wp + cc*amom0
            res12 = cheb12(1)*amom0 + cheb12(2)*amom1
            res24 = cheb24(1)*amom0 + cheb24(2)*amom1
            do k = 3, 13
                amom2 = 2.0_wp*cc*amom1 - amom0
                ak22 = (k - 2)*(k - 2)
                if ((k/2)*2 == k) amom2 = amom2 - 4.0_wp/(ak22 - 1.0_wp)
                res12 = res12 + cheb12(k)*amom2
                res24 = res24 + cheb24(k)*amom2
                amom0 = amom1
                amom1 = amom2
            end do
            do k = 14, 25
                amom2 = 2.0_wp*cc*amom1 - amom0
                ak22 = (k - 2)*(k - 2)
                if ((k/2)*2 == k) amom2 = amom2 - 4.0_wp/(ak22 - 1.0_wp)
                res24 = res24 + cheb24(k)*amom2
                amom0 = amom1
                amom1 = amom2
            end do
            Result = res24
            Abserr = abs(res24 - res12)
        else

            ! apply the 15-point gauss-kronrod scheme.

            ! dqwgtc - external function subprogram defining the weight function

            Krul = Krul - 1
            call dqk15w(f, dqwgtc, c, p2, p3, p4, kp, a, b, Result, Abserr, resabs, &
                        resasc)
            Neval = 15
            if (resasc == Abserr) Krul = Krul + 1
        end if

    end subroutine dqc25c