Implemented Green's function (not working)

dirac.jl

@@ -27,7 +27,7 @@ end
     divs is the number of mesh divisions so solution would be returned as a 2×(1+divs) matrix,
     shooting method divides the interval into two partitions at r_max/2, ensuring convergence at both r=0 and r=r_max,
     the other parameters are the same from dirac!(...)."
-function solveWf(κ, p, E, Φ0, W0, B0, A0, r_max, divs; shooting=true, normalize=true)
+function solveDiracWf(κ, p, E, Φ0, W0, B0, A0, r_max, divs; shooting=true, normalize=true)
     saveat = r_max / divs
 
     if shooting
@@ -124,7 +124,7 @@ function calculateNucleonDensity(N, p, Φ0, W0, B0, A0, r_max, divs, E_min=0, E_
     ρr2 = zeros(2, divs + 1) # ρ×r² values
     
     for (κ, E, occ) in zip(κs, Es, occs)
-        wf = solveWf(κ, p, E, Φ0, W0, B0, A0, r_max, divs; shooting=true, normalize=true)
+        wf = solveDiracWf(κ, p, E, Φ0, W0, B0, A0, r_max, divs; shooting=true, normalize=true)
         wf2 = wf .* wf
         ρr2 += (occ / (4 * pi)) * wf2 # 2j+1 factor is accounted in the occupancy number
     end

mesons.jl

@@ -17,36 +17,47 @@ const λ = -0.003551486718 # LambdaSS
 const ζ = 0.023499504053 # LambdaVV
 const Λv = 0.043376933644 # LambdaVR
 
-const r_reg = 1E-6 # regulator for the centrifugal term in fm
+const r_reg = 1E-9 # regulator for Green's functions in fm
 
-"Coupled radial meson equations that returns ddu=[ddΦ0, ddW0, ddB0, ddA0] in-place where
+"Green's function for Klein-Gordon equation"
-    du=[dΦ0, dW0, dB0, dA0] are the first derivatives of meson fields evaluated at r,
+greensFunctionKG(m, r, rp) = -(1 / m) * rp / (r + r_reg) * sinh(m * min(r, rp)) * exp(-m * max(r, rp))
-    u=[Φ0, W0, B0, A0] are the meson fields evaluated at r,
+
-    κ is the generalized angular momentum,
+"Green's function for Poisson's equation"
-    ρ_sp, ρ_vp are the scalar and vector proton densities as funtions of r in fm,
+greensFunctionP(r, rp) = -rp^2 / (max(r, rp) + r_reg)
-    ρ_sn, ρ_vn are the scalar and vector neutron densities as funtions of r in fm,
+
-    Φ0, W0, B0, A0 are the mean-field potentials (couplings included) in MeV as functions of r in fm,
+"Solve the Klein-Gordon equation (or Poisson's equation when m=0) for a given a source function (as an array) where
+    the other parameters are the same from mesons!(...)."
+function solveKG(m, source, r_max)
+    greensFunction = m == 0 ? greensFunctionP : (r, rp) -> greensFunctionKG(m, r, rp)
+    mesh = range(0, r_max; length=length(source))
+    greenMat = greensFunction.(transpose(mesh), mesh)
+    return greenMat * source
+end
+
+"Iteratively solve meson equations and return the wave functions u(r)=[Φ0(r), W0(r), B0(r), A0(r)] where
+    divs is the number of mesh divisions so the arrays are of length (1+divs),
+    ρ_sp, ρ_vp are the scalar and vector proton densities as arrays,
+    ρ_sn, ρ_vn are the scalar and vector neutron densities as arrays,
+    Φ0, W0, B0, A0 are the mean-field potentials (couplings included) in MeV as as arrays,
     r is the radius in fm.
 Reference: P. Giuliani, K. Godbey, E. Bonilla, F. Viens, and J. Piekarewicz, Frontiers in Physics 10, (2023)"
-function mesons!(ddu, du, u, (ρ_sp, ρ_vp, ρ_sn, ρ_vn), r)
+function solveMesonWfs(ρ_sp, ρ_vp, ρ_sn, ρ_vn, r_max, divs, iterations=3)
-    (dΦ0, dW0, dB0, dA0) = du
+    (Φ0, W0, B0, A0) = (zeros(1 + divs) for _ in 1:4)
-    (Φ0, W0, B0, A0) = u
+    (src_Φ0, src_W0, src_B0, src_A0) = (zeros(1 + divs) for _ in 1:4)
 
-    ddΦ0 = -(2/(r + r_reg)) * dΦ0 + m_s^2 * Φ0 / ħc + g2_s * ((κ/2) * Φ0^2 + (λ/6) * Φ0^3 - (ρ_sp(r) + ρ_sn(r))) / ħc
+    for _ in iterations
-    ddW0 = -(2/(r + r_reg)) * dW0 + m_ω^2 * W0 / ħc + g2_v * ((ζ/6) * W0^3 + 2 * Λv * B0^2 * W0 - (ρ_vp(r) + ρ_vn(r))) / ħc
+        @. src_Φ0 = g2_s * ((κ/2) * Φ0^2 + (λ/6) * Φ0^3 - (ρ_sp + ρ_sn)) / ħc
-    ddB0 = -(2/(r + r_reg)) * dB0 + m_ρ^2 * B0 / ħc + g2_ρ * (2 * Λv * W0^2 * B0 - (ρ_vp(r) - ρ_vn(r)) / 2) / ħc
+        Φ0 .= solveKG(m_s / sqrt(ħc), src_Φ0, r_max)
-    ddA0 = -(2/(r + r_reg)) * dA0 - ρ_vp(r) / ħc # e goes here?
 
-    ddu .= [ddΦ0, ddW0, ddB0, ddA0]
+        @. src_W0 = g2_v * ((ζ/6) * W0^3 + 2 * Λv * B0^2 * W0 - (ρ_vp + ρ_vn)) / ħc
+        W0 .= solveKG(m_ω / sqrt(ħc), src_W0, r_max)
+
+        @. src_B0 = g2_ρ * (2 * Λv * W0^2 * B0 - (ρ_vp - ρ_vn) / 2) / ħc
+        W0 .= solveKG(m_ρ / sqrt(ħc), src_B0, r_max)
+
+        @. src_A0 = ρ_vp / ħc
+        A0 .= solveKG(0, src_A0, r_max) # this doesn't need iterations
     end
 
-"Solve meson equations and return the wave functions u(r)=[Φ0(r), W0(r), B0(r), A0(r)] where
+    return (Φ0, W0, B0, A0)
-    divs is the number of mesh divisions so solution would be returned as a 4×(1+divs) matrix,
-    the other parameters are the same from mesons!(...)."
-function solveWfs(ρ_sp, ρ_vp, ρ_sn, ρ_vn, r_max, divs)
-    prob = SecondOrderODEProblem(mesons!, BigFloat.([0, 0, 0, 0]), BigFloat.([0, 0, 0, 0]), (r_max, 0))
-    sol = solve(prob, Feagin14(), p=(ρ_sp, ρ_vp, ρ_sn, ρ_vn), saveat=r_max/divs)
-    wfs = hcat(sol.u...)
-
-    return wfs
 end

test/Pb208_dens.jl

@@ -25,7 +25,7 @@ groundE = findEs(κ, p, S_interp, V_interp, R_interp, A_interp, r_max, E_min, E_
 println("ground state E = $groundE")
 
 divs = 400
-wf = solveWf(κ, p, groundE, S_interp, V_interp, R_interp, A_interp, r_max, divs)
+wf = solveDiracWf(κ, p, groundE, S_interp, V_interp, R_interp, A_interp, r_max, divs)
 rs = range(0, r_max, length=divs+1)
 gs = wf[1, :]
 fs = wf[2, :]

test/Pb208_flds.jl

@@ -0,0 +1,19 @@
+using DelimitedFiles, Interpolations, Plots
+include("../mesons.jl")
+
+# test data generated from Hartree.f
+# format: x  Rhos(n)  Rhov(n)  Rhot(n)  Rhos(p)  Rhov(p)  Rhot(p)
+test_data = readdlm("test/Pb208DensFSUGarnet.csv")
+xs = test_data[:, 1]
+ρ_sn = test_data[:, 2]
+ρ_vn = test_data[:, 3]
+ρ_sp = test_data[:, 5]
+ρ_vp = test_data[:, 6]
+
+r_max = maximum(xs)
+divs = length(xs) - 1
+
+(Φ0, W0, B0, A0) = solveMesonWfs(ρ_sp, ρ_vp, ρ_sn, ρ_vn, r_max, divs, 3)
+
+plot(xs, hcat(Φ0, W0, B0, A0), layout=4, label=["Φ0" "W0" "B0" "A0"])
+xlabel!("r (fm)")

test/Pb208_wfs.jl

@@ -1,26 +0,0 @@
-using DelimitedFiles, Interpolations, Plots
-include("../mesons.jl")
-
-# test data generated from Hartree.f
-# format: x  Rhos(n)  Rhov(n)  Rhot(n)  Rhos(p)  Rhov(p)  Rhot(p)
-test_data = readdlm("test/Pb208DensFSUGarnet.csv")
-xs = test_data[:, 1]
-rho_sn = test_data[:, 2]
-rho_vn = test_data[:, 3]
-rho_sp = test_data[:, 5]
-rho_vp = test_data[:, 6]
-
-ρ_sn = linear_interpolation(xs, rho_sn)
-ρ_vn = linear_interpolation(xs, rho_vn)
-ρ_sp = linear_interpolation(xs, rho_sp)
-ρ_vp = linear_interpolation(xs, rho_vp)
-
-r_max = maximum(xs)
-divs = 400
-
-wfs = solveWfs(ρ_sp, ρ_vp, ρ_sn, ρ_vn, r_max, divs)
-
-rs = range(0, r_max, length=divs+1)
-
-plot(rs, transpose(wfs))
-xlabel!("r (fm)")