diff --git a/mesons.jl b/mesons.jl index 7c7fba9..79e1903 100644 --- a/mesons.jl +++ b/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, - r_max is the r-cutoff in fm." -function solveKG(m, source, r_max) - Δr = r_max / (length(source) - 1) - rs = range(0, r_max; length=length(source)) - int_measure = ħc .* Δr .* rs .^ 2 + m is the mass of the meson in MeV/c2." +function solveKG(m, source, s::system) + int_measure = ħc .* Δr(s) .* rs(s) .^ 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)) - (Φ0, W0, B0, A0) = initial_sol - (src_Φ0, src_W0, src_B0, src_A0) = (zeros(1 + divs) for _ in 1:4) +function solveMesonFields!(s::system, iterations=10) + (Φ0, W0, B0, A0) = (s.Φ0, s.W0, s.B0, s.A0) + (ρ_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) - W0 .= (solveKG(m_ω, src_W0, r_max) .+ W0) ./ 2 # to supress oscillation - B0 .= (solveKG(m_ρ, src_B0, r_max) .+ B0) ./ 2 # to supress oscillation + Φ0 .= solveKG(m_s, src_Φ0, s) + W0 .= (solveKG(m_ω, src_W0, s) .+ W0) ./ 2 # to supress oscillation + B0 .= (solveKG(m_ρ, src_B0, s) .+ B0) ./ 2 # to supress oscillation end return (Φ0, W0, B0, A0) diff --git a/nucleons.jl b/nucleons.jl index 7abc274..49c4f64 100644 --- a/nucleons.jl +++ b/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()) - Δr = r_max / divs +function solveNucleonWf(κ, p, E, s::system; shooting=true, normalize=true, algo=Tsit5()) + (Φ0, W0, B0, A0) = fields_interp(s) if shooting - @assert divs % 2 == 0 "divs must be an even number when shooting=true" - prob = ODEProblem(dirac!, [0, 1], (r_max, r_max / 2)) - sol = solve(prob, algo, p=(κ, p, E, Φ0, W0, B0, A0), saveat=Δr) + @assert s.divs % 2 == 0 "divs must be an even number when shooting=true" + 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(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)) - sol = solve(prob, algo, p=(κ, p, E, Φ0, W0, B0, A0), saveat=Δr) + prob = ODEProblem(dirac!, [0, 1], (0, next_r_max)) + 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()) - prob = ODEProblem(dirac!, convert.(dtype, [0, 1]), (0, r_max)) - func(E) = solve(prob, algo, p=(κ, p, E, Φ0, W0, B0, A0), saveat=[r_max], save_idxs=[1])[1, 1] +function boundaryValueFunc(κ, p, s::system; dtype=Float64, algo=Tsit5()) + prob = ODEProblem(dirac!, convert.(dtype, [0, 1]), (0, s.r_max)) + (Φ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) - func = boundaryValueFunc(κ, p, Φ0, W0, B0, A0, r_max) +function findEs(κ, p, s::system, E_min=0, E_max=(p ? M_p : M_n), tol_digits=5) + 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" -function fillNucleons(N::Int, κs, Es) +"For a given list of κ values with corresponding energies, attempt to fill Z_or_N lowest lying orbitals and return occupancy numbers" +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) - N -= occ[i] + occ[i] = min(max_occ, Z_or_N) + 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) - ρr2 = zeros(2, divs + 1) # ρ×r² values +function calculateNucleonDensity(κs, Es, occs, p, s::system) + ρ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) - κs, Es = findAllOrbitals(p, Φ0, W0, B0, A0, r_max, E_min, E_max) - occs = fillNucleons(N, κs, Es) - return calculateNucleonDensity(κs, Es, occs, p, Φ0, W0, B0, A0, r_max, divs) +function solveNucleonDensity(p, s::system, E_min=800, E_max=939) + κs, Es = findAllOrbitals(p, s, E_min, E_max) + occs = fillNucleons(Z_or_N(s, p), κs, Es) + return calculateNucleonDensity(κs, Es, occs, p, s) end diff --git a/system.jl b/system.jl index 4b425ef..c3a2531 100644 --- a/system.jl +++ b/system.jl @@ -1,23 +1,41 @@ using Interpolations, PolyLog, Plots -include("nucleons.jl") -include("mesons.jl") -"Defines a nuclear system to be solved" -struct system +"Defines a nuclear system containing relevant parameters and meson/nucleon densities" +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" -function solve_system(s::system, initial_dens=nothing, initial_flds=(zeros(1 + s.divs) for _ in 1:4); monitor_print=true, monitor_plot=false) - if isnothing(initial_dens) - dens_guess = Woods_Saxon.(rs(s)) - ρ_sp = s.Z .* dens_guess - ρ_vp = s.Z .* dens_guess - ρ_sn = s.N .* dens_guess - ρ_vn = s.N .* dens_guess - else - (ρ_sp, ρ_vp, ρ_sn, ρ_vn) = initial_dens - end +"Create linear interpolations for the meson fields" +function fields_interp(s::system) + S_interp = linear_interpolation(rs(s), s.Φ0) + V_interp = linear_interpolation(rs(s), s.W0) + R_interp = linear_interpolation(rs(s), s.B0) + A_interp = linear_interpolation(rs(s), s.A0) + return (S_interp, V_interp, R_interp, A_interp) +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)) - - 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) + # mesons + @time "Meson fields" solveMesonFields!(s, isnan(E_total_previous) ? 50 : 5) # 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 diff --git a/test/Linear_nucleon_spectrum.jl b/test/Linear_nucleon_spectrum.jl index 2f7b92f..5714e9c 100644 --- a/test/Linear_nucleon_spectrum.jl +++ b/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) diff --git a/test/Pb208_meson_flds.jl b/test/Pb208_meson_flds.jl index e360d5f..c438106 100644 --- a/test/Pb208_meson_flds.jl +++ b/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)") \ No newline at end of file diff --git a/test/Pb208_nucleon_basic.jl b/test/Pb208_nucleon_basic.jl index 481d55b..a5f2edf 100644 --- a/test/Pb208_nucleon_basic.jl +++ b/test/Pb208_nucleon_basic.jl @@ -1,5 +1,5 @@ -using DelimitedFiles, Interpolations, Plots -include("../nucleons.jl") +using DelimitedFiles, Plots +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] diff --git a/test/Pb208_nucleon_dens.jl b/test/Pb208_nucleon_dens.jl index c99c23f..9c1c277 100644 --- a/test/Pb208_nucleon_dens.jl +++ b/test/Pb208_nucleon_dens.jl @@ -1,5 +1,5 @@ -using DelimitedFiles, Interpolations, Plots -include("../nucleons.jl") +using DelimitedFiles, Plots +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 -E_max = 939 -divs = 400 +divs = length(xs) - 1 +s = system(N_p, N_n, r_max, divs) -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_vp = plot(rs, ρ_vp, xlabel="r (fm)", label="ρᵥₚ(r) calculated") +p_sp = plot(rs(s), ρ_sp, 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_vn = plot(rs, ρ_vn, xlabel="r (fm)", label="ρᵥₙ(r) calculated") +p_sn = plot(rs(s), ρ_sn, 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) diff --git a/test/Pb208_nucleon_spectrum.jl b/test/Pb208_nucleon_spectrum.jl index 18992eb..c12ea6d 100644 --- a/test/Pb208_nucleon_spectrum.jl +++ b/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) diff --git a/test/Pb208_nucleon_wf.jl b/test/Pb208_nucleon_wf.jl index 0a13905..7898949 100644 --- a/test/Pb208_nucleon_wf.jl +++ b/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_interp, V_interp, R_interp, A_interp, r_max, divs) -rs = range(0, r_max, length=divs+1) +wf = solveNucleonWf(κ, p, groundE, s) gs = wf[1, :] fs = wf[2, :] -plot(rs, gs, label="g(r)") -plot!(rs, fs, label="f(r)") +plot(rs(s), gs, label="g(r)") +plot!(rs(s), fs, label="f(r)") xlabel!("r (fm)")