Rough subroutine for self-consistent solution

Project.toml

@@ -2,4 +2,5 @@
 DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
 Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+PolyLog = "85e3b03c-9856-11eb-0374-4dc1f8670e7f"
 Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"

system.jl

@@ -0,0 +1,70 @@
+using Interpolations, PolyLog
+include("nucleons.jl")
+include("mesons.jl")
+
+"Defines a nuclear system to be solved"
+struct system
+    Z::Int
+    N::Int
+    r_max::Float64
+    divs::Int
+end
+
+"Get mass number of nucleus"
+A(s::system) = s.Z + s.N
+
+"Get r values in the mesh"
+rs(s::system) = range(0, s.r_max, length=s.divs+1)
+
+"Get Δr value for the mesh"
+Δr(s::system) = s.r_max / s.divs
+
+"Create an empty array for the size of the mesh"
+zero_array(s::system) = zeros(1 + s.divs)
+
+"Normalized Woods-Saxon form used for constructing an initial solution"
+Woods_Saxon(r::Float64; R::Float64=7.0, a::Float64=0.5) = -1 / (8π * a^3 * reli3(-exp(R / a)) * (1 + exp((r - R) / a)))
+
+"Run the full program by self-consistent solution of nucleon and meson densities"
+function solve_system(s::system, initial_dens=nothing, initial_flds=(zeros(1 + s.divs) for _ in 1:4))
+    if isnothing(initial_dens)
+        dens_guess = Woods_Saxon.(rs(s))
+        ρ_sp = s.Z .* dens_guess
+        ρ_vp = s.Z .* dens_guess
+        ρ_sn = s.N .* dens_guess
+        ρ_vn = s.N .* dens_guess
+    else
+        (ρ_sp, ρ_vp, ρ_sn, ρ_vn) = initial_dens
+    end
+
+    (Φ0s, W0s, B0s, A0s) = initial_flds
+
+    E_total_previous = NaN
+
+    while true
+        (Φ0s, W0s, B0s, A0s) = solveMesonWfs(ρ_sp, ρ_vp, ρ_sn, ρ_vn, s.r_max, s.divs, 50; initial_sol = (Φ0s, W0s, B0s, A0s))
+
+        S_interp = linear_interpolation(rs(s), Φ0s)
+        V_interp = linear_interpolation(rs(s), W0s)
+        R_interp = linear_interpolation(rs(s), B0s)
+        A_interp = linear_interpolation(rs(s), A0s)
+
+        # protons
+        κs_p, Es_p = findAllOrbitals(true, S_interp, V_interp, R_interp, A_interp, s.r_max)
+        occs_p = fillNucleons(s.Z, κs_p, Es_p)
+        (ρ_sp, ρ_vp) = calculateNucleonDensity(κs_p, Es_p, occs_p, true, S_interp, V_interp, R_interp, A_interp, s.r_max, s.divs)
+
+        # neutrons
+        κs_n, Es_n = findAllOrbitals(false, S_interp, V_interp, R_interp, A_interp, s.r_max)
+        occs_n = fillNucleons(s.N, κs_n, Es_n)
+        (ρ_sp, ρ_vp) = calculateNucleonDensity(κs_n, Es_n, occs_n, false, S_interp, V_interp, R_interp, A_interp, s.r_max, s.divs)
+
+        # total binding energy of nucleons
+        E_total = sum(M_p .- Es_p) + sum(M_n .- Es_n)
+        println("Total binding E per nucleon = $(E_total/A(s))")
+
+        # check convergence
+        abs(E_total_previous - E_total) < 0.1 && break
+        E_total_previous = E_total
+    end
+end

test/Pb208.jl

@@ -0,0 +1,5 @@
+include("../system.jl")
+
+s = system(82, 126, 20.0, 400)
+
+solve_system(s)