diff --git a/nucleons.jl b/nucleons.jl
index 1acdf96..cc82506 100644
--- a/nucleons.jl
+++ b/nucleons.jl
@@ -118,11 +118,9 @@ j_κ(κ::Int) = abs(κ) - 1/2
 l_κ(κ::Int) = abs(κ) - (κ < 0) # since true = 1 and false = 0
 
 "Calculate scalar and vector densities of a nucleon species on [0,r_max] divided into (divs+1) points and returns them as vectors (ρ_s, ρ_v) where
+    the arrays κs, Es, occs tabulate the energies and occupation numbers corresponding to each κ,
     the other parameters are defined above"
-function calculateNucleonDensity(N, p, Φ0, W0, B0, A0, r_max, divs, E_min=0, E_max=(p ? M_p : M_n))
-    κs, Es = findAllOrbitals(p, Φ0, W0, B0, A0, r_max, E_min, E_max)
-    occs = fillNucleons(N, κs, Es)
-    
+function calculateNucleonDensity(κs, Es, occs, p, Φ0, W0, B0, A0, r_max, divs, E_min=0, E_max=(p ? M_p : M_n))
     ρr2 = zeros(2, divs + 1) # ρ×r² values
     
     for (κ, E, occ) in zip(κs, Es, occs)
@@ -139,4 +137,12 @@ function calculateNucleonDensity(N, p, Φ0, W0, B0, A0, r_max, divs, E_min=0, E_
     ρ_v = ρ[1, :] + ρ[2, :] # g^2 + f^2
     
     return (ρ_s, ρ_v)
-end
\ No newline at end of file
+end
+
+"Solve the Dirac equation and calculate scalar and vector densities of a nucleon species where
+    the other parameters are defined above"
+function solveNucleonDensity(N, p, Φ0, W0, B0, A0, r_max, divs, E_min=0, E_max=(p ? M_p : M_n))
+    κs, Es = findAllOrbitals(p, Φ0, W0, B0, A0, r_max, E_min, E_max)
+    occs = fillNucleons(N, κs, Es)
+    return calculateNucleonDensity(κs, Es, occs, p, Φ0, W0, B0, A0, r_max, divs, E_min, E_max)
+end
diff --git a/test/Pb208_nucleon_dens.jl b/test/Pb208_nucleon_dens.jl
index a95c7cb..c99c23f 100644
--- a/test/Pb208_nucleon_dens.jl
+++ b/test/Pb208_nucleon_dens.jl
@@ -24,12 +24,12 @@ divs = 400
 
 rs = range(0, r_max, length=divs+1)
 
-(ρ_sp, ρ_vp) = calculateNucleonDensity(N_p, true, S_interp, V_interp, R_interp, A_interp, r_max, divs, E_min, E_max)
+(ρ_sp, ρ_vp) = solveNucleonDensity(N_p, true, S_interp, V_interp, R_interp, A_interp, r_max, divs, E_min, E_max)
 
 p_sp = plot(rs, ρ_sp, xlabel="r (fm)", label="ρₛₚ(r) calculated")
 p_vp = plot(rs, ρ_vp, xlabel="r (fm)", label="ρᵥₚ(r) calculated")
 
-(ρ_sn, ρ_vn) = calculateNucleonDensity(N_n, false, S_interp, V_interp, R_interp, A_interp, r_max, divs, E_min, E_max)
+(ρ_sn, ρ_vn) = solveNucleonDensity(N_n, false, S_interp, V_interp, R_interp, A_interp, r_max, divs, E_min, E_max)
 
 p_sn = plot(rs, ρ_sn, xlabel="r (fm)", label="ρₛₙ(r) calculated")
 p_vn = plot(rs, ρ_vn, xlabel="r (fm)", label="ρᵥₙ(r) calculated")