From c9ec9c62035fb689c876649eba6e5da4f92f856c Mon Sep 17 00:00:00 2001
From: Nuwan Yapa <yapanuwan@gmail.com>
Date: Thu, 30 Jan 2025 18:36:29 -0500
Subject: [PATCH] Initial solution for mesons

---
 mesons.jl                |  7 ++++---
 test/Pb208_meson_flds.jl | 28 ++++++++++++++--------------
 2 files changed, 18 insertions(+), 17 deletions(-)

diff --git a/mesons.jl b/mesons.jl
index 0295aa1..0a7ec90 100644
--- a/mesons.jl
+++ b/mesons.jl
@@ -44,10 +44,11 @@ end
     ρ_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.
+    r is the radius in fm,
+    An initial guess initial_sol=(Φ0, W0, B0, A0) can be provided to speed up convergence (permuting!).
 Reference: P. Giuliani, K. Godbey, E. Bonilla, F. Viens, and J. Piekarewicz, Frontiers in Physics 10, (2023)"
-function solveMesonWfs(ρ_sp, ρ_vp, ρ_sn, ρ_vn, r_max, divs, iterations=3)
-    (Φ0, W0, B0, A0) = (zeros(1 + divs) for _ in 1:4)
+function solveMesonWfs(ρ_sp, ρ_vp, ρ_sn, ρ_vn, r_max, divs, iterations=3; initial_sol=(zeros(1 + divs) for _ in 1:4))
+    (Φ0, W0, B0, A0) = initial_sol
     (src_Φ0, src_W0, src_B0, src_A0) = (zeros(1 + divs) for _ in 1:4)
 
     # A0 doesn't need iterations
diff --git a/test/Pb208_meson_flds.jl b/test/Pb208_meson_flds.jl
index b8b75b4..f8b11ac 100644
--- a/test/Pb208_meson_flds.jl
+++ b/test/Pb208_meson_flds.jl
@@ -1,6 +1,17 @@
 using DelimitedFiles, Interpolations, Plots
 include("../mesons.jl")
 
+# test data generated from Hartree.f
+# format: x S(x) V(x) R(x) A(x)
+bench_data = readdlm("test/Pb208FldsFSUGarnet.csv") 
+xs_bench = bench_data[:, 1]
+Φ0_bench = bench_data[:, 2]
+W0_bench = bench_data[:, 3]
+B0_bench = bench_data[:, 4]
+A0_bench = bench_data[:, 5]
+
+plot(xs_bench, hcat(Φ0_bench, W0_bench, B0_bench, A0_bench), layout=4, label=["Φ0 benchmark" "W0 benchmark" "B0 benchmark" "A0 benchmark"])
+
 # 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")
@@ -13,18 +24,7 @@ xs = test_data[:, 1]
 r_max = maximum(xs)
 divs = length(xs) - 1
 
-(Φ0, W0, B0, A0) = solveMesonWfs(ρ_sp, ρ_vp, ρ_sn, ρ_vn, r_max, divs, 1)
+(Φ0, W0, B0, A0) = solveMesonWfs(ρ_sp, ρ_vp, ρ_sn, ρ_vn, r_max, divs, 1; initial_sol=(Φ0_bench, W0_bench, B0_bench, A0_bench))
 
-plot(xs, hcat(Φ0, W0, B0, A0), layout=4, label=["Φ0" "W0" "B0" "A0"])
-xlabel!("r (fm)")
-
-# test data generated from Hartree.f
-# format: x S(x) V(x) R(x) A(x)
-test_data = readdlm("test/Pb208FldsFSUGarnet.csv") 
-xs_bench = test_data[:, 1]
-Φ0_bench = test_data[:, 2]
-W0_bench = test_data[:, 3]
-B0_bench = test_data[:, 4]
-A0_bench = test_data[:, 5]
-
-plot!(xs_bench, hcat(Φ0_bench, W0_bench, B0_bench, A0_bench), layout=4, label=["Φ0 benchmark" "W0 benchmark" "B0 benchmark" "A0 benchmark"])
+plot!(xs, hcat(Φ0, W0, B0, A0), layout=4, label=["Φ0" "W0" "B0" "A0"])
+xlabel!("r (fm)")
\ No newline at end of file