51 lines
1.5 KiB
Julia
51 lines
1.5 KiB
Julia
using DelimitedFiles, Plots
|
|
include("../nucleons.jl")
|
|
|
|
# test data generated from Hartree.f
|
|
# format: x S(x) V(x) R(x) A(x)
|
|
test_data = readdlm("test/Pb208FldsFSUGarnet.csv")
|
|
xs = test_data[:, 1]
|
|
Ss = test_data[:, 2]
|
|
Vs = test_data[:, 3]
|
|
Rs = test_data[:, 4]
|
|
As = test_data[:, 5]
|
|
|
|
p = true
|
|
N_p = 82
|
|
N_n = 126
|
|
r_max = maximum(xs)
|
|
divs = length(xs) - 1
|
|
s = system(N_p, N_n, parameters(), r_max, divs)
|
|
|
|
s.Φ0 = Ss
|
|
s.W0 = Vs
|
|
s.B0 = Rs
|
|
s.A0 = As
|
|
|
|
κs, Es = findAllOrbitals(p, s)
|
|
Ebinds = (p ? M_p : M_n) .- Es
|
|
spec = fillNucleons(Z_or_N(s, p), κs, Es)
|
|
|
|
# format: proton, kappa, filling, gnodes, fnodes, Ebind
|
|
bench_data, _ = readdlm("test/Pb208Spectrum.csv", ','; header=true)
|
|
mask = findall(bench_data[:, 1] .== Int(p))
|
|
bench_κs = bench_data[mask, 2]
|
|
bench_Ebinds = bench_data[mask, 6]
|
|
bench_filling = bench_data[mask, 3]
|
|
bench_occ = Int.(bench_filling .* (2 .* j_κ.(Int.(bench_κs)) .+ 1))
|
|
|
|
scatter(κs, Ebinds, label=(p ? "proton" : "neutron") * " spectrum")
|
|
annotate!(κs .+ 0.3, Ebinds .+ 1, [(o, 9) for o in spec.occ])
|
|
|
|
scatter!(bench_κs, bench_Ebinds, label="benchmark"; markershape=:cross, markercolor="red")
|
|
annotate!(bench_κs .+ 0.3, bench_Ebinds .- 1, [(o, 9, :red) for o in bench_occ])
|
|
|
|
# label each κ in spectroscopic notation l_j (GR renders $...$ as LaTeX)
|
|
spectroscopic(κ::Int) = "\$\\mathrm{" * "spdfghiklmno"[l_κ(κ) + 1] * "}_{" * string(Int(2 * j_κ(κ))) * "/2}\$"
|
|
tick_κs = sort(unique(Int.([κs; bench_κs])))
|
|
xticks!(tick_κs, spectroscopic.(tick_κs))
|
|
xlabel!("Orbital")
|
|
ylabel!("Binding energy (MeV)")
|
|
yflip!()
|
|
savefig("temp/Pb208_nucleon_spectrum.pdf")
|