yapanuwan
/
NuclearRMF
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

All functions take s::system

This commit is contained in:
Nuwan Yapa 2025-02-17 19:32:19 -05:00
parent 95c1083693
commit 6850fdf351
9 changed files with 97 additions and 90 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

27
mesons.jl
View File

@ -12,7 +12,7 @@ const g2_s = 110.349189097820 # dimensionless
const g2_v = 187.694676506801 # dimensionless const g2_v = 187.694676506801 # dimensionless
const g2_ρ = 192.927428365698 # dimensionless const g2_ρ = 192.927428365698 # dimensionless
const g2_γ = 0.091701236 # dimensionless # equal to 4πα const g2_γ = 0.091701236 # dimensionless # equal to 4πα
const κ = 3.260178893447 # MeV const κ_ss = 3.260178893447 # MeV
const λ = -0.003551486718 # dimensionless # LambdaSS const λ = -0.003551486718 # dimensionless # LambdaSS
const ζ = 0.023499504053 # dimensionless # LambdaVV const ζ = 0.023499504053 # dimensionless # LambdaVV
const Λv = 0.043376933644 # dimensionless # LambdaVR const Λv = 0.043376933644 # dimensionless # LambdaVR
@ -26,15 +26,12 @@ greensFunctionKG(m, r, rp) = -1 / (m * (r + r_reg) * (rp + r_reg)) * sinh(m * mi
greensFunctionP(r, rp) = -1 / (max(r, rp) + r_reg) 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 "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) 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) return greenMat * (int_measure .* source)
end end
@ -47,21 +44,21 @@ end
r is the radius in fm, r is the radius in fm,
An initial guess initial_sol=(Φ0, W0, B0, A0) can be provided to speed up convergence (permuting!). An initial guess initial_sol=(Φ0, W0, B0, A0) can be provided to speed up convergence (permuting!).
Reference: P. Giuliani, K. Godbey, E. Bonilla, F. Viens, and J. Piekarewicz, Frontiers in Physics 10, (2023)" 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 solveMesonWfs(ρ_sp, ρ_vp, ρ_sn, ρ_vn, s::system, iterations=10; initial_sol=(zeros(1 + divs) for _ in 1:4))
(Φ0, W0, B0, A0) = initial_sol (Φ0, W0, B0, A0) = initial_sol
(src_Φ0, src_W0, src_B0, src_A0) = (zeros(1 + divs) for _ in 1:4) (src_Φ0, src_W0, src_B0, src_A0) = (zero_array(s) for _ in 1:4)
# A0 doesn't need iterations # A0 doesn't need iterations
@. src_A0 = -g2_γ * ρ_vp @. src_A0 = -g2_γ * ρ_vp
A0 .= solveKG(m_γ, src_A0, r_max) A0 .= solveKG(m_γ, src_A0, s)
for _ in 1:iterations 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_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 @. 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 end
return (Φ0, W0, B0, A0) return (Φ0, W0, B0, A0)

62
nucleons.jl
View File

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

24
system.jl
View File

@ -1,6 +1,4 @@
using Interpolations, PolyLog, Plots using Interpolations, PolyLog, Plots
include("nucleons.jl")
include("mesons.jl")
"Defines a nuclear system to be solved" "Defines a nuclear system to be solved"
struct system struct system
@ -11,13 +9,16 @@ struct system
end end
"Get mass number of nucleus" "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" "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" "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" "Create an empty array for the size of the mesh"
zero_array(s::system) = zeros(1 + s.divs) zero_array(s::system) = zeros(1 + s.divs)
@ -25,6 +26,9 @@ zero_array(s::system) = zeros(1 + s.divs)
"Normalized Woods-Saxon form used for constructing an initial solution" "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))) 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)))
include("nucleons.jl")
include("mesons.jl")
"Run the full program by self-consistent solution of nucleon and meson densities" "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) 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) if isnothing(initial_dens)
@ -46,7 +50,7 @@ function solve_system(s::system, initial_dens=nothing, initial_flds=(zeros(1 + s
E_total_previous = NaN E_total_previous = NaN
while true 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)) @time "Meson fields" (Φ0s, W0s, B0s, A0s) = solveMesonWfs(ρ_sp, ρ_vp, ρ_sn, ρ_vn, s, isnan(E_total_previous) ? 50 : 5; initial_sol = (Φ0s, W0s, B0s, A0s))
S_interp = linear_interpolation(rs(s), Φ0s) S_interp = linear_interpolation(rs(s), Φ0s)
V_interp = linear_interpolation(rs(s), W0s) V_interp = linear_interpolation(rs(s), W0s)
@ -54,14 +58,14 @@ function solve_system(s::system, initial_dens=nothing, initial_flds=(zeros(1 + s
A_interp = linear_interpolation(rs(s), A0s) A_interp = linear_interpolation(rs(s), A0s)
# protons # 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_interp, V_interp, R_interp, A_interp, s)
occs_p = fillNucleons(s.Z, κs_p, Es_p) 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" (ρ_sp, ρ_vp) = calculateNucleonDensity(κs_p, Es_p, occs_p, true, S_interp, V_interp, R_interp, A_interp, s)
# neutrons # 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_interp, V_interp, R_interp, A_interp, s)
occs_n = fillNucleons(s.N, κs_n, Es_n) 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" (ρ_sn, ρ_vn) = calculateNucleonDensity(κs_n, Es_n, occs_n, false, S_interp, V_interp, R_interp, A_interp, s)
if monitor_plot if monitor_plot
for s in p.series_list for s in p.series_list

10
test/Linear_nucleon_spectrum.jl
View File

@ -1,5 +1,5 @@
using DelimitedFiles, Interpolations, Plots using DelimitedFiles, Interpolations, Plots
include("../nucleons.jl") include("../system.jl")
# test data generated from Hartree.f # test data generated from Hartree.f
# format: x S(x) V(x) R(x) A(x) # format: x S(x) V(x) R(x) A(x)
@ -16,14 +16,16 @@ R_interp = linear_interpolation(xs, Rs)
A_interp = linear_interpolation(xs, As) A_interp = linear_interpolation(xs, As)
p = false p = false
N = p ? 8 : 8
r_max = maximum(xs) r_max = maximum(xs)
divs = length(xs) - 1
s = system(8, 8, r_max, divs)
E_min = 860 E_min = 860
E_max = 939 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_interp, V_interp, R_interp, A_interp, s, E_min, E_max)
Ebinds = (p ? M_p : M_n) .- Es 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 # format: proton, kappa, filling, gnodes, fnodes, Ebind
bench_data, _ = readdlm("test/LinearSpectrum.csv", ','; header=true) bench_data, _ = readdlm("test/LinearSpectrum.csv", ','; header=true)

5
test/Pb208_meson_flds.jl
View File

@ -1,5 +1,5 @@
using DelimitedFiles, Interpolations, Plots using DelimitedFiles, Interpolations, Plots
include("../mesons.jl") include("../system.jl")
# test data generated from Hartree.f # test data generated from Hartree.f
# format: x S(x) V(x) R(x) A(x) # format: x S(x) V(x) R(x) A(x)
@ -23,8 +23,9 @@ xs = test_data[:, 1]
r_max = maximum(xs) r_max = maximum(xs)
divs = length(xs) - 1 divs = length(xs) - 1
s = system(0, 0, r_max, divs)
(Φ0, W0, B0, A0) = solveMesonWfs(ρ_sp, ρ_vp, ρ_sn, ρ_vn, r_max, divs, 200) # many iterations needed without an initial solution (Φ0, W0, B0, A0) = solveMesonWfs(ρ_sp, ρ_vp, ρ_sn, ρ_vn, s, 200) # many iterations needed without an initial solution
plot!(xs, hcat(Φ0, W0, B0, A0), layout=4, label=["Φ0" "W0" "B0" "A0"]) plot!(xs, hcat(Φ0, W0, B0, A0), layout=4, label=["Φ0" "W0" "B0" "A0"])
xlabel!("r (fm)") xlabel!("r (fm)")

9
test/Pb208_nucleon_basic.jl
View File

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

21
test/Pb208_nucleon_dens.jl
View File

@ -1,5 +1,5 @@
using DelimitedFiles, Interpolations, Plots using DelimitedFiles, Interpolations, Plots
include("../nucleons.jl") include("../system.jl")
# test data generated from Hartree.f # test data generated from Hartree.f
# format: x S(x) V(x) R(x) A(x) # format: x S(x) V(x) R(x) A(x)
@ -18,21 +18,18 @@ A_interp = linear_interpolation(xs, As)
N_p = 82 N_p = 82
N_n = 126 N_n = 126
r_max = maximum(xs) 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) (ρ_sp, ρ_vp) = solveNucleonDensity(true, S_interp, V_interp, R_interp, A_interp, s)
(ρ_sp, ρ_vp) = solveNucleonDensity(N_p, true, S_interp, V_interp, R_interp, A_interp, r_max, divs, E_min, E_max) p_sp = plot(rs(s), ρ_sp, xlabel="r (fm)", label="ρₛₚ(r) calculated")
p_vp = plot(rs(s), ρ_vp, xlabel="r (fm)", label="ρᵥₚ(r) calculated")
p_sp = plot(rs, ρ_sp, xlabel="r (fm)", label="ρₛₚ(r) calculated") (ρ_sn, ρ_vn) = solveNucleonDensity(false, S_interp, V_interp, R_interp, A_interp, s)
p_vp = plot(rs, ρ_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) p_sn = plot(rs(s), ρ_sn, xlabel="r (fm)", label="ρₛₙ(r) calculated")
p_vn = plot(rs(s), ρ_vn, xlabel="r (fm)", label="ρᵥₙ(r) calculated")
p_sn = plot(rs, ρ_sn, xlabel="r (fm)", label="ρₛₙ(r) calculated")
p_vn = plot(rs, ρ_vn, xlabel="r (fm)", label="ρᵥₙ(r) calculated")
# benchmark data generated from Hartree.f # benchmark data generated from Hartree.f
# format: x Rhos(n) Rhov(n) Rhot(n) Rhos(p) Rhov(p) Rhot(p) # format: x Rhos(n) Rhov(n) Rhot(n) Rhos(p) Rhov(p) Rhot(p)

13
test/Pb208_nucleon_spectrum.jl
View File

@ -1,5 +1,5 @@
using DelimitedFiles, Interpolations, Plots using DelimitedFiles, Interpolations, Plots
include("../nucleons.jl") include("../system.jl")
# test data generated from Hartree.f # test data generated from Hartree.f
# format: x S(x) V(x) R(x) A(x) # format: x S(x) V(x) R(x) A(x)
@ -16,14 +16,17 @@ R_interp = linear_interpolation(xs, Rs)
A_interp = linear_interpolation(xs, As) A_interp = linear_interpolation(xs, As)
p = true p = true
N = p ? 82 : 126 N_p = 82
N_n = 126
r_max = maximum(xs) r_max = maximum(xs)
E_min = 860 divs = length(xs) - 1
s = system(N_p, N_n, r_max, divs)
E_min = 800
E_max = 939 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_interp, V_interp, R_interp, A_interp, s, E_min, E_max)
Ebinds = (p ? M_p : M_n) .- Es 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 # format: proton, kappa, filling, gnodes, fnodes, Ebind
bench_data, _ = readdlm("test/Pb208Spectrum.csv", ','; header=true) bench_data, _ = readdlm("test/Pb208Spectrum.csv", ','; header=true)

16
test/Pb208_nucleon_wf.jl
View File

@ -1,5 +1,5 @@
using DelimitedFiles, Interpolations, Plots using DelimitedFiles, Interpolations, Plots
include("../nucleons.jl") include("../system.jl")
# test data generated from Hartree.f # test data generated from Hartree.f
# format: x S(x) V(x) R(x) A(x) # format: x S(x) V(x) R(x) A(x)
@ -18,18 +18,18 @@ A_interp = linear_interpolation(xs, As)
κ = -1 κ = -1
p = true p = true
r_max = maximum(xs) r_max = maximum(xs)
E_min = 880 divs = length(xs) - 1
s = system(0, 0, r_max, divs)
E_min = 800
E_max = 939 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_interp, V_interp, R_interp, A_interp, s, E_min, E_max) |> minimum
println("ground state E = $groundE") println("ground state E = $groundE")
divs = 400 wf = solveNucleonWf(κ, p, groundE, S_interp, V_interp, R_interp, A_interp, 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, :] gs = wf[1, :]
fs = wf[2, :] 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)") xlabel!("r (fm)")