Model parameters are no longer hard-coded

NuclearRMF.jl

@@ -3,45 +3,66 @@ include("nucleons.jl")
 include("mesons.jl")
 
 "Total binding energy of the system"
-total_E(s::system) = -(nucleon_E(s) + meson_E(s))
+function total_E(s::system)
+    E_cm = (3/4) * 41.0 * A(s)^(-1/3) # Center-of-mass correction [Ring and Schuck, Sec. 2.3]
+    return -(nucleon_E(s) + meson_E(s)) + E_cm
+end
 
 "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)))
 
+"Conditionally toggle @time if enable_time is true"
+macro conditional_time(label, expr)
+    quote
+        if $(esc(:enable_time))
+            @time $label $(esc(expr))
+        else
+            $(esc(expr))
+        end
+    end
+end
+
 "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, print_E=true, print_time=false, live_plots=false)
     if reinitialize_densities
-        dens_guess = Woods_Saxon.(s.r_mesh.r)
+        dens_guess = Woods_Saxon.(rs(s))
         @. s.ρ_sp = s.Z * dens_guess
         @. s.ρ_vp = s.Z * dens_guess
         @. s.ρ_sn = s.N * dens_guess
         @. s.ρ_vn = s.N * dens_guess
     end
 
-    if monitor_plot
+    if live_plots
         p = plot(legends=false, size=(1024, 768), layout=(2, 4), title=["ρₛₚ" "ρᵥₚ" "ρₛₙ" "ρᵥₙ" "Φ₀" "W₀" "B₀" "A₀"])
     end
 
     previous_E_per_A = NaN
 
     while true
-        @time "Meson fields" solveMesonFields!(s, isnan(previous_E_per_A) ? 50 : 15)
+        enable_time = print_time
-        @time "Proton densities" solveNucleonDensity!(true, s)
+        @conditional_time "Meson fields" solveMesonFields!(s, isnan(previous_E_per_A) ? 50 : 15)
-        @time "Neutron densities" solveNucleonDensity!(false, s)
+        @conditional_time "Proton densities" solveNucleonDensity!(true, s)
+        @conditional_time "Neutron densities" solveNucleonDensity!(false, s)
 
-        if monitor_plot
+        if live_plots
             for s in p.series_list
                 s.plotattributes[:linecolor] = :gray
             end
-            plot!(p, s.r_mesh.r, hcat(s.ρ_sp, s.ρ_vp, s.ρ_sn, s.ρ_vn, s.Φ0, s.W0, s.B0, s.A0), linecolor=:red)
+            plot!(p, rs(s), hcat(s.ρ_sp, s.ρ_vp, s.ρ_sn, s.ρ_vn, s.Φ0, s.W0, s.B0, s.A0), linecolor=:red)
             display(p)
         end
 
         E_per_A = total_E(s) / A(s)
-        monitor_print && println("Total binding E per nucleon = $E_per_A")
+        print_E && println("Total binding E per nucleon = $E_per_A")
 
         # check convergence
         abs(previous_E_per_A - E_per_A) < 0.0001 && break
         previous_E_per_A = E_per_A
     end
 end
+
+"Calculate RMS radius from density"
+rms_radius(p::Bool, s::system) = 4pi * Δr(s) * sum((rs(s) .^ 4) .* (p ? s.ρ_vp : s.ρ_vn)) / (p ? s.Z : s.N) |> sqrt
+
+"Calculate neutron skin thickness"
+R_skin(s::system) = rms_radius(false, s) - rms_radius(true, s)

Project.toml

@@ -1,6 +1,5 @@
 [deps]
 DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
-FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838"
 Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 PolyLog = "85e3b03c-9856-11eb-0374-4dc1f8670e7f"

common.jl

@@ -1,2 +1,17 @@
 const ħc = 197.33 # MeVfm
 const r_reg = 1E-8 # fm # regulator for R
+
+"Integrate a uniformly discretized function f using Simpson's rule where h is the step size and coefficient is an optional scaling factor."
+function simpsons_integrate(f::AbstractVector{Float64}, h::Float64; coefficient::Float64 = 1.0)
+    @assert length(f) % 2 == 1 "Number of mesh divisions must be even for Simpson's rule"
+    s = sum(enumerate(f)) do (i, fi)
+        if i == 1 || i == length(f)
+            return fi
+        elseif i % 2 == 0
+            return 4fi
+        else
+            return 2fi
+        end
+    end
+    return (h/3) * coefficient * s
+end

mesons.jl

@@ -1,21 +1,6 @@
 include("common.jl")
 include("system.jl")
 
-# Values defined in C. J. Horowitz and J. Piekarewicz, Phys. Rev. Lett. 86, 5647 (2001)
-# Values taken from Hartree.f (FSUGarnet)
-const m_s = 496.939473213388 # MeV/c2
-const m_ω = 782.5 # MeV/c2
-const m_ρ = 763.0 # MeV/c2
-const m_γ = 0.000001000 # MeV/c2 # should be 0?
-const g2_s = 110.349189097820 # dimensionless
-const g2_v = 187.694676506801 # dimensionless
-const g2_ρ = 192.927428365698 # dimensionless
-const g2_γ = 0.091701236 # dimensionless # equal to 4πα
-const κ_ss = 3.260178893447 # MeV
-const λ = -0.003551486718 # dimensionless # LambdaSS
-const ζ = 0.023499504053 # dimensionless # LambdaVV
-const Λv = 0.043376933644 # dimensionless # LambdaVR
-
 "Green's function for Klein-Gordon equation in natural units"
 greensFunctionKG(m, r, rp) = -1 / (m * (r + r_reg) * (rp + r_reg)) * sinh(m * min(r, rp)) * exp(-m * max(r, rp))
 
@@ -25,19 +10,24 @@ greensFunctionP(r, rp) = -1 / (max(r, rp) + r_reg)
 "Solve the Klein-Gordon equation (or Poisson's equation when m=0) and return an array in MeV for a source array given in fm⁻³ where
     m is the mass of the meson in MeV/c2."
 function solveKG(m, source, s::system)
-    int_measure = ħc .* s.r_mesh.w .* (s.r_mesh.r .^ 2)
+
+    @assert s.divs % 2 == 0 "Number of mesh divisions must be even for Simpson's rule"
+    simpsons_weights = (Δr(s)/3) .* [1; repeat([2,4], s.divs ÷ 2)[2:end]; 1]
+    int_measure = ħc .* simpsons_weights .* rs(s) .^ 2
 
     greensFunction = m == 0 ? greensFunctionP : (r, rp) -> greensFunctionKG(m / ħc, r, rp)
-    greenMat = greensFunction.(s.r_mesh.r, transpose(s.r_mesh.r))
+    greenMat = greensFunction.(rs(s), transpose(rs(s)))
 
     return greenMat * (int_measure .* source)
 end
 
 "Iteratively solve meson equations and return the wave functions u(r)=[Φ0(r), W0(r), B0(r), A0(r)] where
+    divs is the number of mesh divisions so the arrays are of length (1+divs),
     r is the radius in fm,
     the inital solutions are read from s and the final solutions are saved in-place.
 Reference: P. Giuliani, K. Godbey, E. Bonilla, F. Viens, and J. Piekarewicz, Frontiers in Physics 10, (2023)"
 function solveMesonFields!(s::system, iterations=15, oscillation_control_parameter=0.3)
+    (m_s, m_ω, m_ρ, m_γ, g2_s, g2_v, g2_ρ, g2_γ, κ_ss, λ, ζ, Λv) = (s.param.m_s, s.param.m_ω, s.param.m_ρ, s.param.m_γ, s.param.g2_s, s.param.g2_v, s.param.g2_ρ, s.param.g2_γ, s.param.κ_ss, s.param.λ, s.param.ζ, s.param.Λv)
     (Φ0, W0, B0, A0) = (s.Φ0, s.W0, s.B0, s.A0)
     (ρ_sp, ρ_vp, ρ_sn, ρ_vn) = (s.ρ_sp, s.ρ_vp, s.ρ_sn, s.ρ_vn)
 
@@ -67,14 +57,16 @@ end
 
 "Calculate the total energy associated with meson fields"
 function meson_E(s::system)
-    int = 0.0
+    (κ_ss, λ, ζ, Λv) = (s.param.κ_ss, s.param.λ, s.param.ζ, s.param.Λv)
-    for (r, w, Φ0, W0, B0, A0, ρ_sp, ρ_vp, ρ_sn, ρ_vn) in zip(s.r_mesh.r, s.r_mesh.w, s.Φ0, s.W0, s.B0, s.A0, s.ρ_sp, s.ρ_vp, s.ρ_sn, s.ρ_vn)
+
+    E_densities = map(zip(s.Φ0, s.W0, s.B0, s.A0, s.ρ_sp, s.ρ_vp, s.ρ_sn, s.ρ_vn)) do (Φ0, W0, B0, A0, ρ_sp, ρ_vp, ρ_sn, ρ_vn)
         E_σ = (1/2) * (Φ0/ħc) * (ρ_sp + ρ_sn) - ((κ_ss/ħc)/12 * (Φ0/ħc)^3 + (λ/24) * (Φ0/ħc)^4)
         E_ω = -(1/2) * (W0/ħc) * (ρ_vp + ρ_vn) + (ζ/24) * (W0/ħc)^4
         E_ρ = -(1/4) * (2B0/ħc) * (ρ_vp - ρ_vn)
         E_γ = -(1/2) * (A0/ħc) * ρ_vp
         E_ωρ = Λv * (W0/ħc)^2 * (2B0/ħc)^2
-        int += (E_σ + E_ω + E_ρ + E_γ + E_ωρ) * r^2 * w
+        E_σ + E_ω + E_ρ + E_γ + E_ωρ
     end
-    return 4π * int * ħc
+    
+    return simpsons_integrate(E_densities .* rs(s).^2, Δr(s); coefficient=4π * ħc)
 end

nucleons.jl

@@ -33,8 +33,8 @@ function optimized_dirac_potentials(p::Bool, s::system)
     @. f1s =  M - s.Φ0 - s.W0 - (p - 0.5) * 2 * s.B0 - p * s.A0
     @. f2s = -M + s.Φ0 - s.W0 - (p - 0.5) * 2 * s.B0 - p * s.A0
 
-    f1 = linear_interpolation(s.r_mesh.r, f1s, extrapolation_bc = Flat())
+    f1 = linear_interpolation(rs(s), f1s)
-    f2 = linear_interpolation(s.r_mesh.r, f2s, extrapolation_bc = Flat())
+    f2 = linear_interpolation(rs(s), f2s)
 
     return (f1, f2)
 end
@@ -55,17 +55,16 @@ end
 init_BC() = [1.0, 1.0] # TODO: Why not [0.0, 1.0]?
 
 "Solve the Dirac equation and return the wave function u(r)=[g(r), f(r)] where
-    the solution would be returned as a 2×mesh_size matrix,
+    divs is the number of mesh divisions so solution would be returned as a 2×(1+divs) matrix,
-    shooting method divides the interval into two partitions at r_max/2, ensuring convergence at both r=0 and r=r_max,
     the other parameters are the same from dirac!(...)."
 function solveNucleonWf(κ, p::Bool, E, s::system; normalize=true, algo=Vern9())
     (f1, f2) = optimized_dirac_potentials(p, s)
 
     # partitioning
-    mid_idx = length(s.r_mesh.r) ÷ 2
+    mid_idx = s.divs ÷ 2
-    r_mid = s.r_mesh.r[mid_idx]
+    r_mid = rs(s)[mid_idx]
-    left_r = s.r_mesh.r[1:mid_idx]
+    left_r = rs(s)[1:mid_idx]
-    right_r = s.r_mesh.r[mid_idx:end]
+    right_r = rs(s)[mid_idx:end]
 
     # left partition
     prob = ODEProblem(dirac!, init_BC(), (0, r_mid))
@@ -73,7 +72,7 @@ function solveNucleonWf(κ, p::Bool, E, s::system; normalize=true, algo=Vern9())
     wf_left = hcat(sol.u...)
 
     # right partition
-    prob = ODEProblem(dirac!, asymp_BC(κ, p, E, s.r_mesh.r_max), (s.r_mesh.r_max, r_mid))
+    prob = ODEProblem(dirac!, asymp_BC(κ, p, E, s.r_max), (s.r_max, r_mid))
     sol = solve(prob, algo, p=(κ, E, f1, f2), saveat=right_r)
     wf_right = reverse(hcat(sol.u...); dims=2)    
     
@@ -90,7 +89,9 @@ function solveNucleonWf(κ, p::Bool, E, s::system; normalize=true, algo=Vern9())
     wf = hcat(wf_left[:, 1:(end - 1)], wf_right)
 
     if normalize
-        wf ./= sqrt(sum(wf .* wf .* transpose(s.r_mesh.w)))
+        g2_int = simpsons_integrate(wf[1, :] .^ 2, Δr(s))
+        f2_int = simpsons_integrate(wf[2, :] .^ 2, Δr(s))
+        wf ./= sqrt(g2_int + f2_int)
     end
 
     return wf
@@ -99,11 +100,11 @@ end
 "Returns a function that solves the Dirac equation in two partitions and returns 
     the determinant of [g_left(r) g_right(r); f_left(r) f_right(r)],
     where is r is in fm."
-function determinantFunc(κ, p::Bool, s::system, r::Float64=s.r_mesh.r_max/2, algo=Vern9())
+function determinantFunc(κ, p::Bool, s::system, r::Float64=s.r_max/2, algo=Vern9())
     (f1, f2) = optimized_dirac_potentials(p, s)
     prob_left  = ODEProblem(dirac!, init_BC(), (0, r))
     function func(E)
-        prob_right = ODEProblem(dirac!, asymp_BC(κ, p, E, s.r_mesh.r_max), (s.r_mesh.r_max, r))
+        prob_right = ODEProblem(dirac!, asymp_BC(κ, p, E, s.r_max), (s.r_max, r))
         u_left  = solve(prob_left,  algo, p=(κ, E, f1, f2), saveat=[r])
         u_right = solve(prob_right, algo, p=(κ, E, f1, f2), saveat=[r])
         return u_left[1, 1] * u_right[2, 1] - u_right[1, 1] * u_left[2, 1]
@@ -160,14 +161,14 @@ j_κ(κ::Int) = abs(κ) - 1/2
 "Orbital angular momentum l for a given κ value"
 l_κ(κ::Int) = abs(κ) - (κ < 0) # since true = 1 and false = 0
 
-"Calculate scalar and vector densities of a nucleon species evaluated at mesh points and returns them as vectors (ρ_s, ρ_v) where
+"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)::Tuple{Vector{Float64}, Vector{Float64}}
     spectrum = p ? s.p_spectrum : s.n_spectrum
     (κs, Es, occs) = (spectrum.κ, spectrum.E, spectrum.occ)
     
-    ρr2 = zeros(2, length(s.r_mesh.r)) # ρ×r² values
+    ρr2 = zeros(2, s.divs + 1) # ρ×r² values
     
     for (κ, E, occ) in zip(κs, Es, occs)
         wf = solveNucleonWf(κ, p, E, s; normalize=true)
@@ -175,11 +176,9 @@ function calculateNucleonDensity(p::Bool, s::system)::Tuple{Vector{Float64}, Vec
         ρr2 += (occ / (4 * pi)) * wf2 # 2j+1 factor is accounted in the occupancy number
     end
 
-    r2s = s.r_mesh.r .^ 2
+    r2s = rs(s).^2
     ρ = ρr2 ./ transpose(r2s)
-    if !all(isfinite.(ρ[:, 1]))
+    ρ[:, 1] .= ρ[:, 2] # dirty fix for NaN at r=0
-        ρ[:, 1] .= ρ[:, 2] # dirty fix for NaN at r=0
-    end
     
     ρ_s = ρ[1, :] - ρ[2, :] # g^2 - f^2
     ρ_v = ρ[1, :] + ρ[2, :] # g^2 + f^2

system.jl

@@ -1,4 +1,29 @@
-using FastGaussQuadrature
+"Stores a set of coupling constants and meson masses that goes into the Lagrangian"
+struct parameters
+    # Values defined in C. J. Horowitz and J. Piekarewicz, Phys. Rev. Lett. 86, 5647 (2001)
+    m_s::Float64 # MeV/c2
+    m_ω::Float64 # MeV/c2
+    m_ρ::Float64 # MeV/c2
+    m_γ::Float64 # MeV/c2
+    g2_s::Float64 # dimensionless
+    g2_v::Float64 # dimensionless
+    g2_ρ::Float64 # dimensionless
+    g2_γ::Float64 # dimensionless
+    κ_ss::Float64 # MeV
+    λ::Float64 # dimensionless, aka LambdaSS
+    ζ::Float64 # dimensionless, aka LambdaVV
+    Λv::Float64 # dimensionless, aka LambdaVR
+
+    "Dummy struct when parameters are not needed (for testing)"
+    parameters() = new(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)
+
+    "Initialize parameters from a string with values provided in order of struct definition separated by commas"
+    function parameters(s::String)
+        values = parse.(Float64, strip.(split(s, ',')))
+        @assert length(values) == 12 "String must contain exactly 12 values separated by commas"
+        return new(values...)
+    end
+end
 
 "Tabulates a nucleon spectrum (protons or neutrons) containing κ and occupancy"
 struct spectrum
@@ -10,35 +35,14 @@ end
 "Initializes an unfilled spectrum"
 unfilled_spectrum() = spectrum(Int[], Float64[], Int[])
 
-"Defines a mesh"
-struct mesh
-    r_max::Float64
-    r::Vector{Float64}
-    w::Vector{Float64}
-end
-
-"Create a uniform mesh"
-function uniform_mesh(r_max::Float64, divs::Int)::mesh
-    r = range(0, r_max, length=divs+1) |> collect
-    w = fill(r_max / divs, divs+1)
-    return mesh(r_max, r, w)
-end
-
-"Create a Gauss-Legendre mesh"
-function gauleg_mesh(r_max::Float64, n::Int)::mesh
-    a, b = (0.0, r_max)
-    nodes, weights = gausslegendre(n)
-    transformed_nodes = (b - a) / 2 .* nodes .+ (b + a) / 2
-    transformed_weights = (b - a) / 2 .* weights
-    return mesh(r_max, transformed_nodes, transformed_weights)
-end
-
 "Defines a nuclear system containing relevant parameters and meson/nucleon densities"
 mutable struct system
     Z::Int
     N::Int
 
-    r_mesh::mesh
+    param::parameters
+    r_max::Float64
+    divs::Int
 
     p_spectrum::spectrum
     n_spectrum::spectrum
@@ -53,21 +57,24 @@ mutable struct system
     ρ_sn::Vector{Float64}
     ρ_vn::Vector{Float64}
 
-    "Initialize an unsolved system with a uniform r-mesh"
+    "Initialize an unsolved system"
-    system(Z::Int, N::Int, r_max::Float64, divs::Int) = new(Z, N, uniform_mesh(r_max, divs), unfilled_spectrum(), unfilled_spectrum(), [zeros(1 + divs) for _ in 1:8]...)
+    system(Z, N, parameters, r_max, divs) = new(Z, N, parameters, r_max, divs, unfilled_spectrum(), unfilled_spectrum(), [zeros(1 + divs) for _ in 1:8]...)
 
-    "Dummy struct to define the mesh"
+    "Dummy struct to define the mesh (for testing)"
-    system(r_max, divs) = system(0, 0, r_max, divs)
+    system(r_max, divs) = system(0, 0, parameters(), r_max, divs)
-
-    "Initialize an unsolved system with a Gauss-Legendre r-mesh"
-    system(Z::Int, N::Int, mesh_size::Int, r_max::Float64) = new(Z, N, gauleg_mesh(r_max, mesh_size), unfilled_spectrum(), unfilled_spectrum(), [zeros(mesh_size) for _ in 1:8]...)
 end
 
 "Get mass number of nucleus"
 A(s::system)::Int = s.Z + s.N
 
+"Get r values in the mesh"
+rs(s::system)::StepRangeLen = range(0, s.r_max, length=s.divs+1)
+
+"Get Δr value for the mesh"
+Δr(s::system)::Float64 = s.r_max / s.divs
+
 "Get the number of protons or neutrons in the system"
 Z_or_N(s::system, p::Bool)::Int = p ? s.Z : s.N
 
 "Create an empty array for the size of the mesh"
-zero_array(s::system) = zeros(length(s.r_mesh.r))
+zero_array(s::system) = zeros(1 + s.divs)

test/Linear_nucleon_spectrum.jl

@@ -13,7 +13,7 @@ As = test_data[:, 5]
 p = false
 r_max = maximum(xs)
 divs = length(xs) - 1
-s = system(8, 8, r_max, divs)
+s = system(8, 8, parameters(), r_max, divs)
 
 s.Φ0 = Ss
 s.W0 = Vs

test/Pb208.jl

@@ -2,9 +2,13 @@ include("../NuclearRMF.jl")
 
 Z = 82
 N = 126
+
+# Parameter values calibrated by FSUGarnet for Pb-208
+param = parameters("496.939473213388, 782.5, 763.0, 0.000001000, 110.349189097820, 187.694676506801, 192.927428365698, 0.091701236, 3.260178893447, -0.003551486718, 0.023499504053, 0.043376933644")
+
 r_max = 20.0
-mesh_points = 400
+divs = 400
 
-s = system(Z, N, mesh_points, r_max)
+s = system(Z, N, param, r_max, divs)
 
-solve_system!(s; monitor_plot=true)
+solve_system!(s; live_plots=false)

test/Pb208ParamsFSUGarnet.txt

@@ -0,0 +1 @@
+496.939473213388, 782.5, 763.0, 0.000001000, 110.349189097820, 187.694676506801, 192.927428365698, 0.091701236, 3.260178893447, -0.003551486718, 0.023499504053, 0.043376933644

test/Pb208_meson_flds.jl

@@ -17,9 +17,12 @@ plot(xs_bench, hcat(Φ0_bench, W0_bench, B0_bench, A0_bench), layout=4, label=["
 test_data = readdlm("test/Pb208DensFSUGarnet.csv")
 xs = test_data[:, 1]
 
+N_p = 82
+N_n = 126
+param = parameters(read("test/Pb208ParamsFSUGarnet.txt", String))
 r_max = maximum(xs)
 divs = length(xs) - 1
-s = system(r_max, divs)
+s = system(N_p, N_n, param, r_max, divs)
 
 s.ρ_sn = test_data[:, 2]
 s.ρ_vn = test_data[:, 3]

test/Pb208_nucleon_dens.jl

@@ -14,7 +14,7 @@ N_p = 82
 N_n = 126
 r_max = maximum(xs)
 divs = length(xs) - 1
-s = system(N_p, N_n, r_max, divs)
+s = system(N_p, N_n, parameters(), r_max, divs)
 
 s.Φ0 = Ss
 s.W0 = Vs
@@ -23,13 +23,13 @@ s.A0 = As
 
 solveNucleonDensity!(true, s)
 
-p_sp = plot(s.r_mesh.r, s.ρ_sp, xlabel="r (fm)", label="ρₛₚ(r) calculated")
+p_sp = plot(rs(s), s.ρ_sp, xlabel="r (fm)", label="ρₛₚ(r) calculated")
-p_vp = plot(s.r_mesh.r, s.ρ_vp, xlabel="r (fm)", label="ρᵥₚ(r) calculated")
+p_vp = plot(rs(s), s.ρ_vp, xlabel="r (fm)", label="ρᵥₚ(r) calculated")
 
 solveNucleonDensity!(false, s)
 
-p_sn = plot(s.r_mesh.r, s.ρ_sn, xlabel="r (fm)", label="ρₛₙ(r) calculated")
+p_sn = plot(rs(s), s.ρ_sn, xlabel="r (fm)", label="ρₛₙ(r) calculated")
-p_vn = plot(s.r_mesh.r, s.ρ_vn, xlabel="r (fm)", label="ρᵥₙ(r) calculated")
+p_vn = plot(rs(s), s.ρ_vn, xlabel="r (fm)", label="ρᵥₙ(r) calculated")
 
 # benchmark data generated from Hartree.f
 # format: x  Rhos(n)  Rhov(n)  Rhot(n)  Rhos(p)  Rhov(p)  Rhot(p)

test/Pb208_nucleon_spectrum.jl

@@ -15,7 +15,7 @@ N_p = 82
 N_n = 126
 r_max = maximum(xs)
 divs = length(xs) - 1
-s = system(N_p, N_n, r_max, divs)
+s = system(N_p, N_n, parameters(), r_max, divs)
 
 s.Φ0 = Ss
 s.W0 = Vs

test/Pb208_nucleon_wf.jl

@@ -28,6 +28,6 @@ wf = solveNucleonWf(κ, p, groundE, s)
 gs = wf[1, :]
 fs = wf[2, :]
 
-plot(s.r_mesh.r, gs, label="g(r)")
+plot(rs(s), gs, label="g(r)")
-plot!(s.r_mesh.r, fs, label="f(r)")
+plot!(rs(s), fs, label="f(r)")
 xlabel!("r (fm)")