Single particle energies matched with Hartree.f

mesons.jl

@@ -1,6 +1,6 @@
 using DifferentialEquations
 
-const ħc = 197.327 # MeVfm
+const ħc = 197.33 # MeVfm
 
 # Values defined in C. J. Horowitz and J. Piekarewicz, Phys. Rev. Lett. 86, 5647 (2001)
 # Values taken from Hartree.f (FSUGarnet)
@@ -17,7 +17,7 @@ const λ = -0.003551486718 # dimensionless # LambdaSS
 const ζ = 0.023499504053 # dimensionless # LambdaVV
 const Λv = 0.043376933644 # dimensionless # LambdaVR
 
-const r_reg = 1E-9 # fm # regulator for Green's functions
+const r_reg = 1E-8 # fm # regulator for Green's functions
 
 "Green's function for Klein-Gordon equation in natural units"
 greensFunctionKG(m, r, rp) = -1 / (m * (r + r_reg) * (rp + r_reg)) * sinh(m * min(r, rp)) * exp(-m * max(r, rp))

nucleons.jl

@@ -1,10 +1,10 @@
 using DifferentialEquations, Roots
 
-const ħc = 197.327 # MeVfm
-const M_n = 939.5654133 # MeV/c2
-const M_p = 938.2720813 # MeV/c2
+const ħc = 197.33 # MeVfm
+const M_n = 939.0 # MeV/c2
+const M_p = 939.0 # MeV/c2
 
-const r_reg = 1E-6 # fm # regulator for the centrifugal term
+const r_reg = 1E-8 # fm # regulator for the centrifugal term
 
 "The spherical Dirac equation that returns du=[dg, df] in-place where
     u=[g, f] are the reduced radial components evaluated at r,

test/Pb208_nucleon_basic.jl

@@ -16,13 +16,15 @@ R_interp = linear_interpolation(xs, Rs)
 A_interp = linear_interpolation(xs, As)
 
 κ = -1
-p = true
+p = false
 r_max = maximum(xs)
-E_min = 880
+E_min = 850
 E_max = 939
 
 boundEs = findEs(κ, p, S_interp, V_interp, R_interp, A_interp, r_max, E_min, E_max)
-println("bound E = $boundEs")
+
+binding_Es = round.((p ? M_p : M_n) .- boundEs; digits=3) |> unique
+println("binding energies = $binding_Es")
 
 Es = collect(E_min:0.5:E_max)
 boundaryVals = [boundaryValue(κ, p, E, S_interp, V_interp, R_interp, A_interp, r_max)^2 for E in Es]