Fixed some units

mesons.jl

@@ -1,34 +1,35 @@
 using DifferentialEquations
 
-const ħc = 197.327 # ħc in MeVfm
+const ħc = 197.327 # MeVfm
 
 # Values defined in C. J. Horowitz and J. Piekarewicz, Phys. Rev. Lett. 86, 5647 (2001)
 # Values taken from Hartree.f (FSUGarnet) and
-const m_ρ = 763.0
-const m_ω = 782.5
-const m_s = 496.939473213388
-const m_γ = 0.000001000 # not defined in paper
-const g2_s = 110.349189097820
-const g2_v = 187.694676506801
-const g2_ρ = 192.927428365698
-const g2_g = 0.091701236 # not defined in paper
-const κ = 3.260178893447
-const λ = -0.003551486718 # LambdaSS
-const ζ = 0.023499504053 # LambdaVV
-const Λv = 0.043376933644 # LambdaVR
+const m_ρ = 763.0 # MeV/c2
+const m_ω = 782.5 # MeV/c2
+const m_s = 496.939473213388 # MeV/c2
+const m_γ = 0.000001000 # MeV/c2 # not defined in paper
+const g2_s = 110.349189097820 # units?
+const g2_v = 187.694676506801 # units?
+const g2_ρ = 192.927428365698 # units?
+const g2_g = 0.091701236 # units? # not defined in paper
+const κ = 3.260178893447 # units?
+const λ = -0.003551486718 # units? # LambdaSS
+const ζ = 0.023499504053 # units? # LambdaVV
+const Λv = 0.043376933644 # units? # LambdaVR
 
-const r_reg = 1E-9 # regulator for Green's functions in fm
+const r_reg = 1E-9 # fm # regulator for Green's functions
 
-"Green's function for Klein-Gordon equation"
+"Green's function for Klein-Gordon equation in natural units"
 greensFunctionKG(m, r, rp) = -(1 / m) * rp / (r + r_reg) * sinh(m * min(r, rp)) * exp(-m * max(r, rp))
 
-"Green's function for Poisson's equation"
+"Green's function for Poisson's equation in natural units"
 greensFunctionP(r, rp) = -rp^2 / (max(r, rp) + r_reg)
 
 "Solve the Klein-Gordon equation (or Poisson's equation when m=0) for a given a source function (as an array) where
-    the other parameters are the same from mesons!(...)."
+    m is the mass of the meson in MeV/c2,
+    r_max is the r-cutoff in fm."
 function solveKG(m, source, r_max)
-    greensFunction = m == 0 ? greensFunctionP : (r, rp) -> greensFunctionKG(m, r, rp)
+    greensFunction = m == 0 ? greensFunctionP : (r, rp) -> greensFunctionKG(m / ħc, r, rp)
     mesh = range(0, r_max; length=length(source))
     greenMat = greensFunction.(transpose(mesh), mesh)
     return greenMat * source
@@ -50,14 +51,14 @@ function solveMesonWfs(ρ_sp, ρ_vp, ρ_sn, ρ_vn, r_max, divs, iterations=3)
     (src_Φ0, src_W0, src_B0) = (zeros(1 + divs) for _ in 1:3)
 
     for _ in 1:iterations
-        @. src_Φ0 = g2_s * ((κ/2) * Φ0^2 + (λ/6) * Φ0^3 - (ρ_sp + ρ_sn)) / ħc
-        Φ0 .= solveKG(m_s / sqrt(ħc), src_Φ0, r_max)
+        @. src_Φ0 = g2_s * ((κ/2) * Φ0^2 + (λ/6) * Φ0^3 - (ρ_sp + ρ_sn))
+        Φ0 .= solveKG(m_s, src_Φ0, r_max)
 
-        @. src_W0 = g2_v * ((ζ/6) * W0^3 + 2 * Λv * B0^2 * W0 - (ρ_vp + ρ_vn)) / ħc
-        W0 .= solveKG(m_ω / sqrt(ħc), src_W0, r_max)
+        @. src_W0 = g2_v * ((ζ/6) * W0^3 + 2 * Λv * B0^2 * W0 - (ρ_vp + ρ_vn))
+        W0 .= solveKG(m_ω, src_W0, r_max)
 
-        @. src_B0 = g2_ρ * (2 * Λv * W0^2 * B0 - (ρ_vp - ρ_vn) / 2) / ħc
-        B0 .= solveKG(m_ρ / sqrt(ħc), src_B0, r_max)
+        @. src_B0 = g2_ρ * (2 * Λv * W0^2 * B0 - (ρ_vp - ρ_vn) / 2)
+        B0 .= solveKG(m_ρ, src_B0, r_max)
     end
 
     return (Φ0, W0, B0, A0)

nucleons.jl

@@ -1,10 +1,10 @@
 using DifferentialEquations, Roots
 
-const ħc = 197.327 # ħc in MeVfm
-const M_n = 939.5654133 # Neutron mass in MeV/c2
-const M_p = 938.2720813 # Proton mass in MeV/c2
+const ħc = 197.327 # MeVfm
+const M_n = 939.5654133 # MeV/c2
+const M_p = 938.2720813 # MeV/c2
 
-const r_reg = 1E-6 # regulator for the centrifugal term in fm
+const r_reg = 1E-6 # 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_meson_flds.jl

@@ -13,7 +13,7 @@ xs = test_data[:, 1]
 r_max = maximum(xs)
 divs = length(xs) - 1
 
-(Φ0, W0, B0, A0) = solveMesonWfs(ρ_sp, ρ_vp, ρ_sn, ρ_vn, r_max, divs, 3)
+(Φ0, W0, B0, A0) = solveMesonWfs(ρ_sp, ρ_vp, ρ_sn, ρ_vn, r_max, divs, 1)
 
 plot(xs, hcat(Φ0, W0, B0, A0), layout=4, label=["Φ0" "W0" "B0" "A0"])
 xlabel!("r (fm)")