Total E calculation (meson part wrong)

mesons.jl

@@ -62,3 +62,17 @@ function solveMesonFields!(s::system, iterations=10)
 
     return (Φ0, W0, B0, A0)
 end
+
+"Calculate the total energy associated with meson fields"
+function meson_E(s::system)
+    int = 0.0
+    for (r, Φ0, W0, B0, A0, ρ_sp, ρ_vp, ρ_sn, ρ_vn) in zip(rs(s), s.Φ0, s.W0, s.B0, s.A0, s.ρ_sp, s.ρ_vp, s.ρ_sn, s.ρ_vn)
+        E_σ = (1/2) * Φ0 * (ρ_sp + ρ_sn) - ((κ_ss/ħc)/12 * (Φ0/ħc)^3 + (λ/24) * (Φ0/ħc)^4)
+        E_ω = -(1/2) * W0 * (ρ_vp + ρ_vn) + (ζ/24) * (W0/ħc)^4
+        E_ρ = -(1/4) * (2B0/ħc) * (ρ_vp - ρ_vn)
+        E_γ = -(1/2) * g2_γ * A0 * ρ_vp
+        E_ωρ = Λv * (W0/ħc)^2 * (2B0/ħc)^2
+        int += (E_σ + E_ω + E_ρ + E_γ + E_ωρ) * r^2
+    end
+    return 4π * int * Δr(s) * ħc
+end

system.jl

@@ -69,11 +69,11 @@ function solve_system(s::system; reinitialize_densities=true, monitor_print=true
         p = plot(legends=false, size=(1024, 768), layout=(2, 4), title=["ρₛₚ" "ρᵥₚ" "ρₛₙ" "ρᵥₙ" "Φ₀" "W₀" "B₀" "A₀"])
     end
 
-    E_total_previous = NaN
+    previous_E_per_A = NaN
 
     while true
         # mesons
-        @time "Meson fields" solveMesonFields!(s, isnan(E_total_previous) ? 50 : 5)
+        @time "Meson fields" solveMesonFields!(s, isnan(previous_E_per_A) ? 50 : 5)
 
         # protons
         @time "Proton spectrum" (κs_p, Es_p) = findAllOrbitals(true, s)
@@ -93,12 +93,14 @@ function solve_system(s::system; reinitialize_densities=true, monitor_print=true
             display(p)
         end
 
-        # total binding energy of nucleons
-        E_total = sum(M_p .- Es_p) + sum(M_n .- Es_n)
-        monitor_print && println("Total binding E per nucleon = $(E_total/A(s))")
+        E_per_A = total_E(s, Es_p, occs_p, Es_n, occs_n) / A(s)
+        monitor_print && println("Total binding E per nucleon = $E_per_A")
 
         # check convergence
-        abs(E_total_previous - E_total) < 0.1 && break
-        E_total_previous = E_total
+        abs(previous_E_per_A - E_per_A) < 0.1 && break
+        previous_E_per_A = E_per_A
     end
 end
+
+"Total binding energy of the system"
+total_E(s::system, Es_p, occs_p, Es_n, occs_n) = sum(occs_p .* (M_p .- Es_p)) + sum(occs_n .* (M_n .- Es_n)) - meson_E(s)