From 3f1c8803bb9b515270e6dc75555c7a62458a1418 Mon Sep 17 00:00:00 2001
From: Nuwan Yapa <yapanuwan@gmail.com>
Date: Mon, 3 Feb 2025 20:31:37 -0500
Subject: [PATCH] Nuclear spectrum benchmarking

---
 test/Pb208Spectrum.csv         | 39 ++++++++++++++++++++++++++++++++++
 test/Pb208_nucleon_spectrum.jl | 19 ++++++++++++++---
 2 files changed, 55 insertions(+), 3 deletions(-)
 create mode 100644 test/Pb208Spectrum.csv

diff --git a/test/Pb208Spectrum.csv b/test/Pb208Spectrum.csv
new file mode 100644
index 0000000..be25139
--- /dev/null
+++ b/test/Pb208Spectrum.csv
@@ -0,0 +1,39 @@
+proton,kappa,filling,gnodes,fnodes,Ebind
+0,-1,1,0,0,63.80958
+0,-2,1,0,0,57.09607
+0,1,1,0,1,56.51421
+0,-3,1,0,0,48.96921
+0,2,1,0,1,47.53564
+0,-1,1,1,1,44.05962
+0,-4,1,0,0,39.82046
+0,3,1,0,1,37.11666
+0,-2,1,1,1,32.52852
+0,1,1,1,2,31.35759
+0,-5,1,0,0,30.03271
+0,4,1,0,1,25.74641
+0,-3,1,1,1,21.54695
+0,-6,1,0,0,19.94124
+0,2,1,1,2,19.70747
+0,-1,1,2,2,19.05405
+0,5,1,0,1,13.96937
+0,-4,1,1,1,11.33989
+0,-7,1,0,0,9.84013
+0,3,1,1,2,9.04831
+0,-2,1,2,2,8.57808
+0,1,1,2,3,7.67567
+1,-1,1,0,0,50.48167
+1,-2,1,0,0,44.63156
+1,1,1,0,1,43.88757
+1,-3,1,0,0,37.28703
+1,2,1,0,1,35.553
+1,-1,1,1,1,31.53505
+1,-4,1,0,0,28.92962
+1,3,1,0,1,25.80862
+1,-2,1,1,1,21.04461
+1,-5,1,0,0,19.94885
+1,1,1,1,2,19.8522
+1,4,1,0,1,15.19352
+1,-3,1,1,1,11.01997
+1,-6,1,0,0,10.64656
+1,2,1,1,2,9.22511
+1,-1,1,2,2,8.31272
\ No newline at end of file
diff --git a/test/Pb208_nucleon_spectrum.jl b/test/Pb208_nucleon_spectrum.jl
index 282015b..10c0965 100644
--- a/test/Pb208_nucleon_spectrum.jl
+++ b/test/Pb208_nucleon_spectrum.jl
@@ -22,9 +22,22 @@ E_min = 880
 E_max = 939
 
 κs, Es = findAllOrbitals(p, S_interp, V_interp, R_interp, A_interp, r_max, E_min, E_max)
+Ebinds = (p ? M_p : M_n) .- Es
 occ = fillNucleons(N, κs, Es)
 
-scatter(κs, Es, label=(p ? "proton" : "neutron") * " spectrum")
-annotate!(κs .+ 0.3, Es, occ)
+# 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 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])
+
 xlabel!("κ")
-ylabel!("E (MeV)")
\ No newline at end of file
+ylabel!("E (MeV)")