Compare commits
No commits in common. "main" and "simpsons" have entirely different histories.
|
|
@ -3,27 +3,13 @@ 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
|
||||
total_E(s::system) = -(nucleon_E(s) + meson_E(s))
|
||||
|
||||
"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)
|
||||
function solve_system!(s::system; reinitialize_densities=true, monitor_print=true, monitor_plot=false)
|
||||
if reinitialize_densities
|
||||
dens_guess = Woods_Saxon.(rs(s))
|
||||
@. s.ρ_sp = s.Z * dens_guess
|
||||
|
|
@ -32,19 +18,18 @@ function solve_system!(s::system; reinitialize_densities=true, print_E=true, pri
|
|||
@. s.ρ_vn = s.N * dens_guess
|
||||
end
|
||||
|
||||
if live_plots
|
||||
if monitor_plot
|
||||
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)
|
||||
@time "Meson fields" solveMesonFields!(s, isnan(previous_E_per_A) ? 50 : 15)
|
||||
@time "Proton densities" solveNucleonDensity!(true, s)
|
||||
@time "Neutron densities" solveNucleonDensity!(false, s)
|
||||
|
||||
if live_plots
|
||||
if monitor_plot
|
||||
for s in p.series_list
|
||||
s.plotattributes[:linecolor] = :gray
|
||||
end
|
||||
|
|
@ -53,16 +38,10 @@ function solve_system!(s::system; reinitialize_densities=true, print_E=true, pri
|
|||
end
|
||||
|
||||
E_per_A = total_E(s) / A(s)
|
||||
print_E && println("Total binding E per nucleon = $E_per_A")
|
||||
monitor_print && 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)
|
||||
|
|
|
|||
15
common.jl
15
common.jl
|
|
@ -1,17 +1,2 @@
|
|||
const ħc = 197.33 # MeVfm
|
||||
const r_reg = 1E-8 # fm # regulator for R
|
||||
|
||||
"Integrate a uniformly discretized function f using Simpson's rule where h is the step size and coefficient is an optional scaling factor."
|
||||
function simpsons_integrate(f::AbstractVector{Float64}, h::Float64; coefficient::Float64 = 1.0)
|
||||
@assert length(f) % 2 == 1 "Number of mesh divisions must be even for Simpson's rule"
|
||||
s = sum(enumerate(f)) do (i, fi)
|
||||
if i == 1 || i == length(f)
|
||||
return fi
|
||||
elseif i % 2 == 0
|
||||
return 4fi
|
||||
else
|
||||
return 2fi
|
||||
end
|
||||
end
|
||||
return (h/3) * coefficient * s
|
||||
end
|
||||
|
|
|
|||
26
mesons.jl
26
mesons.jl
|
|
@ -1,6 +1,21 @@
|
|||
include("common.jl")
|
||||
include("system.jl")
|
||||
|
||||
# Values defined in C. J. Horowitz and J. Piekarewicz, Phys. Rev. Lett. 86, 5647 (2001)
|
||||
# Values taken from Hartree.f (FSUGarnet)
|
||||
const m_s = 496.939473213388 # MeV/c2
|
||||
const m_ω = 782.5 # MeV/c2
|
||||
const m_ρ = 763.0 # MeV/c2
|
||||
const m_γ = 0.000001000 # MeV/c2 # should be 0?
|
||||
const g2_s = 110.349189097820 # dimensionless
|
||||
const g2_v = 187.694676506801 # dimensionless
|
||||
const g2_ρ = 192.927428365698 # dimensionless
|
||||
const g2_γ = 0.091701236 # dimensionless # equal to 4πα
|
||||
const κ_ss = 3.260178893447 # MeV
|
||||
const λ = -0.003551486718 # dimensionless # LambdaSS
|
||||
const ζ = 0.023499504053 # dimensionless # LambdaVV
|
||||
const Λv = 0.043376933644 # dimensionless # LambdaVR
|
||||
|
||||
"Green's function for Klein-Gordon equation in natural units"
|
||||
greensFunctionKG(m, r, rp) = -1 / (m * (r + r_reg) * (rp + r_reg)) * sinh(m * min(r, rp)) * exp(-m * max(r, rp))
|
||||
|
||||
|
|
@ -27,7 +42,6 @@ end
|
|||
the inital solutions are read from s and the final solutions are saved in-place.
|
||||
Reference: P. Giuliani, K. Godbey, E. Bonilla, F. Viens, and J. Piekarewicz, Frontiers in Physics 10, (2023)"
|
||||
function solveMesonFields!(s::system, iterations=15, oscillation_control_parameter=0.3)
|
||||
(m_s, m_ω, m_ρ, m_γ, g2_s, g2_v, g2_ρ, g2_γ, κ_ss, λ, ζ, Λv) = (s.param.m_s, s.param.m_ω, s.param.m_ρ, s.param.m_γ, s.param.g2_s, s.param.g2_v, s.param.g2_ρ, s.param.g2_γ, s.param.κ_ss, s.param.λ, s.param.ζ, s.param.Λv)
|
||||
(Φ0, W0, B0, A0) = (s.Φ0, s.W0, s.B0, s.A0)
|
||||
(ρ_sp, ρ_vp, ρ_sn, ρ_vn) = (s.ρ_sp, s.ρ_vp, s.ρ_sn, s.ρ_vn)
|
||||
|
||||
|
|
@ -57,16 +71,14 @@ end
|
|||
|
||||
"Calculate the total energy associated with meson fields"
|
||||
function meson_E(s::system)
|
||||
(κ_ss, λ, ζ, Λv) = (s.param.κ_ss, s.param.λ, s.param.ζ, s.param.Λv)
|
||||
|
||||
E_densities = map(zip(s.Φ0, s.W0, s.B0, s.A0, s.ρ_sp, s.ρ_vp, s.ρ_sn, s.ρ_vn)) do (Φ0, W0, B0, A0, ρ_sp, ρ_vp, ρ_sn, ρ_vn)
|
||||
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/ħc) * (ρ_sp + ρ_sn) - ((κ_ss/ħc)/12 * (Φ0/ħc)^3 + (λ/24) * (Φ0/ħc)^4)
|
||||
E_ω = -(1/2) * (W0/ħc) * (ρ_vp + ρ_vn) + (ζ/24) * (W0/ħc)^4
|
||||
E_ρ = -(1/4) * (2B0/ħc) * (ρ_vp - ρ_vn)
|
||||
E_γ = -(1/2) * (A0/ħc) * ρ_vp
|
||||
E_ωρ = Λv * (W0/ħc)^2 * (2B0/ħc)^2
|
||||
E_σ + E_ω + E_ρ + E_γ + E_ωρ
|
||||
int += (E_σ + E_ω + E_ρ + E_γ + E_ωρ) * r^2
|
||||
end
|
||||
|
||||
return simpsons_integrate(E_densities .* rs(s).^2, Δr(s); coefficient=4π * ħc)
|
||||
return 4π * int * Δr(s) * ħc
|
||||
end
|
||||
|
|
|
|||
|
|
@ -89,9 +89,7 @@ function solveNucleonWf(κ, p::Bool, E, s::system; normalize=true, algo=Vern9())
|
|||
wf = hcat(wf_left[:, 1:(end - 1)], wf_right)
|
||||
|
||||
if normalize
|
||||
g2_int = simpsons_integrate(wf[1, :] .^ 2, Δr(s))
|
||||
f2_int = simpsons_integrate(wf[2, :] .^ 2, Δr(s))
|
||||
wf ./= sqrt(g2_int + f2_int)
|
||||
wf ./= norm(wf) * sqrt(Δr(s)) # integration by Reimann sum
|
||||
end
|
||||
|
||||
return wf
|
||||
|
|
|
|||
34
system.jl
34
system.jl
|
|
@ -1,30 +1,3 @@
|
|||
"Stores a set of coupling constants and meson masses that goes into the Lagrangian"
|
||||
struct parameters
|
||||
# Values defined in C. J. Horowitz and J. Piekarewicz, Phys. Rev. Lett. 86, 5647 (2001)
|
||||
m_s::Float64 # MeV/c2
|
||||
m_ω::Float64 # MeV/c2
|
||||
m_ρ::Float64 # MeV/c2
|
||||
m_γ::Float64 # MeV/c2
|
||||
g2_s::Float64 # dimensionless
|
||||
g2_v::Float64 # dimensionless
|
||||
g2_ρ::Float64 # dimensionless
|
||||
g2_γ::Float64 # dimensionless
|
||||
κ_ss::Float64 # MeV
|
||||
λ::Float64 # dimensionless, aka LambdaSS
|
||||
ζ::Float64 # dimensionless, aka LambdaVV
|
||||
Λv::Float64 # dimensionless, aka LambdaVR
|
||||
|
||||
"Dummy struct when parameters are not needed (for testing)"
|
||||
parameters() = new(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)
|
||||
|
||||
"Initialize parameters from a string with values provided in order of struct definition separated by commas"
|
||||
function parameters(s::String)
|
||||
values = parse.(Float64, strip.(split(s, ',')))
|
||||
@assert length(values) == 12 "String must contain exactly 12 values separated by commas"
|
||||
return new(values...)
|
||||
end
|
||||
end
|
||||
|
||||
"Tabulates a nucleon spectrum (protons or neutrons) containing κ and occupancy"
|
||||
struct spectrum
|
||||
κ::Vector{Int}
|
||||
|
|
@ -40,7 +13,6 @@ mutable struct system
|
|||
Z::Int
|
||||
N::Int
|
||||
|
||||
param::parameters
|
||||
r_max::Float64
|
||||
divs::Int
|
||||
|
||||
|
|
@ -58,10 +30,10 @@ mutable struct system
|
|||
ρ_vn::Vector{Float64}
|
||||
|
||||
"Initialize an unsolved system"
|
||||
system(Z, N, parameters, r_max, divs) = new(Z, N, parameters, r_max, divs, unfilled_spectrum(), unfilled_spectrum(), [zeros(1 + divs) for _ in 1:8]...)
|
||||
system(Z, N, r_max, divs) = new(Z, N, r_max, divs, unfilled_spectrum(), unfilled_spectrum(), [zeros(1 + divs) for _ in 1:8]...)
|
||||
|
||||
"Dummy struct to define the mesh (for testing)"
|
||||
system(r_max, divs) = system(0, 0, parameters(), r_max, divs)
|
||||
"Dummy struct to define the mesh"
|
||||
system(r_max, divs) = system(0, 0, r_max, divs)
|
||||
end
|
||||
|
||||
"Get mass number of nucleus"
|
||||
|
|
|
|||
|
|
@ -13,7 +13,7 @@ As = test_data[:, 5]
|
|||
p = false
|
||||
r_max = maximum(xs)
|
||||
divs = length(xs) - 1
|
||||
s = system(8, 8, parameters(), r_max, divs)
|
||||
s = system(8, 8, r_max, divs)
|
||||
|
||||
s.Φ0 = Ss
|
||||
s.W0 = Vs
|
||||
|
|
|
|||
|
|
@ -2,13 +2,9 @@ include("../NuclearRMF.jl")
|
|||
|
||||
Z = 82
|
||||
N = 126
|
||||
|
||||
# Parameter values calibrated by FSUGarnet for Pb-208
|
||||
param = parameters("496.939473213388, 782.5, 763.0, 0.000001000, 110.349189097820, 187.694676506801, 192.927428365698, 0.091701236, 3.260178893447, -0.003551486718, 0.023499504053, 0.043376933644")
|
||||
|
||||
r_max = 20.0
|
||||
divs = 400
|
||||
|
||||
s = system(Z, N, param, r_max, divs)
|
||||
s = system(Z, N, r_max, divs)
|
||||
|
||||
solve_system!(s; live_plots=false)
|
||||
solve_system!(s; monitor_plot=true)
|
||||
|
|
|
|||
|
|
@ -1 +0,0 @@
|
|||
496.939473213388, 782.5, 763.0, 0.000001000, 110.349189097820, 187.694676506801, 192.927428365698, 0.091701236, 3.260178893447, -0.003551486718, 0.023499504053, 0.043376933644
|
||||
|
|
@ -17,12 +17,9 @@ plot(xs_bench, hcat(Φ0_bench, W0_bench, B0_bench, A0_bench), layout=4, label=["
|
|||
test_data = readdlm("test/Pb208DensFSUGarnet.csv")
|
||||
xs = test_data[:, 1]
|
||||
|
||||
N_p = 82
|
||||
N_n = 126
|
||||
param = parameters(read("test/Pb208ParamsFSUGarnet.txt", String))
|
||||
r_max = maximum(xs)
|
||||
divs = length(xs) - 1
|
||||
s = system(N_p, N_n, param, r_max, divs)
|
||||
s = system(r_max, divs)
|
||||
|
||||
s.ρ_sn = test_data[:, 2]
|
||||
s.ρ_vn = test_data[:, 3]
|
||||
|
|
|
|||
|
|
@ -14,7 +14,7 @@ N_p = 82
|
|||
N_n = 126
|
||||
r_max = maximum(xs)
|
||||
divs = length(xs) - 1
|
||||
s = system(N_p, N_n, parameters(), r_max, divs)
|
||||
s = system(N_p, N_n, r_max, divs)
|
||||
|
||||
s.Φ0 = Ss
|
||||
s.W0 = Vs
|
||||
|
|
|
|||
|
|
@ -15,7 +15,7 @@ N_p = 82
|
|||
N_n = 126
|
||||
r_max = maximum(xs)
|
||||
divs = length(xs) - 1
|
||||
s = system(N_p, N_n, parameters(), r_max, divs)
|
||||
s = system(N_p, N_n, r_max, divs)
|
||||
|
||||
s.Φ0 = Ss
|
||||
s.W0 = Vs
|
||||
|
|
|
|||
Loading…
Reference in New Issue