From bee055dcf1a2b62aeb9e70e17a9d0546aa0c211c Mon Sep 17 00:00:00 2001
From: Nuwan Yapa <yapanuwan@gmail.com>
Date: Thu, 2 Apr 2026 21:49:41 -0400
Subject: [PATCH] Beautify KG solving routine (back to old structure)

---
 mesons.jl | 57 ++++++++++++++++++-------------------------------------
 1 file changed, 18 insertions(+), 39 deletions(-)

diff --git a/mesons.jl b/mesons.jl
index e997ad7..d7410da 100644
--- a/mesons.jl
+++ b/mesons.jl
@@ -1,56 +1,35 @@
 include("common.jl")
 include("system.jl")
 
-"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.
+"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) 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)
-    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)
+    r_grid = range(0, s.r_max, length=s.divs+1)
 
     # Cubic-spline the source for evaluation at arbitrary quadrature points
     src_spline = cubic_spline_interpolation(r_grid, source)
 
-    result = zeros(N + 1)
-    ε = 1e-10
-
-    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
+    if m > 1E-10    # Massive field (Klein-Gordon equation)
+        greensFunction = (r, rp) -> greensFunctionKG(m / ħc, r, rp)
     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, ε)
+        greensFunction = (r, rp) -> greensFunctionP(r, rp)
+    end
 
-            I1 = x > ε ? adaptive_simps(0.0, x) do r
-                r^2 * src_spline(r)
-            end : 0.0
+    result = zeros(s.divs + 1)
 
-            I2 = adaptive_simps(x, r_max) do r
-                r * src_spline(r)
-            end
-
-            result[i+1] = -ħc * (I1 / x + I2)
+    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