69 lines
2.3 KiB
Julia
69 lines
2.3 KiB
Julia
using PolyLog, Plots
|
||
include("nucleons.jl")
|
||
include("mesons.jl")
|
||
|
||
"Total binding energy of the system"
|
||
function total_E(s::system)
|
||
E_cm = (3/4) * 41.0 * A(s)^(-1/3) # Center-of-mass correction [Ring and Schuck, Sec. 2.3]
|
||
return -(nucleon_E(s) + meson_E(s)) + E_cm
|
||
end
|
||
|
||
"Normalized Woods-Saxon form used for constructing an initial solution"
|
||
Woods_Saxon(r::Float64; R::Float64=7.0, a::Float64=0.5) = -1 / (8π * a^3 * reli3(-exp(R / a)) * (1 + exp((r - R) / a)))
|
||
|
||
"Conditionally toggle @time if enable_time is true"
|
||
macro conditional_time(label, expr)
|
||
quote
|
||
if $(esc(:enable_time))
|
||
@time $label $(esc(expr))
|
||
else
|
||
$(esc(expr))
|
||
end
|
||
end
|
||
end
|
||
|
||
"Run the full program by self-consistent solution of nucleon and meson densities"
|
||
function solve_system!(s::system; reinitialize_densities=true, print_E=true, print_time=false, live_plots=false)
|
||
if reinitialize_densities
|
||
dens_guess = Woods_Saxon.(rs(s))
|
||
@. s.ρ_sp = s.Z * dens_guess
|
||
@. s.ρ_vp = s.Z * dens_guess
|
||
@. s.ρ_sn = s.N * dens_guess
|
||
@. s.ρ_vn = s.N * dens_guess
|
||
end
|
||
|
||
if live_plots
|
||
p = plot(legends=false, size=(1024, 768), layout=(2, 4), title=["ρₛₚ" "ρᵥₚ" "ρₛₙ" "ρᵥₙ" "Φ₀" "W₀" "B₀" "A₀"])
|
||
end
|
||
|
||
previous_E_per_A = NaN
|
||
|
||
while true
|
||
enable_time = print_time
|
||
@conditional_time "Meson fields" solveMesonFields!(s, isnan(previous_E_per_A) ? 50 : 15)
|
||
@conditional_time "Proton densities" solveNucleonDensity!(true, s)
|
||
@conditional_time "Neutron densities" solveNucleonDensity!(false, s)
|
||
|
||
if live_plots
|
||
for s in p.series_list
|
||
s.plotattributes[:linecolor] = :gray
|
||
end
|
||
plot!(p, rs(s), hcat(s.ρ_sp, s.ρ_vp, s.ρ_sn, s.ρ_vn, s.Φ0, s.W0, s.B0, s.A0), linecolor=:red)
|
||
display(p)
|
||
end
|
||
|
||
E_per_A = total_E(s) / A(s)
|
||
print_E && println("Total binding E per nucleon = $E_per_A")
|
||
|
||
# check convergence
|
||
abs(previous_E_per_A - E_per_A) < 0.0001 && break
|
||
previous_E_per_A = E_per_A
|
||
end
|
||
end
|
||
|
||
"Calculate RMS radius from density"
|
||
rms_radius(p::Bool, s::system) = 4pi * Δr(s) * sum((rs(s) .^ 4) .* (p ? s.ρ_vp : s.ρ_vn)) / (p ? s.Z : s.N) |> sqrt
|
||
|
||
"Calculate neutron skin thickness"
|
||
R_skin(s::system) = rms_radius(false, s) - rms_radius(true, s)
|