common.jl

@@ -15,46 +15,3 @@ 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,32 +7,18 @@ 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) 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
+"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/c²,
+    m is the mass of the meson in MeV/c2."
-    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)
 
-    # Cubic-spline the source for evaluation at arbitrary quadrature points
+    @assert s.divs % 2 == 0 "Number of mesh divisions must be even for Simpson's rule"
-    src_spline = cubic_spline_interpolation(r_grid, source)
+    simpsons_weights = (Δr(s)/3) .* [1; repeat([2,4], s.divs ÷ 2)[2:end]; 1]
+    int_measure = ħc .* simpsons_weights .* rs(s) .^ 2
 
-    if m > 1E-10    # Massive field (Klein-Gordon equation)
+    greensFunction = m == 0 ? greensFunctionP : (r, rp) -> greensFunctionKG(m / ħc, r, rp)
-        greensFunction = (r, rp) -> greensFunctionKG(m / ħc, r, rp)
+    greenMat = greensFunction.(rs(s), transpose(rs(s)))
-    else    # Massless field (Poisson/Coulomb equation)
-        greensFunction = (r, rp) -> greensFunctionP(r, rp)
-    end
 
-    result = zeros(s.divs + 1)
+    return greenMat * (int_measure .* source)
-
-    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