diff --git a/dirac.jl b/dirac.jl
index 4e84fc9..d711057 100644
--- a/dirac.jl
+++ b/dirac.jl
@@ -38,3 +38,26 @@ function findEs(κ, p, Φ0, W0, B0, A0, r_max, E_min=0, E_max=(p ? M_p : M_n))
     f(E) = boundaryValue(κ, p, E, Φ0, W0, B0, A0, r_max)
     return find_zeros(f, (E_min, E_max))
 end
+
+"Find N lowest lying orbitals and return two lists containing κ values and corresponding energies for a single species where
+    the other parameters are defined above"
+function fillNucleons(N, p, Φ0, W0, B0, A0, r_max, E_min=0, E_max=(p ? M_p : M_n))
+    κs = Int[]
+    Es = Float64[]
+    # start from κ=-1 and go both up and down
+    for direction in [-1, 1]
+        for κ in direction * (1:100) # cutoff is 100
+            new_Es = findEs(κ, p, Φ0, W0, B0, A0, r_max, E_min, E_max)
+            if isempty(new_Es); break; end
+            append!(Es, new_Es)
+            append!(κs, fill(κ, length(new_Es)))
+        end
+    end
+    if length(Es) < N
+        @warn "Only found $(length(Es)) orbitals for $N nucleons"
+        return (κs, Es)
+    else
+        idxs = sortperm(Es)[1:N]
+        return (κs[idxs], Es[idxs])
+    end
+end
\ No newline at end of file
diff --git a/test/Pb208.jl b/test/Pb208_basic.jl
similarity index 100%
rename from test/Pb208.jl
rename to test/Pb208_basic.jl
diff --git a/test/Pb208_spectrum.jl b/test/Pb208_spectrum.jl
new file mode 100644
index 0000000..f2aeefc
--- /dev/null
+++ b/test/Pb208_spectrum.jl
@@ -0,0 +1,28 @@
+using DelimitedFiles, Interpolations, Plots
+include("../dirac.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]
+
+S_interp = linear_interpolation(xs, Ss)
+V_interp = linear_interpolation(xs, Vs)
+R_interp = linear_interpolation(xs, Rs)
+A_interp = linear_interpolation(xs, As)
+
+p = true
+N = 82
+r_max = maximum(xs)
+E_min = 880
+E_max = 939
+
+κs, Es = fillNucleons(N, p, S_interp, V_interp, R_interp, A_interp, r_max, E_min, E_max)
+
+scatter(κs, Es, label=(p ? "proton" : "neutron") * " spectrum")
+xlabel!("κ")
+ylabel!("E (MeV)")
\ No newline at end of file