Avoid dirty trick if unnecessary

NuclearRMF.jl

@@ -11,7 +11,7 @@ Woods_Saxon(r::Float64; R::Float64=7.0, a::Float64=0.5) = -1 / (8π * a^3 * reli
 "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)
     if reinitialize_densities
-        dens_guess = Woods_Saxon.(rs(s))
+        dens_guess = Woods_Saxon.(s.r_mesh.r)
         @. s.ρ_sp = s.Z * dens_guess
         @. s.ρ_vp = s.Z * dens_guess
         @. s.ρ_sn = s.N * dens_guess
@@ -33,7 +33,7 @@ function solve_system!(s::system; reinitialize_densities=true, monitor_print=tru
             for s in p.series_list
                 s.plotattributes[:linecolor] = :gray
             end
-            plot!(p, rs(s), hcat(s.ρ_sp, s.ρ_vp, s.ρ_sn, s.ρ_vn, s.Φ0, s.W0, s.B0, s.A0), linecolor=:red)
+            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)
             display(p)
         end
 

Project.toml

@@ -1,5 +1,6 @@
 [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"

mesons.jl

@@ -25,16 +25,15 @@ 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 .* Δr(s) .* rs(s) .^ 2
+    int_measure = ħc .* s.r_mesh.w .* (s.r_mesh.r .^ 2)
 
     greensFunction = m == 0 ? greensFunctionP : (r, rp) -> greensFunctionKG(m / ħc, r, rp)
-    greenMat = greensFunction.(rs(s), transpose(rs(s)))
+    greenMat = greensFunction.(s.r_mesh.r, transpose(s.r_mesh.r))
 
     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)"
@@ -69,13 +68,13 @@ end
 "Calculate the total energy associated with meson fields"
 function meson_E(s::system)
     int = 0.0
-    for (r, Φ0, W0, B0, A0, ρ_sp, ρ_vp, ρ_sn, ρ_vn) in zip(rs(s), s.Φ0, s.W0, s.B0, s.A0, s.ρ_sp, s.ρ_vp, s.ρ_sn, s.ρ_vn)
+    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_σ = (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
+        int += (E_σ + E_ω + E_ρ + E_γ + E_ωρ) * r^2 * w
     end
-    return 4π * int * Δr(s) * ħc
+    return 4π * int * ħ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(rs(s), f1s)
+    f1 = linear_interpolation(s.r_mesh.r, f1s, extrapolation_bc = Flat())
-    f2 = linear_interpolation(rs(s), f2s)
+    f2 = linear_interpolation(s.r_mesh.r, f2s, extrapolation_bc = Flat())
 
     return (f1, f2)
 end
@@ -55,16 +55,17 @@ 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
-    divs is the number of mesh divisions so solution would be returned as a 2×(1+divs) matrix,
+    the solution would be returned as a 2×mesh_size 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 = s.divs ÷ 2
+    mid_idx = length(s.r_mesh.r) ÷ 2
-    r_mid = rs(s)[mid_idx]
+    r_mid = s.r_mesh.r[mid_idx]
-    left_r = rs(s)[1:mid_idx]
+    left_r = s.r_mesh.r[1:mid_idx]
-    right_r = rs(s)[mid_idx:end]
+    right_r = s.r_mesh.r[mid_idx:end]
 
     # left partition
     prob = ODEProblem(dirac!, init_BC(), (0, r_mid))
@@ -72,7 +73,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_max), (s.r_max, r_mid))
+    prob = ODEProblem(dirac!, asymp_BC(κ, p, E, s.r_mesh.r_max), (s.r_mesh.r_max, r_mid))
     sol = solve(prob, algo, p=(κ, E, f1, f2), saveat=right_r)
     wf_right = reverse(hcat(sol.u...); dims=2)
     
@@ -89,7 +90,7 @@ function solveNucleonWf(κ, p::Bool, E, s::system; normalize=true, algo=Vern9())
     wf = hcat(wf_left[:, 1:(end - 1)], wf_right)
 
     if normalize
-        wf ./= norm(wf) * sqrt(Δr(s)) # integration by Reimann sum
+        wf ./= sqrt(sum(wf .* wf .* transpose(s.r_mesh.w)))
     end
 
     return wf
@@ -98,11 +99,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_max/2, algo=Vern9())
+function determinantFunc(κ, p::Bool, s::system, r::Float64=s.r_mesh.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_max), (s.r_max, r))
+        prob_right = ODEProblem(dirac!, asymp_BC(κ, p, E, s.r_mesh.r_max), (s.r_mesh.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]
@@ -159,14 +160,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 on [0,r_max] divided into (divs+1) points and returns them as vectors (ρ_s, ρ_v) where
+"Calculate scalar and vector densities of a nucleon species evaluated at mesh 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, s.divs + 1) # ρ×r² values
+    ρr2 = zeros(2, length(s.r_mesh.r)) # ρ×r² values
     
     for (κ, E, occ) in zip(κs, Es, occs)
         wf = solveNucleonWf(κ, p, E, s; normalize=true)
@@ -174,9 +175,11 @@ 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 = rs(s).^2
+    r2s = s.r_mesh.r .^ 2
     ρ = ρr2 ./ transpose(r2s)
-    ρ[:, 1] .= ρ[:, 2] # dirty fix for NaN at r=0
+    if !all(isfinite.(ρ[:, 1]))
+        ρ[:, 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,3 +1,5 @@
+using FastGaussQuadrature
+
 "Tabulates a nucleon spectrum (protons or neutrons) containing κ and occupancy"
 struct spectrum
     κ::Vector{Int}
@@ -8,13 +10,35 @@ 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_max::Float64
+    r_mesh::mesh
-    divs::Int
 
     p_spectrum::spectrum
     n_spectrum::spectrum
@@ -29,24 +53,21 @@ mutable struct system
     ρ_sn::Vector{Float64}
     ρ_vn::Vector{Float64}
 
-    "Initialize an unsolved system"
+    "Initialize an unsolved system with a uniform r-mesh"
-    system(Z, N, r_max, divs) = new(Z, N, r_max, divs, unfilled_spectrum(), unfilled_spectrum(), [zeros(1 + divs) for _ in 1:8]...)
+    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]...)
 
     "Dummy struct to define the mesh"
     system(r_max, divs) = system(0, 0, 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(1 + s.divs)
+zero_array(s::system) = zeros(length(s.r_mesh.r))

test/Pb208.jl

@@ -3,8 +3,8 @@ include("../NuclearRMF.jl")
 Z = 82
 N = 126
 r_max = 20.0
-divs = 400
+mesh_points = 400
 
-s = system(Z, N, r_max, divs)
+s = system(Z, N, mesh_points, r_max)
 
 solve_system!(s; monitor_plot=true)

test/Pb208_nucleon_dens.jl

@@ -23,13 +23,13 @@ s.A0 = As
 
 solveNucleonDensity!(true, s)
 
-p_sp = plot(rs(s), s.ρ_sp, xlabel="r (fm)", label="ρₛₚ(r) calculated")
+p_sp = plot(s.r_mesh.r, s.ρ_sp, xlabel="r (fm)", label="ρₛₚ(r) calculated")
-p_vp = plot(rs(s), s.ρ_vp, xlabel="r (fm)", label="ρᵥₚ(r) calculated")
+p_vp = plot(s.r_mesh.r, s.ρ_vp, xlabel="r (fm)", label="ρᵥₚ(r) calculated")
 
 solveNucleonDensity!(false, s)
 
-p_sn = plot(rs(s), s.ρ_sn, xlabel="r (fm)", label="ρₛₙ(r) calculated")
+p_sn = plot(s.r_mesh.r, s.ρ_sn, xlabel="r (fm)", label="ρₛₙ(r) calculated")
-p_vn = plot(rs(s), s.ρ_vn, xlabel="r (fm)", label="ρᵥₙ(r) calculated")
+p_vn = plot(s.r_mesh.r, 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_wf.jl

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