Fields and densities read off from s::system

mesons.jl

@@ -38,14 +38,13 @@ 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,
-    An initial guess initial_sol=(Φ0, W0, B0, A0) can be provided to speed up convergence (permuting!).
+    the inital solutions are read from s and the final solutions are saved in-place.
 Reference: P. Giuliani, K. Godbey, E. Bonilla, F. Viens, and J. Piekarewicz, Frontiers in Physics 10, (2023)"
-function solveMesonWfs(ρ_sp, ρ_vp, ρ_sn, ρ_vn, s::system, iterations=10; initial_sol=(zeros(1 + divs) for _ in 1:4))
+function solveMesonFields!(s::system, iterations=10)
-    (Φ0, W0, B0, A0) = initial_sol
+    (Φ0, W0, B0, A0) = (s.Φ0, s.W0, s.B0, s.A0)
+    (ρ_sp, ρ_vp, ρ_sn, ρ_vn) = (s.ρ_sp, s.ρ_vp, s.ρ_sn, s.ρ_vn)
+
     (src_Φ0, src_W0, src_B0, src_A0) = (zero_array(s) for _ in 1:4)
 
     # A0 doesn't need iterations

nucleons.jl

@@ -27,7 +27,9 @@ 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 solveNucleonWf(κ, p, E, Φ0, W0, B0, A0, s::system; shooting=true, normalize=true)
+function solveNucleonWf(κ, p, E, s::system; shooting=true, normalize=true)
+    (Φ0, W0, B0, A0) = fields_interp(s)
+
     if shooting
         @assert s.divs % 2 == 0 "divs must be an even number when shooting=true"
         prob = ODEProblem(dirac!, [0, 1], (s.r_max, s.r_max / 2))
@@ -62,8 +64,9 @@ end
 "Returns a function that solves the Dirac equation and returns g(r=r_max) where
     r_max is the outer boundary in fm,
     the other parameters are the same from dirac!(...)."
-function boundaryValueFunc(κ, p, Φ0, W0, B0, A0, s::system; dtype=Float64, algo=Tsit5())
+function boundaryValueFunc(κ, p, s::system; dtype=Float64, algo=Tsit5())
     prob = ODEProblem(dirac!, convert.(dtype, [0, 1]), (0, s.r_max))
+    (Φ0, W0, B0, A0) = fields_interp(s)
     func(E) = solve(prob, algo, p=(κ, p, E, Φ0, W0, B0, A0), saveat=[s.r_max], save_idxs=[1])[1, 1]
     return func
 end
@@ -71,21 +74,21 @@ end
 "Find all bound energies between E_min (=0) and E_max (=mass) where
     tol_digits determines the precision for root finding and the threshold for identifying duplicates,
     the other parameters are the same from dirac!(...)."
-function findEs(κ, p, Φ0, W0, B0, A0, s::system, E_min=0, E_max=(p ? M_p : M_n), tol_digits=5)
+function findEs(κ, p, s::system, E_min=0, E_max=(p ? M_p : M_n), tol_digits=5)
-    func = boundaryValueFunc(κ, p, Φ0, W0, B0, A0, s)
+    func = boundaryValueFunc(κ, p, s)
     Es = find_zeros(func, (E_min, E_max); xatol=1/10^tol_digits)
     return unique(E -> round(E; digits=tol_digits), Es)
 end
 
 "Find all orbitals and return two lists containing κ values and corresponding energies for a single species where
     the other parameters are defined above"
-function findAllOrbitals(p, Φ0, W0, B0, A0, s::system, E_min=0, E_max=(p ? M_p : M_n))
+function findAllOrbitals(p, s::system, E_min=0, E_max=(p ? M_p : M_n))
     κs = Int[]
     Es = Float64[]
     # start from κ=-1 and go both up and down
     for direction in [-1, 1]
         for κ in direction * (1:100) # cutoff is 100
-            new_Es = findEs(κ, p, Φ0, W0, B0, A0, s, E_min, E_max)
+            new_Es = findEs(κ, p, s, E_min, E_max)
             if isempty(new_Es); break; end
             append!(Es, new_Es)
             append!(κs, fill(κ, length(new_Es)))
@@ -120,11 +123,11 @@ l_κ(κ::Int) = abs(κ) - (κ < 0) # since true = 1 and false = 0
 "Calculate scalar and vector densities of a nucleon species on [0,r_max] divided into (divs+1) points and returns them as vectors (ρ_s, ρ_v) where
     the arrays κs, Es, occs tabulate the energies and occupation numbers corresponding to each κ,
     the other parameters are defined above"
-function calculateNucleonDensity(κs, Es, occs, p, Φ0, W0, B0, A0, s::system)
+function calculateNucleonDensity(κs, Es, occs, p, s::system)
     ρr2 = zeros(2, s.divs + 1) # ρ×r² values
     
     for (κ, E, occ) in zip(κs, Es, occs)
-        wf = solveNucleonWf(κ, p, E, Φ0, W0, B0, A0, s; shooting=true, normalize=true)
+        wf = solveNucleonWf(κ, p, E, s; shooting=true, normalize=true)
         wf2 = wf .* wf
         ρr2 += (occ / (4 * pi)) * wf2 # 2j+1 factor is accounted in the occupancy number
     end
@@ -141,8 +144,8 @@ end
 
 "Solve the Dirac equation and calculate scalar and vector densities of a nucleon species where
     the other parameters are defined above"
-function solveNucleonDensity(p, Φ0, W0, B0, A0, s::system, E_min=800, E_max=939)
+function solveNucleonDensity(p, s::system, E_min=800, E_max=939)
-    κs, Es = findAllOrbitals(p, Φ0, W0, B0, A0, s, E_min, E_max)
+    κs, Es = findAllOrbitals(p, s, E_min, E_max)
     occs = fillNucleons(Z_or_N(s, p), κs, Es)
-    return calculateNucleonDensity(κs, Es, occs, p, Φ0, W0, B0, A0, s)
+    return calculateNucleonDensity(κs, Es, occs, p, s)
 end

system.jl

@@ -18,11 +18,11 @@ mutable struct system
     ρ_sn::Vector{Float64}
     ρ_vn::Vector{Float64}
 
-    "Dummy struct to define the mesh"
-    system(r_max, divs) = new(0, 0, r_max, divs, [Float64[] for _ in 1:8]...)
-
     "Initialize an unsolved system"
     system(Z, N, r_max, divs) = new(Z, N, r_max, divs, [zeros(1 + divs) for _ in 1:8]...)
+
+    "Dummy struct to define the mesh"
+    system(r_max, divs) = system(0, 0, r_max, divs)
 end
 
 "Get mass number of nucleus"
@@ -73,18 +73,17 @@ function solve_system(s::system; reinitialize_densities=true, monitor_print=true
 
     while true
         # mesons
-        @time "Meson fields" solveMesonWfs(s.ρ_sp, s.ρ_vp, s.ρ_sn, s.ρ_vn, s, isnan(E_total_previous) ? 50 : 5; initial_sol = (s.Φ0, s.W0, s.B0, s.A0))
+        @time "Meson fields" solveMesonFields!(s, isnan(E_total_previous) ? 50 : 5)
-        (S_interp, V_interp, R_interp, A_interp) = fields_interp(s)
 
         # protons
-        @time "Proton spectrum" (κs_p, Es_p) = findAllOrbitals(true, S_interp, V_interp, R_interp, A_interp, s)
+        @time "Proton spectrum" (κs_p, Es_p) = findAllOrbitals(true, s)
         occs_p = fillNucleons(s.Z, κs_p, Es_p)
-        @time "Proton densities" (s.ρ_sp, s.ρ_vp) = calculateNucleonDensity(κs_p, Es_p, occs_p, true, S_interp, V_interp, R_interp, A_interp, s)
+        @time "Proton densities" (s.ρ_sp, s.ρ_vp) = calculateNucleonDensity(κs_p, Es_p, occs_p, true, s)
 
         # neutrons
-        @time "Neutron spectrum" (κs_n, Es_n) = findAllOrbitals(false, S_interp, V_interp, R_interp, A_interp, s)
+        @time "Neutron spectrum" (κs_n, Es_n) = findAllOrbitals(false, s)
         occs_n = fillNucleons(s.N, κs_n, Es_n)
-        @time "Neutron densities" (s.ρ_sn, s.ρ_vn) = calculateNucleonDensity(κs_n, Es_n, occs_n, false, S_interp, V_interp, R_interp, A_interp, s)
+        @time "Neutron densities" (s.ρ_sn, s.ρ_vn) = calculateNucleonDensity(κs_n, Es_n, occs_n, false, s)
 
         if monitor_plot
             for s in p.series_list

test/Pb208_meson_flds.jl

@@ -16,16 +16,17 @@ plot(xs_bench, hcat(Φ0_bench, W0_bench, B0_bench, A0_bench), layout=4, label=["
 # 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
 s = system(r_max, divs)
 
-(Φ0, W0, B0, A0) = solveMesonWfs(ρ_sp, ρ_vp, ρ_sn, ρ_vn, s, 200) # many iterations needed without an initial solution
+s.ρ_sn = test_data[:, 2]
+s.ρ_vn = test_data[:, 3]
+s.ρ_sp = test_data[:, 5]
+s.ρ_vp = test_data[:, 6]
+
+(Φ0, W0, B0, A0) = solveMesonFields!(s, 200) # many iterations needed without an initial solution
 
 plot!(xs, hcat(Φ0, W0, B0, A0), layout=4, label=["Φ0" "W0" "B0" "A0"])
 xlabel!("r (fm)")

test/Pb208_nucleon_basic.jl

@@ -21,17 +21,15 @@ s.W0 = Vs
 s.B0 = Rs
 s.A0 = As
 
-(S_interp, V_interp, R_interp, A_interp) = fields_interp(s)
-
 E_min = 850
 E_max = 939
 
-boundEs = findEs(κ, p, S_interp, V_interp, R_interp, A_interp, s, E_min, E_max)
+boundEs = findEs(κ, p, s, E_min, E_max)
 
 binding_Es = (p ? M_p : M_n) .- boundEs
 println("binding energies = $binding_Es")
 
-func = boundaryValueFunc(κ, p, S_interp, V_interp, R_interp, A_interp, s)
+func = boundaryValueFunc(κ, p, s)
 Es = collect(E_min:0.5:E_max)
 boundaryVals = [func(E)^2 for E in Es]
 

test/Pb208_nucleon_dens.jl

@@ -21,14 +21,12 @@ s.W0 = Vs
 s.B0 = Rs
 s.A0 = As
 
-(S_interp, V_interp, R_interp, A_interp) = fields_interp(s)
+(ρ_sp, ρ_vp) = solveNucleonDensity(true, s)
-
-(ρ_sp, ρ_vp) = solveNucleonDensity(true, S_interp, V_interp, R_interp, A_interp, s)
 
 p_sp = plot(rs(s), ρ_sp, xlabel="r (fm)", label="ρₛₚ(r) calculated")
 p_vp = plot(rs(s), ρ_vp, xlabel="r (fm)", label="ρᵥₚ(r) calculated")
 
-(ρ_sn, ρ_vn) = solveNucleonDensity(false, S_interp, V_interp, R_interp, A_interp, s)
+(ρ_sn, ρ_vn) = solveNucleonDensity(false, s)
 
 p_sn = plot(rs(s), ρ_sn, xlabel="r (fm)", label="ρₛₙ(r) calculated")
 p_vn = plot(rs(s), ρ_vn, xlabel="r (fm)", label="ρᵥₙ(r) calculated")

test/Pb208_nucleon_spectrum.jl

@@ -22,12 +22,10 @@ s.W0 = Vs
 s.B0 = Rs
 s.A0 = As
 
-(S_interp, V_interp, R_interp, A_interp) = fields_interp(s)
-
 E_min = 800
 E_max = 939
 
-κs, Es = findAllOrbitals(p, S_interp, V_interp, R_interp, A_interp, s, E_min, E_max)
+κs, Es = findAllOrbitals(p, s, E_min, E_max)
 Ebinds = (p ? M_p : M_n) .- Es
 occ = fillNucleons(Z_or_N(s, p), κs, Es)
 

test/Pb208_nucleon_wf.jl

@@ -21,15 +21,13 @@ s.W0 = Vs
 s.B0 = Rs
 s.A0 = As
 
-(S_interp, V_interp, R_interp, A_interp) = fields_interp(s)
-
 E_min = 800
 E_max = 939
 
-groundE = findEs(κ, p, S_interp, V_interp, R_interp, A_interp, s, E_min, E_max) |> minimum
+groundE = findEs(κ, p, s, E_min, E_max) |> minimum
 println("ground state E = $groundE")
 
-wf = solveNucleonWf(κ, p, groundE, S_interp, V_interp, R_interp, A_interp, s)
+wf = solveNucleonWf(κ, p, groundE, s)
 gs = wf[1, :]
 fs = wf[2, :]