const ħc = 197.33 # MeVfm
const r_reg = 1E-8 # fm # regulator for R

"Integrate a uniformly discretized function f using Simpson's rule where h is the step size and coefficient is an optional scaling factor."
function simpsons_integrate(f::AbstractVector{Float64}, h::Float64; coefficient::Float64 = 1.0)
    @assert length(f) % 2 == 1 "Number of mesh divisions must be even for Simpson's rule"
    s = sum(enumerate(f)) do (i, fi)
        if i == 1 || i == length(f)
            return fi
        elseif i % 2 == 0
            return 4fi
        else
            return 2fi
        end
    end
    return (h/3) * coefficient * s
end

"Adaptive Simpson's quadrature for a function f on [a,b].
    Recursively subdivides into three sub-panels until |S_new - S_old| < eps×area,
    matching the FORTRAN 'simps' routine by K. Wehrberger. Returns the integral estimate."
function adaptive_simps(f, a::Float64, b::Float64; tol::Float64=1e-5, max_depth::Int=20)
    a ≥ b && return 0.0
    fa = f(a)
    fm = 4.0 * f((a + b) / 2)
    fb = f(b)
    area = 1.0
    est = 1.0
    return _simps_recurse(f, a, b - a, fa, fm, fb, area, est, tol, 0, max_depth)
end

"Internal recursive worker for adaptive_simps (3-panel subdivision matching FORTRAN simps)."
function _simps_recurse(f, a, da, fa, fm, fb, area, est, eps, depth, max_depth)
    depth ≥ max_depth && return est  # bail out at max depth

    dx = da / 3
    x1 = a + dx
    x2 = x1 + dx
    f1 = 4.0 * f(a + 0.5 * dx)
    f2 = f(x1)
    f3 = f(x2)
    f4 = 4.0 * f(a + 2.5 * dx)
    dx6 = dx / 6

    est1 = (fa + f1 + f2) * dx6
    est2 = (f2 + fm + f3) * dx6
    est3 = (f3 + f4 + fb) * dx6
    area = area - abs(est) + abs(est1) + abs(est2) + abs(est3)
    s = est1 + est2 + est3

    if abs(est - s) ≤ eps * area && est != 1.0
        return s
    end

    eps_child = eps / 1.7
    s1 = _simps_recurse(f, a,  dx, fa, f1, f2, area, est1, eps_child, depth + 1, max_depth)
    s2 = _simps_recurse(f, x1, dx, f2, fm, f3, area, est2, eps_child, depth + 1, max_depth)
    s3 = _simps_recurse(f, x2, dx, f3, f4, fb, area, est3, eps_child, depth + 1, max_depth)
    return s1 + s2 + s3
end