Distinguish in-place functions

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)
-
+        @time "Proton densities" solveNucleonDensity!(true, s)
-        # protons
+        @time "Neutron densities" solveNucleonDensity!(false, s)
-        @time "Proton densities" (s.ρ_sp, s.ρ_vp) = solveNucleonDensity(true, s)
-
-        # neutrons
-        @time "Neutron densities" (s.ρ_sn, s.ρ_vn) = solveNucleonDensity(false, s)
 
         if monitor_plot
             for s in p.series_list

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"

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)

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")