From 5680b7e489673c10e46925bf7b676e4cd29ccee4 Mon Sep 17 00:00:00 2001 From: Nuwan Yapa Date: Tue, 25 Feb 2025 16:37:07 -0500 Subject: [PATCH] Distinguish in-place functions --- NuclearRMF.jl | 11 +++-------- nucleons.jl | 15 +++++++++++---- test/Pb208.jl | 2 +- test/Pb208_nucleon_dens.jl | 4 ++-- 4 files changed, 17 insertions(+), 15 deletions(-) diff --git a/NuclearRMF.jl b/NuclearRMF.jl index 646aa45..9877aad 100644 --- a/NuclearRMF.jl +++ b/NuclearRMF.jl @@ -9,7 +9,7 @@ total_E(s::system) = -(nucleon_E(s) + meson_E(s)) 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; reinitialize_densities=true, monitor_print=true, monitor_plot=false) +function solve_system!(s::system; reinitialize_densities=true, monitor_print=true, monitor_plot=false) if reinitialize_densities dens_guess = Woods_Saxon.(rs(s)) @. s.ρ_sp = s.Z * dens_guess @@ -25,14 +25,9 @@ function solve_system(s::system; reinitialize_densities=true, monitor_print=true previous_E_per_A = NaN while true - # mesons @time "Meson fields" solveMesonFields!(s, isnan(previous_E_per_A) ? 50 : 15) - - # protons - @time "Proton densities" (s.ρ_sp, s.ρ_vp) = solveNucleonDensity(true, s) - - # neutrons - @time "Neutron densities" (s.ρ_sn, s.ρ_vn) = solveNucleonDensity(false, s) + @time "Proton densities" solveNucleonDensity!(true, s) + @time "Neutron densities" solveNucleonDensity!(false, s) if monitor_plot for s in p.series_list diff --git a/nucleons.jl b/nucleons.jl index 5e7a73f..e57e0f5 100644 --- a/nucleons.jl +++ b/nucleons.jl @@ -144,7 +144,7 @@ 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(p::Bool, s::system) +function calculateNucleonDensity(p::Bool, s::system)::Tuple{Vector{Float64}, Vector{Float64}} spectrum = p ? s.p_spectrum : s.n_spectrum (κs, Es, occs) = (spectrum.κ, spectrum.E, spectrum.occ) @@ -166,9 +166,9 @@ function calculateNucleonDensity(p::Bool, s::system) return (ρ_s, ρ_v) end -"Solve the Dirac equation and calculate scalar and vector densities of a nucleon species where +"For a nucleon species, solve the Dirac equation and save the spectrum & densities in-place where the other parameters are defined above" -function solveNucleonDensity(p::Bool, s::system, E_min=850.0, E_max=938.0) +function solveNucleonDensity!(p::Bool, s::system, E_min=850.0, E_max=938.0) κs, Es = findAllOrbitals(p, s, E_min, E_max) spec = fillNucleons(Z_or_N(s, p), κs, Es) if p @@ -176,7 +176,14 @@ function solveNucleonDensity(p::Bool, s::system, E_min=850.0, E_max=938.0) else s.n_spectrum = spec end - return calculateNucleonDensity(p, s) + (ρ_s, ρ_v) = calculateNucleonDensity(p, s) + if p + s.ρ_sp = ρ_s + s.ρ_vp = ρ_v + else + s.ρ_sn = ρ_s + s.ρ_vn = ρ_v + end end "Total energy of filled nucleons in the system" diff --git a/test/Pb208.jl b/test/Pb208.jl index 442bb53..cae50b9 100644 --- a/test/Pb208.jl +++ b/test/Pb208.jl @@ -7,4 +7,4 @@ divs = 400 s = system(Z, N, r_max, divs) -solve_system(s; monitor_plot=true) +solve_system!(s; monitor_plot=true) diff --git a/test/Pb208_nucleon_dens.jl b/test/Pb208_nucleon_dens.jl index d6887fc..1e11111 100644 --- a/test/Pb208_nucleon_dens.jl +++ b/test/Pb208_nucleon_dens.jl @@ -21,12 +21,12 @@ s.W0 = Vs s.B0 = Rs s.A0 = As -(s.ρ_sp, s.ρ_vp) = solveNucleonDensity(true, s) +solveNucleonDensity!(true, s) p_sp = plot(rs(s), s.ρ_sp, xlabel="r (fm)", label="ρₛₚ(r) calculated") p_vp = plot(rs(s), s.ρ_vp, xlabel="r (fm)", label="ρᵥₚ(r) calculated") -(s.ρ_sn, s.ρ_vn) = solveNucleonDensity(false, s) +solveNucleonDensity!(false, s) p_sn = plot(rs(s), s.ρ_sn, xlabel="r (fm)", label="ρₛₙ(r) calculated") p_vn = plot(rs(s), s.ρ_vn, xlabel="r (fm)", label="ρᵥₙ(r) calculated")