Beautify KG solving routine (back to old structure)

common.jl

@@ -15,3 +15,46 @@ function simpsons_integrate(f::AbstractVector{Float64}, h::Float64; coefficient:
     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

mesons.jl

@@ -7,18 +7,32 @@ greensFunctionKG(m, r, rp) = -1 / (m * (r + r_reg) * (rp + r_reg)) * sinh(m * mi
 "Green's function for Poisson's equation in natural units"
 greensFunctionP(r, rp) = -1 / (max(r, rp) + r_reg)
 
-"Solve the Klein-Gordon equation (or Poisson's equation when m=0) and return an array in MeV for a source array given in fm⁻³ where
-    m is the mass of the meson in MeV/c2."
+"Solve the Klein-Gordon equation (or Poisson's equation when m=0) using the Green's function method with cubic-spline interpolation of the source and adaptive Simpson's integration for each grid point. Returns an array of field values in MeV, where
+    m is the mass of the meson in MeV/c²,
+    source is the source density array in fm⁻³."
 function solveKG(m, source, s::system)
+    r_grid = range(0, s.r_max, length=s.divs+1)
 
-    @assert s.divs % 2 == 0 "Number of mesh divisions must be even for Simpson's rule"
-    simpsons_weights = (Δr(s)/3) .* [1; repeat([2,4], s.divs ÷ 2)[2:end]; 1]
-    int_measure = ħc .* simpsons_weights .* rs(s) .^ 2
+    # Cubic-spline the source for evaluation at arbitrary quadrature points
+    src_spline = cubic_spline_interpolation(r_grid, source)
 
-    greensFunction = m == 0 ? greensFunctionP : (r, rp) -> greensFunctionKG(m / ħc, r, rp)
-    greenMat = greensFunction.(rs(s), transpose(rs(s)))
+    if m > 1E-10    # Massive field (Klein-Gordon equation)
+        greensFunction = (r, rp) -> greensFunctionKG(m / ħc, r, rp)
+    else    # Massless field (Poisson/Coulomb equation)
+        greensFunction = (r, rp) -> greensFunctionP(r, rp)
+    end
 
-    return greenMat * (int_measure .* source)
+    result = zeros(s.divs + 1)
+
+    for i in 0:s.divs
+        r = r_grid[i + 1]
+        I = adaptive_simps(0.0, s.r_max) do rp
+            rp^2 * greensFunction(r, rp) * src_spline(rp)
+        end
+        result[i+1] = ħc * I
+    end
+
+    return result
 end
 
 "Iteratively solve meson equations and return the wave functions u(r)=[Φ0(r), W0(r), B0(r), A0(r)] where