using DelimitedFiles, Plots
include("../system.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]

N_p = 82
N_n = 126
r_max = maximum(xs)
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

(s.ρ_sp, s.ρ_vp) = solveNucleonDensity(true, s)

p_sp = plot(rs(s), s.ρ_sp, xlabel="r (fm)", label="ρₛₚ(r) calculated")
p_vp = plot(rs(s), s.ρ_vp, xlabel="r (fm)", label="ρᵥₚ(r) calculated")

(s.ρ_sn, s.ρ_vn) = solveNucleonDensity(false, s)

p_sn = plot(rs(s), s.ρ_sn, xlabel="r (fm)", label="ρₛₙ(r) calculated")
p_vn = plot(rs(s), 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)
bench_data = readdlm("test/Pb208DensFSUGarnet.csv")
xs = bench_data[:, 1]
rho_sn = bench_data[:, 2]
rho_vn = bench_data[:, 3]
rho_sp = bench_data[:, 5]
rho_vp = bench_data[:, 6]

plot!(p_sp, xs, rho_sp, label="ρₛₚ(r) benchmark")
plot!(p_vp, xs, rho_vp, label="ρᵥₚ(r) benchmark")
plot!(p_sn, xs, rho_sn, label="ρₛₙ(r) benchmark")
plot!(p_vn, xs, rho_vn, label="ρᵥₙ(r) benchmark")

plot(p_sp, p_vp, p_sn, p_vn, layout=(2, 2))