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/LinearFldsFSUGarnet.csv") 
xs = test_data[:, 1]
Ss = test_data[:, 2]
Vs = test_data[:, 3]
Rs = test_data[:, 4]
As = test_data[:, 5]

p = false
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

κs, Es = findAllOrbitals(p, s)
Ebinds = (p ? M_p : M_n) .- 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)
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!("Binding energy (MeV)")
yflip!()