61 lines
2.1 KiB
Julia
61 lines
2.1 KiB
Julia
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
|