From 43e91bf78e676fa3e64c51ec65d2a4f4d5834444 Mon Sep 17 00:00:00 2001 From: Nuwan Yapa Date: Tue, 25 Feb 2025 16:12:36 -0500 Subject: [PATCH] struct spectrum --- NuclearRMF.jl | 8 ++------ nucleons.jl | 24 ++++++++++++++++-------- system.jl | 23 +++++++++++++---------- test/Linear_nucleon_spectrum.jl | 4 ++-- test/Pb208_nucleon_spectrum.jl | 4 ++-- 5 files changed, 35 insertions(+), 28 deletions(-) diff --git a/NuclearRMF.jl b/NuclearRMF.jl index ab951fa..646aa45 100644 --- a/NuclearRMF.jl +++ b/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) - 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) + @time "Proton densities" (s.ρ_sp, s.ρ_vp) = solveNucleonDensity(true, s) # neutrons - @time "Neutron spectrum" (s.κ_n, s.E_n) = findAllOrbitals(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) + @time "Neutron densities" (s.ρ_sn, s.ρ_vn) = solveNucleonDensity(false, s) if monitor_plot for s in p.series_list diff --git a/nucleons.jl b/nucleons.jl index 095cc2b..5e7a73f 100644 --- a/nucleons.jl +++ b/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" -function fillNucleons(Z_or_N::Int, κs, Es) +"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)::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" - return occ + @assert Z_or_N == 0 "All orbitals could not be filled" + 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) - return calculateNucleonDensity(κs, Es, occs, p, s) + spec = fillNucleons(Z_or_N(s, p), κs, Es) + 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)) diff --git a/system.jl b/system.jl index 6df5e62..c557325 100644 --- a/system.jl +++ b/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::Vector{Int} - 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} + p_spectrum::spectrum + n_spectrum::spectrum Φ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) diff --git a/test/Linear_nucleon_spectrum.jl b/test/Linear_nucleon_spectrum.jl index 12cb119..0c7d036 100644 --- a/test/Linear_nucleon_spectrum.jl +++ b/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]) diff --git a/test/Pb208_nucleon_spectrum.jl b/test/Pb208_nucleon_spectrum.jl index 538ec19..da55a59 100644 --- a/test/Pb208_nucleon_spectrum.jl +++ b/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])