Merge branch 'structs'

mesons.jl

@@ -12,7 +12,7 @@ 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 κ = 3.260178893447 # MeV
+const κ_ss = 3.260178893447 # MeV
 const λ = -0.003551486718 # dimensionless # LambdaSS
 const ζ = 0.023499504053 # dimensionless # LambdaVV
 const Λv = 0.043376933644 # dimensionless # LambdaVR
@@ -26,42 +26,38 @@ greensFunctionKG(m, r, rp) = -1 / (m * (r + r_reg) * (rp + r_reg)) * sinh(m * mi
 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,
+    m is the mass of the meson in MeV/c2."
-    r_max is the r-cutoff in fm."
+function solveKG(m, source, s::system)
-function solveKG(m, source, r_max)
+    int_measure = ħc .* Δr(s) .* rs(s) .^ 2
-    Δr = r_max / (length(source) - 1)
-    rs = range(0, r_max; length=length(source))
-    int_measure = ħc .* Δr .* rs .^ 2
 
     greensFunction = m == 0 ? greensFunctionP : (r, rp) -> greensFunctionKG(m / ħc, r, rp)
-    greenMat = greensFunction.(rs, transpose(rs))
+    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),
-    ρ_sp, ρ_vp are the scalar and vector proton densities as arrays,
-    ρ_sn, ρ_vn are the scalar and vector neutron densities as arrays,
-    Φ0, W0, B0, A0 are the mean-field potentials (couplings included) in MeV as as arrays,
     r is the radius in fm,
-    An initial guess initial_sol=(Φ0, W0, B0, A0) can be provided to speed up convergence (permuting!).
+    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 solveMesonWfs(ρ_sp, ρ_vp, ρ_sn, ρ_vn, r_max, divs, iterations=10; initial_sol=(zeros(1 + divs) for _ in 1:4))
+function solveMesonFields!(s::system, iterations=10)
-    (Φ0, W0, B0, A0) = initial_sol
+    (Φ0, W0, B0, A0) = (s.Φ0, s.W0, s.B0, s.A0)
-    (src_Φ0, src_W0, src_B0, src_A0) = (zeros(1 + divs) for _ in 1:4)
+    (ρ_sp, ρ_vp, ρ_sn, ρ_vn) = (s.ρ_sp, s.ρ_vp, s.ρ_sn, s.ρ_vn)
+
+    (src_Φ0, src_W0, src_B0, src_A0) = (zero_array(s) for _ in 1:4)
 
     # A0 doesn't need iterations
     @. src_A0 = -g2_γ * ρ_vp
-    A0 .= solveKG(m_γ, src_A0, r_max)
+    A0 .= solveKG(m_γ, src_A0, s)
 
     for _ in 1:iterations
-        @. src_Φ0 = g2_s * ((κ/ħc)/2 * (Φ0/ħc)^2 + (λ/6) * (Φ0/ħc)^3) - g2_s * (ρ_sp + ρ_sn)
+        @. src_Φ0 = g2_s * ((κ_ss/ħc)/2 * (Φ0/ħc)^2 + (λ/6) * (Φ0/ħc)^3) - g2_s * (ρ_sp + ρ_sn)
         @. src_W0 = g2_v * ((ζ/6) * (W0/ħc)^3 + 2Λv * (2B0/ħc)^2 * (W0/ħc)) - g2_v * (ρ_vp + ρ_vn)
         @. src_B0 = (2Λv * g2_ρ * (W0/ħc)^2 * (2B0/ħc) - g2_ρ/2 * (ρ_vp - ρ_vn)) / 2
-        Φ0 .= solveKG(m_s, src_Φ0, r_max)
+        Φ0 .= solveKG(m_s, src_Φ0, s)
-        W0 .= (solveKG(m_ω, src_W0, r_max) .+ W0) ./ 2 # to supress oscillation
+        W0 .= (solveKG(m_ω, src_W0, s) .+ W0) ./ 2 # to supress oscillation
-        B0 .= (solveKG(m_ρ, src_B0, r_max) .+ B0) ./ 2 # to supress oscillation
+        B0 .= (solveKG(m_ρ, src_B0, s) .+ B0) ./ 2 # to supress oscillation
     end
 
     return (Φ0, W0, B0, A0)

nucleons.jl

@@ -27,19 +27,21 @@ end
     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, E, Φ0, W0, B0, A0, r_max, divs; shooting=true, normalize=true, algo=Tsit5())
+function solveNucleonWf(κ, p, E, s::system; shooting=true, normalize=true, algo=Tsit5())
-    Δr = r_max / divs
+    (Φ0, W0, B0, A0) = fields_interp(s)
 
     if shooting
-        @assert divs % 2 == 0 "divs must be an even number when shooting=true"
+        @assert s.divs % 2 == 0 "divs must be an even number when shooting=true"
-        prob = ODEProblem(dirac!, [0, 1], (r_max, r_max / 2))
+        prob = ODEProblem(dirac!, [0, 1], (s.r_max, s.r_max / 2))
-        sol = solve(prob, algo, p=(κ, p, E, Φ0, W0, B0, A0), saveat=Δr)
+        sol = solve(prob, algo, p=(κ, p, E, Φ0, W0, B0, A0), saveat=Δr(s))
         wf_right = reverse(hcat(sol.u...); dims=2)
-        r_max = r_max / 2 # for the next step
+        next_r_max = s.r_max / 2 # for the next step
+    else
+        next_r_max = s.r_max
     end
     
-    prob = ODEProblem(dirac!, [0, 1], (0, r_max))
+    prob = ODEProblem(dirac!, [0, 1], (0, next_r_max))
-    sol = solve(prob, algo, p=(κ, p, E, Φ0, W0, B0, A0), saveat=Δr)
+    sol = solve(prob, algo, p=(κ, p, E, Φ0, W0, B0, A0), saveat=Δr(s))
     wf = hcat(sol.u...)
 
     if shooting # join two segments
@@ -52,7 +54,7 @@ function solveNucleonWf(κ, p, E, Φ0, W0, B0, A0, r_max, divs; shooting=true, n
     end
 
     if normalize
-        norm = sum(wf .* wf) * Δr # integration by Reimann sum
+        norm = sum(wf .* wf) * Δr(s) # integration by Reimann sum
         wf = wf ./ sqrt(norm)
     end
 
@@ -62,30 +64,31 @@ end
 "Returns a function that solves the Dirac equation and returns g(r=r_max) where
     r_max is the outer boundary in fm,
     the other parameters are the same from dirac!(...)."
-function boundaryValueFunc(κ, p, Φ0, W0, B0, A0, r_max; dtype=Float64, algo=Tsit5())
+function boundaryValueFunc(κ, p, s::system; dtype=Float64, algo=Tsit5())
-    prob = ODEProblem(dirac!, convert.(dtype, [0, 1]), (0, r_max))
+    prob = ODEProblem(dirac!, convert.(dtype, [0, 1]), (0, s.r_max))
-    func(E) = solve(prob, algo, p=(κ, p, E, Φ0, W0, B0, A0), saveat=[r_max], save_idxs=[1])[1, 1]
+    (Φ0, W0, B0, A0) = fields_interp(s)
+    func(E) = solve(prob, algo, p=(κ, p, E, Φ0, W0, B0, A0), saveat=[s.r_max], save_idxs=[1])[1, 1]
     return func
 end
 
 "Find all bound energies between E_min (=0) and E_max (=mass) where
     tol_digits determines the precision for root finding and the threshold for identifying duplicates,
     the other parameters are the same from dirac!(...)."
-function findEs(κ, p, Φ0, W0, B0, A0, r_max, E_min=0, E_max=(p ? M_p : M_n), tol_digits=5)
+function findEs(κ, p, s::system, E_min=0, E_max=(p ? M_p : M_n), tol_digits=5)
-    func = boundaryValueFunc(κ, p, Φ0, W0, B0, A0, r_max)
+    func = boundaryValueFunc(κ, p, s)
     Es = find_zeros(func, (E_min, E_max); xatol=1/10^tol_digits)
     return unique(E -> round(E; digits=tol_digits), Es)
 end
 
 "Find all orbitals and return two lists containing κ values and corresponding energies for a single species where
     the other parameters are defined above"
-function findAllOrbitals(p, Φ0, W0, B0, A0, r_max, E_min=0, E_max=(p ? M_p : M_n))
+function findAllOrbitals(p, s::system, E_min=0, E_max=(p ? M_p : M_n))
     κs = Int[]
     Es = Float64[]
     # start from κ=-1 and go both up and down
     for direction in [-1, 1]
         for κ in direction * (1:100) # cutoff is 100
-            new_Es = findEs(κ, p, Φ0, W0, B0, A0, r_max, E_min, E_max)
+            new_Es = findEs(κ, p, s, E_min, E_max)
             if isempty(new_Es); break; end
             append!(Es, new_Es)
             append!(κs, fill(κ, length(new_Es)))
@@ -94,20 +97,20 @@ function findAllOrbitals(p, Φ0, W0, B0, A0, r_max, E_min=0, E_max=(p ? M_p : M_
     return (κs, Es)
 end
 
-"For a given list of κ values with corresponding energies, attempt to fill N lowest lying orbitals and return occupancy numbers"
+"For a given list of κ values with corresponding energies, attempt to fill Z_or_N lowest lying orbitals and return occupancy numbers"
-function fillNucleons(N::Int, κs, Es)
+function fillNucleons(Z_or_N::Int, κs, Es)
     sort_i = sortperm(Es)
 
     occ = zeros(Int, length(κs))
 
     for i in sort_i
-        if N ≤ 0; break; end;
+        if Z_or_N ≤ 0; break; end;
         max_occ = 2 * j_κ(κs[i]) + 1
-        occ[i] = min(max_occ, N)
+        occ[i] = min(max_occ, Z_or_N)
-        N -= occ[i]
+        Z_or_N -= occ[i]
     end
 
-    N == 0 || @warn "All orbitals could not be filled"
+    Z_or_N == 0 || @warn "All orbitals could not be filled"
     return occ
 end
 
@@ -120,16 +123,16 @@ 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, Φ0, W0, B0, A0, r_max, divs)
+function calculateNucleonDensity(κs, Es, occs, p, s::system)
-    ρr2 = zeros(2, divs + 1) # ρ×r² values
+    ρr2 = zeros(2, s.divs + 1) # ρ×r² values
     
     for (κ, E, occ) in zip(κs, Es, occs)
-        wf = solveNucleonWf(κ, p, E, Φ0, W0, B0, A0, r_max, divs; shooting=true, normalize=true)
+        wf = solveNucleonWf(κ, p, E, s; shooting=true, normalize=true)
         wf2 = wf .* wf
         ρr2 += (occ / (4 * pi)) * wf2 # 2j+1 factor is accounted in the occupancy number
     end
 
-    r2s = (collect ∘ range)(0, r_max, length=divs+1).^2
+    r2s = rs(s).^2
     ρ = ρr2 ./ transpose(r2s)
     ρ[:, 1] .= ρ[:, 2] # dirty fix for NaN at r=0
     
@@ -141,8 +144,8 @@ end
 
 "Solve the Dirac equation and calculate scalar and vector densities of a nucleon species where
     the other parameters are defined above"
-function solveNucleonDensity(N, p, Φ0, W0, B0, A0, r_max, divs, E_min=800, E_max=939)
+function solveNucleonDensity(p, s::system, E_min=800, E_max=939)
-    κs, Es = findAllOrbitals(p, Φ0, W0, B0, A0, r_max, E_min, E_max)
+    κs, Es = findAllOrbitals(p, s, E_min, E_max)
-    occs = fillNucleons(N, κs, Es)
+    occs = fillNucleons(Z_or_N(s, p), κs, Es)
-    return calculateNucleonDensity(κs, Es, occs, p, Φ0, W0, B0, A0, r_max, divs)
+    return calculateNucleonDensity(κs, Es, occs, p, s)
 end

system.jl

@@ -1,23 +1,41 @@
 using Interpolations, PolyLog, Plots
-include("nucleons.jl")
-include("mesons.jl")
 
-"Defines a nuclear system to be solved"
+"Defines a nuclear system containing relevant parameters and meson/nucleon densities"
-struct system
+mutable struct system
     Z::Int
     N::Int
+
     r_max::Float64
     divs::Int
+
+    Φ0::Vector{Float64}
+    W0::Vector{Float64}
+    B0::Vector{Float64}
+    A0::Vector{Float64}
+
+    ρ_sp::Vector{Float64}
+    ρ_vp::Vector{Float64}
+    ρ_sn::Vector{Float64}
+    ρ_vn::Vector{Float64}
+
+    "Initialize an unsolved system"
+    system(Z, N, r_max, divs) = new(Z, N, r_max, divs, [zeros(1 + divs) for _ in 1:8]...)
+
+    "Dummy struct to define the mesh"
+    system(r_max, divs) = system(0, 0, r_max, divs)
 end
 
 "Get mass number of nucleus"
-A(s::system) = s.Z + s.N
+A(s::system)::Int = s.Z + s.N
 
 "Get r values in the mesh"
-rs(s::system) = range(0, s.r_max, length=s.divs+1)
+rs(s::system)::StepRangeLen = range(0, s.r_max, length=s.divs+1)
 
 "Get Δr value for the mesh"
-Δr(s::system) = s.r_max / s.divs
+Δ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)
@@ -25,19 +43,27 @@ zero_array(s::system) = zeros(1 + s.divs)
 "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)))
 
-"Run the full program by self-consistent solution of nucleon and meson densities"
+"Create linear interpolations for the meson fields"
-function solve_system(s::system, initial_dens=nothing, initial_flds=(zeros(1 + s.divs) for _ in 1:4); monitor_print=true, monitor_plot=false)
+function fields_interp(s::system)
-    if isnothing(initial_dens)
+    S_interp = linear_interpolation(rs(s), s.Φ0)
-        dens_guess = Woods_Saxon.(rs(s))
+    V_interp = linear_interpolation(rs(s), s.W0)
-        ρ_sp = s.Z .* dens_guess
+    R_interp = linear_interpolation(rs(s), s.B0)
-        ρ_vp = s.Z .* dens_guess
+    A_interp = linear_interpolation(rs(s), s.A0)
-        ρ_sn = s.N .* dens_guess
+    return (S_interp, V_interp, R_interp, A_interp)
-        ρ_vn = s.N .* dens_guess
+end
-    else
-        (ρ_sp, ρ_vp, ρ_sn, ρ_vn) = initial_dens
-    end
 
-    (Φ0s, W0s, B0s, A0s) = initial_flds
+include("nucleons.jl")
+include("mesons.jl")
+
+"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))
+        @. 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
         p = plot(legends=false, size=(1024, 768), layout=(2, 4), title=["ρₛₚ" "ρᵥₚ" "ρₛₙ" "ρᵥₙ" "Φ₀" "W₀" "B₀" "A₀"])
@@ -46,28 +72,24 @@ function solve_system(s::system, initial_dens=nothing, initial_flds=(zeros(1 + s
     E_total_previous = NaN
 
     while true
-        @time "Meson fields" (Φ0s, W0s, B0s, A0s) = solveMesonWfs(ρ_sp, ρ_vp, ρ_sn, ρ_vn, s.r_max, s.divs, isnan(E_total_previous) ? 50 : 5; initial_sol = (Φ0s, W0s, B0s, A0s))
+        # mesons
-
+        @time "Meson fields" solveMesonFields!(s, isnan(E_total_previous) ? 50 : 5)
-        S_interp = linear_interpolation(rs(s), Φ0s)
-        V_interp = linear_interpolation(rs(s), W0s)
-        R_interp = linear_interpolation(rs(s), B0s)
-        A_interp = linear_interpolation(rs(s), A0s)
 
         # protons
-        @time "Proton spectrum" (κs_p, Es_p) = findAllOrbitals(true, S_interp, V_interp, R_interp, A_interp, s.r_max)
+        @time "Proton spectrum" (κs_p, Es_p) = findAllOrbitals(true, s)
         occs_p = fillNucleons(s.Z, κs_p, Es_p)
-        @time "Proton densities" (ρ_sp, ρ_vp) = calculateNucleonDensity(κs_p, Es_p, occs_p, true, S_interp, V_interp, R_interp, A_interp, s.r_max, s.divs)
+        @time "Proton densities" (s.ρ_sp, s.ρ_vp) = calculateNucleonDensity(κs_p, Es_p, occs_p, true, s)
 
         # neutrons
-        @time "Neutron spectrum" (κs_n, Es_n) = findAllOrbitals(false, S_interp, V_interp, R_interp, A_interp, s.r_max)
+        @time "Neutron spectrum" (κs_n, Es_n) = findAllOrbitals(false, s)
         occs_n = fillNucleons(s.N, κs_n, Es_n)
-        @time "Neutron densities" (ρ_sn, ρ_vn) = calculateNucleonDensity(κs_n, Es_n, occs_n, false, S_interp, V_interp, R_interp, A_interp, s.r_max, s.divs)
+        @time "Neutron densities" (s.ρ_sn, s.ρ_vn) = calculateNucleonDensity(κs_n, Es_n, occs_n, false, s)
 
         if monitor_plot
             for s in p.series_list
                 s.plotattributes[:linecolor] = :gray
             end
-            plot!(p, rs(s), hcat(ρ_sp, ρ_vp, ρ_sn, ρ_vn, Φ0s, W0s, B0s, A0s), 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
 

test/Linear_nucleon_spectrum.jl

@@ -1,5 +1,5 @@
 using DelimitedFiles, Interpolations, Plots
-include("../nucleons.jl")
+include("../system.jl")
 
 # test data generated from Hartree.f
 # format: x S(x) V(x) R(x) A(x)
@@ -10,20 +10,22 @@ Vs = test_data[:, 3]
 Rs = test_data[:, 4]
 As = test_data[:, 5]
 
-S_interp = linear_interpolation(xs, Ss)
-V_interp = linear_interpolation(xs, Vs)
-R_interp = linear_interpolation(xs, Rs)
-A_interp = linear_interpolation(xs, As)
-
 p = false
-N = p ? 8 : 8
 r_max = maximum(xs)
+divs = length(xs) - 1
+s = system(8, 8, r_max, divs)
+
+s.Φ0 = Ss
+s.W0 = Vs
+s.B0 = Rs
+s.A0 = As
+
 E_min = 860
 E_max = 939
 
-κs, Es = findAllOrbitals(p, S_interp, V_interp, R_interp, A_interp, r_max, E_min, E_max)
+κs, Es = findAllOrbitals(p, s, E_min, E_max)
 Ebinds = (p ? M_p : M_n) .- Es
-occ = fillNucleons(N, κs, Es)
+occ = fillNucleons(Z_or_N(s, p), κs, Es)
 
 # format: proton, kappa, filling, gnodes, fnodes, Ebind
 bench_data, _ = readdlm("test/LinearSpectrum.csv", ','; header=true)

test/Pb208_meson_flds.jl

@@ -1,5 +1,5 @@
 using DelimitedFiles, Interpolations, Plots
-include("../mesons.jl")
+include("../system.jl")
 
 # test data generated from Hartree.f
 # format: x S(x) V(x) R(x) A(x)
@@ -16,15 +16,17 @@ plot(xs_bench, hcat(Φ0_bench, W0_bench, B0_bench, A0_bench), layout=4, label=["
 # format: x  Rhos(n)  Rhov(n)  Rhot(n)  Rhos(p)  Rhov(p)  Rhot(p)
 test_data = readdlm("test/Pb208DensFSUGarnet.csv")
 xs = test_data[:, 1]
-ρ_sn = test_data[:, 2]
-ρ_vn = test_data[:, 3]
-ρ_sp = test_data[:, 5]
-ρ_vp = test_data[:, 6]
 
 r_max = maximum(xs)
 divs = length(xs) - 1
+s = system(r_max, divs)
 
-(Φ0, W0, B0, A0) = solveMesonWfs(ρ_sp, ρ_vp, ρ_sn, ρ_vn, r_max, divs, 200) # many iterations needed without an initial solution
+s.ρ_sn = test_data[:, 2]
+s.ρ_vn = test_data[:, 3]
+s.ρ_sp = test_data[:, 5]
+s.ρ_vp = test_data[:, 6]
+
+(Φ0, W0, B0, A0) = solveMesonFields!(s, 200) # many iterations needed without an initial solution
 
 plot!(xs, hcat(Φ0, W0, B0, A0), layout=4, label=["Φ0" "W0" "B0" "A0"])
 xlabel!("r (fm)")

test/Pb208_nucleon_basic.jl

@@ -1,5 +1,5 @@
-using DelimitedFiles, Interpolations, Plots
+using DelimitedFiles, Plots
-include("../nucleons.jl")
+include("../system.jl")
 
 # test data generated from Hartree.f
 # format: x S(x) V(x) R(x) A(x)
@@ -10,23 +10,26 @@ Vs = test_data[:, 3]
 Rs = test_data[:, 4]
 As = test_data[:, 5]
 
-S_interp = linear_interpolation(xs, Ss)
-V_interp = linear_interpolation(xs, Vs)
-R_interp = linear_interpolation(xs, Rs)
-A_interp = linear_interpolation(xs, As)
-
 κ = -1
 p = false
 r_max = maximum(xs)
+divs = length(xs) - 1
+s = system(r_max, divs)
+
+s.Φ0 = Ss
+s.W0 = Vs
+s.B0 = Rs
+s.A0 = As
+
 E_min = 850
 E_max = 939
 
-boundEs = findEs(κ, p, S_interp, V_interp, R_interp, A_interp, r_max, E_min, E_max)
+boundEs = findEs(κ, p, s, E_min, E_max)
 
 binding_Es = (p ? M_p : M_n) .- boundEs
 println("binding energies = $binding_Es")
 
-func = boundaryValueFunc(κ, p, S_interp, V_interp, R_interp, A_interp, r_max)
+func = boundaryValueFunc(κ, p, s)
 Es = collect(E_min:0.5:E_max)
 boundaryVals = [func(E)^2 for E in Es]
 

test/Pb208_nucleon_dens.jl

@@ -1,5 +1,5 @@
-using DelimitedFiles, Interpolations, Plots
+using DelimitedFiles, Plots
-include("../nucleons.jl")
+include("../system.jl")
 
 # test data generated from Hartree.f
 # format: x S(x) V(x) R(x) A(x)
@@ -10,29 +10,26 @@ Vs = test_data[:, 3]
 Rs = test_data[:, 4]
 As = test_data[:, 5]
 
-S_interp = linear_interpolation(xs, Ss)
-V_interp = linear_interpolation(xs, Vs)
-R_interp = linear_interpolation(xs, Rs)
-A_interp = linear_interpolation(xs, As)
-
 N_p = 82
 N_n = 126
 r_max = maximum(xs)
-E_min = 860
+divs = length(xs) - 1
-E_max = 939
+s = system(N_p, N_n, r_max, divs)
-divs = 400
 
-rs = range(0, r_max, length=divs+1)
+s.Φ0 = Ss
+s.W0 = Vs
+s.B0 = Rs
+s.A0 = As
 
-(ρ_sp, ρ_vp) = solveNucleonDensity(N_p, true, S_interp, V_interp, R_interp, A_interp, r_max, divs, E_min, E_max)
+(ρ_sp, ρ_vp) = solveNucleonDensity(true, s)
 
-p_sp = plot(rs, ρ_sp, xlabel="r (fm)", label="ρₛₚ(r) calculated")
+p_sp = plot(rs(s), ρ_sp, xlabel="r (fm)", label="ρₛₚ(r) calculated")
-p_vp = plot(rs, ρ_vp, xlabel="r (fm)", label="ρᵥₚ(r) calculated")
+p_vp = plot(rs(s), ρ_vp, xlabel="r (fm)", label="ρᵥₚ(r) calculated")
 
-(ρ_sn, ρ_vn) = solveNucleonDensity(N_n, false, S_interp, V_interp, R_interp, A_interp, r_max, divs, E_min, E_max)
+(ρ_sn, ρ_vn) = solveNucleonDensity(false, s)
 
-p_sn = plot(rs, ρ_sn, xlabel="r (fm)", label="ρₛₙ(r) calculated")
+p_sn = plot(rs(s), ρ_sn, xlabel="r (fm)", label="ρₛₙ(r) calculated")
-p_vn = plot(rs, ρ_vn, xlabel="r (fm)", label="ρᵥₙ(r) calculated")
+p_vn = plot(rs(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

@@ -1,5 +1,5 @@
 using DelimitedFiles, Interpolations, Plots
-include("../nucleons.jl")
+include("../system.jl")
 
 # test data generated from Hartree.f
 # format: x S(x) V(x) R(x) A(x)
@@ -10,20 +10,24 @@ Vs = test_data[:, 3]
 Rs = test_data[:, 4]
 As = test_data[:, 5]
 
-S_interp = linear_interpolation(xs, Ss)
-V_interp = linear_interpolation(xs, Vs)
-R_interp = linear_interpolation(xs, Rs)
-A_interp = linear_interpolation(xs, As)
-
 p = true
-N = p ? 82 : 126
+N_p = 82
+N_n = 126
 r_max = maximum(xs)
-E_min = 860
+divs = length(xs) - 1
+s = system(N_p, N_n, r_max, divs)
+
+s.Φ0 = Ss
+s.W0 = Vs
+s.B0 = Rs
+s.A0 = As
+
+E_min = 800
 E_max = 939
 
-κs, Es = findAllOrbitals(p, S_interp, V_interp, R_interp, A_interp, r_max, E_min, E_max)
+κs, Es = findAllOrbitals(p, s, E_min, E_max)
 Ebinds = (p ? M_p : M_n) .- Es
-occ = fillNucleons(N, κs, Es)
+occ = fillNucleons(Z_or_N(s, p), κs, Es)
 
 # format: proton, kappa, filling, gnodes, fnodes, Ebind
 bench_data, _ = readdlm("test/Pb208Spectrum.csv", ','; header=true)

test/Pb208_nucleon_wf.jl

@@ -1,5 +1,5 @@
 using DelimitedFiles, Interpolations, Plots
-include("../nucleons.jl")
+include("../system.jl")
 
 # test data generated from Hartree.f
 # format: x S(x) V(x) R(x) A(x)
@@ -10,26 +10,27 @@ Vs = test_data[:, 3]
 Rs = test_data[:, 4]
 As = test_data[:, 5]
 
-S_interp = linear_interpolation(xs, Ss)
-V_interp = linear_interpolation(xs, Vs)
-R_interp = linear_interpolation(xs, Rs)
-A_interp = linear_interpolation(xs, As)
-
 κ = -1
 p = true
 r_max = maximum(xs)
-E_min = 880
+divs = length(xs) - 1
+s = system(r_max, divs)
+
+s.Φ0 = Ss
+s.W0 = Vs
+s.B0 = Rs
+s.A0 = As
+
+E_min = 800
 E_max = 939
 
-groundE = findEs(κ, p, S_interp, V_interp, R_interp, A_interp, r_max, E_min, E_max) |> minimum
+groundE = findEs(κ, p, s, E_min, E_max) |> minimum
 println("ground state E = $groundE")
 
-divs = 400
+wf = solveNucleonWf(κ, p, groundE, s)
-wf = solveNucleonWf(κ, p, groundE, S_interp, V_interp, R_interp, A_interp, r_max, divs)
-rs = range(0, r_max, length=divs+1)
 gs = wf[1, :]
 fs = wf[2, :]
 
-plot(rs, gs, label="g(r)")
+plot(rs(s), gs, label="g(r)")
-plot!(rs, fs, label="f(r)")
+plot!(rs(s), fs, label="f(r)")
 xlabel!("r (fm)")