From 912c90838e829be22e23275e6987bf4a916ffac4 Mon Sep 17 00:00:00 2001
From: Nuwan Yapa <yapanuwan@gmail.com>
Date: Thu, 2 Apr 2026 21:16:28 -0400
Subject: [PATCH] Adaptive Simpson's for KG eqns (finally solves discrepancy)

---
 common.jl | 43 +++++++++++++++++++++++++++++++++++++
 mesons.jl | 63 ++++++++++++++++++++++++++++++++++++++++++-------------
 2 files changed, 92 insertions(+), 14 deletions(-)

diff --git a/common.jl b/common.jl
index c0b7503..ab281fe 100644
--- a/common.jl
+++ b/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
diff --git a/mesons.jl b/mesons.jl
index b3888e0..e997ad7 100644
--- a/mesons.jl
+++ b/mesons.jl
@@ -1,24 +1,59 @@
 include("common.jl")
 include("system.jl")
 
-"Green's function for Klein-Gordon equation in natural units"
-greensFunctionKG(m, r, rp) = -1 / (m * (r + r_reg) * (rp + r_reg)) * sinh(m * min(r, rp)) * exp(-m * max(r, rp))
-
-"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 (matching the approach in Hartree.f). Returns an array of field values in MeV.
+    m is the mass of the meson in MeV/c²,
+    source is the source density array in fm⁻³."
 function solveKG(m, source, s::system)
+    N = s.divs
+    dr = Δr(s)
+    r_max = s.r_max
+    m_nat = m / ħc    # mass in fm⁻¹
+    r_grid = range(0, r_max, length=N+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)))
+    result = zeros(N + 1)
+    ε = 1e-10
 
-    return greenMat * (int_measure .* source)
+    if m_nat > ε    # Massive field (Klein-Gordon equation)
+        for i in 0:N
+            x = max(i * dr, ε)  # regularize at r=0
+
+            # ∫₀ˣ r' × sinh(m×r') × exp(-m×x) × S(r') dr'
+            I1 = x > ε ? adaptive_simps(0.0, x) do r
+                r * sinh(m_nat * r) * exp(-m_nat * x) * src_spline(r)
+            end : 0.0
+
+            # ∫ₓ^r_max r' × sinh(m×x) × exp(-m×r') × S(r') dr'
+            I2 = adaptive_simps(x, r_max) do r
+                r * sinh(m_nat * x) * exp(-m_nat * r) * src_spline(r)
+            end
+
+            result[i+1] = -ħc * (I1 + I2) / (m_nat * x)
+        end
+    else    # Massless field (Poisson/Coulomb equation)
+        # In the m→0 limit: sinh(m r')/m → r' and exp(-m r) → 1,
+        # so field(x) = (1/x) ∫₀ˣ r'² S(r') dr' + ∫ₓ^∞ r' S(r') dr'
+        for i in 0:N
+            x = max(i * dr, ε)
+
+            I1 = x > ε ? adaptive_simps(0.0, x) do r
+                r^2 * src_spline(r)
+            end : 0.0
+
+            I2 = adaptive_simps(x, r_max) do r
+                r * src_spline(r)
+            end
+
+            result[i+1] = -ħc * (I1 / x + I2)
+        end
+    end
+
+    return result
 end
 
 "Iteratively solve meson equations and return the wave functions u(r)=[Φ0(r), W0(r), B0(r), A0(r)] where