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

Saving single particle spectrum in struct system

This commit is contained in:
Nuwan Yapa 2025-02-24 20:10:28 -05:00
parent f1553f548f
commit 3a0cd0a26a
3 changed files with 22 additions and 9 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

16
NuclearRMF.jl
View File

@ -3,7 +3,7 @@ include("nucleons.jl")
include("mesons.jl") include("mesons.jl")
"Total binding energy of the system" "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) total_E(s::system) = -(nucleon_E(s) + meson_E(s))
"Normalized Woods-Saxon form used for constructing an initial solution" "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))) 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)))
@ -29,14 +29,14 @@ function solve_system(s::system; reinitialize_densities=true, monitor_print=true
@time "Meson fields" solveMesonFields!(s, isnan(previous_E_per_A) ? 50 : 15) @time "Meson fields" solveMesonFields!(s, isnan(previous_E_per_A) ? 50 : 15)
# protons # protons
@time "Proton spectrum" (κs_p, Es_p) = findAllOrbitals(true, s) @time "Proton spectrum" (s.κ_p, s.E_p) = findAllOrbitals(true, s)
occs_p = fillNucleons(s.Z, κs_p, Es_p) s.occ_p = fillNucleons(s.Z, s.κ_p, s.E_p)
@time "Proton densities" (s.ρ_sp, s.ρ_vp) = calculateNucleonDensity(κs_p, Es_p, occs_p, true, s) @time "Proton densities" (s.ρ_sp, s.ρ_vp) = calculateNucleonDensity(s.κ_p, s.E_p, s.occ_p, true, s)
# neutrons # neutrons
@time "Neutron spectrum" (κs_n, Es_n) = findAllOrbitals(false, s) @time "Neutron spectrum" (s.κ_n, s.E_n) = findAllOrbitals(false, s)
occs_n = fillNucleons(s.N, κs_n, Es_n) s.occ_n = fillNucleons(s.N, s.κ_n, s.E_n)
@time "Neutron densities" (s.ρ_sn, s.ρ_vn) = calculateNucleonDensity(κs_n, Es_n, occs_n, false, s) @time "Neutron densities" (s.ρ_sn, s.ρ_vn) = calculateNucleonDensity(s.κ_n, s.E_n, s.occ_n, false, s)
if monitor_plot if monitor_plot
for s in p.series_list for s in p.series_list
@ -46,7 +46,7 @@ function solve_system(s::system; reinitialize_densities=true, monitor_print=true
display(p) display(p)
end end
E_per_A = total_E(s, Es_p, occs_p, Es_n, occs_n) / A(s) E_per_A = total_E(s) / A(s)
monitor_print && println("Total binding E per nucleon = $E_per_A") monitor_print && println("Total binding E per nucleon = $E_per_A")
# check convergence # check convergence

3
nucleons.jl
View File

@ -170,3 +170,6 @@ function solveNucleonDensity(p::Bool, s::system, E_min=850.0, E_max=938.0)
occs = fillNucleons(Z_or_N(s, p), κs, Es) occs = fillNucleons(Z_or_N(s, p), κs, Es)
return calculateNucleonDensity(κs, Es, occs, p, s) return calculateNucleonDensity(κs, Es, occs, p, s)
end end
"Total energy of filled nucleons in the system"
nucleon_E(s::system) = sum(s.occ_p .* (s.E_p .- M_p)) + sum(s.occ_n .* (s.E_n .- M_n))

12
system.jl
View File

@ -8,6 +8,16 @@ mutable struct system
r_max::Float64 r_max::Float64
divs::Int divs::Int
# single particle energies and corresponding occupancy numbers for protons
κ_p::Vector{Int}
E_p::Vector{Float64}
occ_p::Vector{Int}
# single particle energies and corresponding occupancy numbers for neutrons
κ_n::Vector{Int}
E_n::Vector{Float64}
occ_n::Vector{Int}
Φ0::Vector{Float64} Φ0::Vector{Float64}
W0::Vector{Float64} W0::Vector{Float64}
B0::Vector{Float64} B0::Vector{Float64}
@ -19,7 +29,7 @@ mutable struct system
ρ_vn::Vector{Float64} ρ_vn::Vector{Float64}
"Initialize an unsolved system" "Initialize an unsolved system"
system(Z, N, r_max, divs) = new(Z, N, r_max, divs, [zeros(1 + divs) for _ in 1:8]...) system(Z, N, r_max, divs) = new(Z, N, r_max, divs, Int[], Float64[], Int[], Int[], Float64[], Int[], [zeros(1 + divs) for _ in 1:8]...)
"Dummy struct to define the mesh" "Dummy struct to define the mesh"
system(r_max, divs) = system(0, 0, r_max, divs) system(r_max, divs) = system(0, 0, r_max, divs)