Beautify KG solving routine (back to old structure)

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
+"Green's function for Klein-Gordon equation in natural units"
-method with cubic-spline interpolation of the source and adaptive Simpson's integration for
+greensFunctionKG(m, r, rp) = -1 / (m * (r + r_reg) * (rp + r_reg)) * sinh(m * min(r, rp)) * exp(-m * max(r, rp))
-each grid point (matching the approach in Hartree.f). Returns an array of field values in MeV.
+
+"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
+    r_grid = range(0, s.r_max, length=s.divs+1)
-    dr = Δr(s)
-    r_max = s.r_max
-    m_nat = m / ħc    # mass in fm⁻¹
-    r_grid = range(0, r_max, length=N+1)
 
     # Cubic-spline the source for evaluation at arbitrary quadrature points
     src_spline = cubic_spline_interpolation(r_grid, source)
 
-    result = zeros(N + 1)
+    if m > 1E-10    # Massive field (Klein-Gordon equation)
-    ε = 1e-10
+        greensFunction = (r, rp) -> greensFunctionKG(m / ħc, r, rp)
-
-    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,
+        greensFunction = (r, rp) -> greensFunctionP(r, rp)
-        # 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)
+    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