yapanuwan
/
NuclearRMF
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Monitor via plots

This commit is contained in:
Nuwan Yapa 2025-02-10 16:40:22 -05:00
parent 31457c9dc7
commit 664e38d624
2 changed files with 6 additions and 4 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
system.jl
View File

@ -1,4 +1,4 @@
using Interpolations, PolyLog using Interpolations, PolyLog, Plots
include("nucleons.jl") include("nucleons.jl")
include("mesons.jl") include("mesons.jl")
@ -26,7 +26,7 @@ zero_array(s::system) = zeros(1 + s.divs)
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))) 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)))
"Run the full program by self-consistent solution of nucleon and meson densities" "Run the full program by self-consistent solution of nucleon and meson densities"
function solve_system(s::system, initial_dens=nothing, initial_flds=(zeros(1 + s.divs) for _ in 1:4)) function solve_system(s::system, initial_dens=nothing, initial_flds=(zeros(1 + s.divs) for _ in 1:4); monitor_print=true, monitor_plot=false)
if isnothing(initial_dens) if isnothing(initial_dens)
dens_guess = Woods_Saxon.(rs(s)) dens_guess = Woods_Saxon.(rs(s))
ρ_sp = s.Z .* dens_guess ρ_sp = s.Z .* dens_guess
@ -59,9 +59,11 @@ function solve_system(s::system, initial_dens=nothing, initial_flds=(zeros(1 + s
occs_n = fillNucleons(s.N, κs_n, Es_n) occs_n = fillNucleons(s.N, κs_n, Es_n)
(ρ_sp, ρ_vp) = calculateNucleonDensity(κs_n, Es_n, occs_n, false, S_interp, V_interp, R_interp, A_interp, s.r_max, s.divs) (ρ_sp, ρ_vp) = calculateNucleonDensity(κs_n, Es_n, occs_n, false, S_interp, V_interp, R_interp, A_interp, s.r_max, s.divs)
monitor_plot && plot(rs(s), hcat(ρ_sp, ρ_vp, ρ_sn, ρ_vn), layout=(2, 2), title=["ρ_sp" "ρ_vp" "ρ_sn" "ρ_vn"]) |> display
# total binding energy of nucleons # total binding energy of nucleons
E_total = sum(M_p .- Es_p) + sum(M_n .- Es_n) E_total = sum(M_p .- Es_p) + sum(M_n .- Es_n)
println("Total binding E per nucleon = $(E_total/A(s))") monitor_print && println("Total binding E per nucleon = $(E_total/A(s))")
# check convergence # check convergence
abs(E_total_previous - E_total) < 0.1 && break abs(E_total_previous - E_total) < 0.1 && break

2
test/Pb208.jl
View File

@ -2,4 +2,4 @@ include("../system.jl")
s = system(82, 126, 20.0, 400) s = system(82, 126, 20.0, 400)
solve_system(s) solve_system(s; monitor_plot=true)