struct spectrum

NuclearRMF.jl

@@ -29,14 +29,10 @@ function solve_system(s::system; reinitialize_densities=true, monitor_print=true
         @time "Meson fields" solveMesonFields!(s, isnan(previous_E_per_A) ? 50 : 15)
 
         # protons
-        @time "Proton spectrum" (s.κ_p, s.E_p) = findAllOrbitals(true, s)
+        @time "Proton densities" (s.ρ_sp, s.ρ_vp) = solveNucleonDensity(true, s)
-        s.occ_p = fillNucleons(s.Z, s.κ_p, s.E_p)
-        @time "Proton densities" (s.ρ_sp, s.ρ_vp) = calculateNucleonDensity(s.κ_p, s.E_p, s.occ_p, true, s)
 
         # neutrons
-        @time "Neutron spectrum" (s.κ_n, s.E_n) = findAllOrbitals(false, s)
+        @time "Neutron densities" (s.ρ_sn, s.ρ_vn) = solveNucleonDensity(false, s)
-        s.occ_n = fillNucleons(s.N, s.κ_n, s.E_n)
-        @time "Neutron densities" (s.ρ_sn, s.ρ_vn) = calculateNucleonDensity(s.κ_n, s.E_n, s.occ_n, false, s)
 
         if monitor_plot
             for s in p.series_list

nucleons.jl

@@ -118,8 +118,8 @@ function findAllOrbitals(p::Bool, s::system, E_min=850.0, E_max=938.0)
     return (κs, Es)
 end
 
-"For a given list of κ values with corresponding energies, attempt to fill Z_or_N lowest lying orbitals and return occupancy numbers"
+"For a given list of κ values with corresponding energies, attempt to fill Z_or_N lowest lying orbitals and return the spectrum"
-function fillNucleons(Z_or_N::Int, κs, Es)
+function fillNucleons(Z_or_N::Int, κs, Es)::spectrum
     sort_i = sortperm(Es)
 
     occ = zeros(Int, length(κs))
@@ -131,8 +131,8 @@ function fillNucleons(Z_or_N::Int, κs, Es)
         Z_or_N -= occ[i]
     end
 
-    Z_or_N == 0 || @warn "All orbitals could not be filled"
+    @assert Z_or_N == 0 "All orbitals could not be filled"
-    return occ
+    return spectrum(κs, Es, occ)
 end
 
 "Total angular momentum j for a given κ value"
@@ -144,7 +144,10 @@ 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::Bool, s::system)
+function calculateNucleonDensity(p::Bool, s::system)
+    spectrum = p ? s.p_spectrum : s.n_spectrum
+    (κs, Es, occs) = (spectrum.κ, spectrum.E, spectrum.occ)
+    
     ρr2 = zeros(2, s.divs + 1) # ρ×r² values
     
     for (κ, E, occ) in zip(κs, Es, occs)
@@ -167,9 +170,14 @@ end
     the other parameters are defined above"
 function solveNucleonDensity(p::Bool, s::system, E_min=850.0, E_max=938.0)
     κs, Es = findAllOrbitals(p, s, E_min, E_max)
-    occs = fillNucleons(Z_or_N(s, p), κs, Es)
+    spec = fillNucleons(Z_or_N(s, p), κs, Es)
-    return calculateNucleonDensity(κs, Es, occs, p, s)
+    if p
+        s.p_spectrum = spec
+    else
+        s.n_spectrum = spec
+    end
+    return calculateNucleonDensity(p, s)
 end
 
 "Total energy of filled nucleons in the system"
-nucleon_E(s::system) = sum(s.occ_p .* (s.E_p .- M_p)) + sum(s.occ_n .* (s.E_n .- M_n))
+nucleon_E(s::system) = sum(s.p_spectrum.occ .* (s.p_spectrum.E .- M_p)) + sum(s.n_spectrum.occ .* (s.n_spectrum.E .- M_n))

system.jl

@@ -1,5 +1,15 @@
 using Interpolations
 
+"Tabulates a nucleon spectrum (protons or neutrons) containing κ and occupancy"
+struct spectrum
+    κ::Vector{Int}
+    E::Vector{Float64}
+    occ::Vector{Int}
+end
+
+"Initializes an unfilled spectrum"
+unfilled_spectrum() = spectrum(Int[], Float64[], Int[])
+
 "Defines a nuclear system containing relevant parameters and meson/nucleon densities"
 mutable struct system
     Z::Int
@@ -8,15 +18,8 @@ mutable struct system
     r_max::Float64
     divs::Int
 
-    # single particle energies and corresponding occupancy numbers for protons
+    p_spectrum::spectrum
-    κ_p::Vector{Int}
+    n_spectrum::spectrum
-    E_p::Vector{Float64}
-    occ_p::Vector{Int}
-
-    # single particle energies and corresponding occupancy numbers for neutrons
-    κ_n::Vector{Int}
-    E_n::Vector{Float64}
-    occ_n::Vector{Int}
 
     Φ0::Vector{Float64}
     W0::Vector{Float64}
@@ -29,7 +32,7 @@ mutable struct system
     ρ_vn::Vector{Float64}
 
     "Initialize an unsolved system"
-    system(Z, N, r_max, divs) = new(Z, N, r_max, divs, Int[], Float64[], Int[], Int[], Float64[], Int[], [zeros(1 + divs) for _ in 1:8]...)
+    system(Z, N, r_max, divs) = new(Z, N, r_max, divs, unfilled_spectrum(), unfilled_spectrum(), [zeros(1 + divs) for _ in 1:8]...)
 
     "Dummy struct to define the mesh"
     system(r_max, divs) = system(0, 0, r_max, divs)

test/Linear_nucleon_spectrum.jl

@@ -22,7 +22,7 @@ s.A0 = As
 
 κs, Es = findAllOrbitals(p, s)
 Ebinds = (p ? M_p : M_n) .- Es
-occ = fillNucleons(Z_or_N(s, p), κs, Es)
+spec = fillNucleons(Z_or_N(s, p), κs, Es)
 
 # format: proton, kappa, filling, gnodes, fnodes, Ebind
 bench_data, _ = readdlm("test/LinearSpectrum.csv", ','; header=true)
@@ -33,7 +33,7 @@ bench_filling = bench_data[mask, 3]
 bench_occ = Int.(bench_filling .* (2 .* j_κ.(Int.(bench_κs)) .+ 1))
 
 scatter(κs, Ebinds, label=(p ? "proton" : "neutron") * " spectrum")
-annotate!(κs .+ 0.3, Ebinds .+ 1, [(o, 9) for o in occ])
+annotate!(κs .+ 0.3, Ebinds .+ 1, [(o, 9) for o in spec.occ])
 
 scatter!(bench_κs, bench_Ebinds, label="benchmark"; markershape=:cross, markercolor="red")
 annotate!(bench_κs .+ 0.3, bench_Ebinds .- 1, [(o, 9, :red) for o in bench_occ])

test/Pb208_nucleon_spectrum.jl

@@ -24,7 +24,7 @@ s.A0 = As
 
 κs, Es = findAllOrbitals(p, s)
 Ebinds = (p ? M_p : M_n) .- Es
-occ = fillNucleons(Z_or_N(s, p), κs, Es)
+spec = fillNucleons(Z_or_N(s, p), κs, Es)
 
 # format: proton, kappa, filling, gnodes, fnodes, Ebind
 bench_data, _ = readdlm("test/Pb208Spectrum.csv", ','; header=true)
@@ -35,7 +35,7 @@ bench_filling = bench_data[mask, 3]
 bench_occ = Int.(bench_filling .* (2 .* j_κ.(Int.(bench_κs)) .+ 1))
 
 scatter(κs, Ebinds, label=(p ? "proton" : "neutron") * " spectrum")
-annotate!(κs .+ 0.3, Ebinds .+ 1, [(o, 9) for o in occ])
+annotate!(κs .+ 0.3, Ebinds .+ 1, [(o, 9) for o in spec.occ])
 
 scatter!(bench_κs, bench_Ebinds, label="benchmark"; markershape=:cross, markercolor="red")
 annotate!(bench_κs .+ 0.3, bench_Ebinds .- 1, [(o, 9, :red) for o in bench_occ])