Compare commits
No commits in common. "main" and "structs" have entirely different histories.
|
|
@ -1,68 +0,0 @@
|
|||
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)
|
||||
|
|
@ -3,3 +3,4 @@ DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
|
|||
Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
|
||||
Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
|
||||
PolyLog = "85e3b03c-9856-11eb-0374-4dc1f8670e7f"
|
||||
Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
|
||||
|
|
|
|||
62
bisection.jl
62
bisection.jl
|
|
@ -1,62 +0,0 @@
|
|||
"Bisection method: Finds a zero in the interval [a, b]"
|
||||
function bisection(f::Function, a::Float64, b::Float64; tol::Float64 = 1e-8, max_iter::Int = 1000)::Float64
|
||||
fa::Float64 = f(a)
|
||||
fb::Float64 = f(b)
|
||||
|
||||
# Ensure the endpoints bracket a root.
|
||||
if fa * fb > 0.0
|
||||
error("No sign change detected in the interval: f(a)*f(b) > 0")
|
||||
end
|
||||
|
||||
c::Float64 = a
|
||||
for _ in 1:max_iter
|
||||
c = (a + b) / 2.0
|
||||
fc::Float64 = f(c)
|
||||
|
||||
# Check for convergence using the function value or the interval width.
|
||||
if abs(fc) < tol || abs(b - a) < tol
|
||||
return c
|
||||
elseif fa * fc < 0.0
|
||||
b = c
|
||||
fb = fc
|
||||
else
|
||||
a = c
|
||||
fa = fc
|
||||
end
|
||||
end
|
||||
return c # Return the current approximation if max iterations reached.
|
||||
end
|
||||
|
||||
"Function to find all zeros in [x_start, x_end] by scanning for sign changes."
|
||||
function find_all_zeros(f::Function, x_start::Float64, x_end::Float64; partitions::Int = 1000, tol::Float64 = 1e-8)::Vector{Float64}
|
||||
zeros_list::Vector{Float64} = Float64[]
|
||||
Δ::Float64 = (x_end - x_start) / partitions
|
||||
x_prev::Float64 = x_start
|
||||
f_prev::Float64 = f(x_prev)
|
||||
|
||||
for i in 1:partitions
|
||||
x_curr::Float64 = x_start + i * Δ
|
||||
f_curr::Float64 = f(x_curr)
|
||||
|
||||
# Check for a sign change or an exact zero.
|
||||
if f_prev * f_curr < 0.0 || f_prev == 0.0
|
||||
try
|
||||
root::Float64 = bisection(f, x_prev, x_curr; tol=tol)
|
||||
# Add the root if it's new (avoid duplicates from neighboring intervals)
|
||||
if isempty(zeros_list) || abs(root - zeros_list[end]) > tol
|
||||
push!(zeros_list, root)
|
||||
end
|
||||
catch err
|
||||
# Skip this interval if bisection fails (should not occur if a sign change exists)
|
||||
end
|
||||
elseif f_curr == 0.0
|
||||
if isempty(zeros_list) || abs(x_curr - zeros_list[end]) > tol
|
||||
push!(zeros_list, x_curr)
|
||||
end
|
||||
end
|
||||
|
||||
x_prev = x_curr
|
||||
f_prev = f_curr
|
||||
end
|
||||
return zeros_list
|
||||
end
|
||||
17
common.jl
17
common.jl
|
|
@ -1,17 +0,0 @@
|
|||
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
|
||||
56
mesons.jl
56
mesons.jl
|
|
@ -1,5 +1,23 @@
|
|||
include("common.jl")
|
||||
include("system.jl")
|
||||
using DifferentialEquations
|
||||
|
||||
const ħc = 197.33 # MeVfm
|
||||
|
||||
# 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
|
||||
|
||||
const r_reg = 1E-8 # fm # regulator for Green's functions
|
||||
|
||||
"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))
|
||||
|
|
@ -10,10 +28,7 @@ greensFunctionP(r, rp) = -1 / (max(r, rp) + r_reg)
|
|||
"Solve the Klein-Gordon equation (or Poisson's equation when m=0) and return an array in MeV for a source array given in fm⁻³ where
|
||||
m is the mass of the meson in MeV/c2."
|
||||
function solveKG(m, source, s::system)
|
||||
|
||||
@assert s.divs % 2 == 0 "Number of mesh divisions must be even for Simpson's rule"
|
||||
simpsons_weights = (Δr(s)/3) .* [1; repeat([2,4], s.divs ÷ 2)[2:end]; 1]
|
||||
int_measure = ħc .* simpsons_weights .* rs(s) .^ 2
|
||||
int_measure = ħc .* Δr(s) .* rs(s) .^ 2
|
||||
|
||||
greensFunction = m == 0 ? greensFunctionP : (r, rp) -> greensFunctionKG(m / ħc, r, rp)
|
||||
greenMat = greensFunction.(rs(s), transpose(rs(s)))
|
||||
|
|
@ -26,8 +41,7 @@ end
|
|||
r is the radius in fm,
|
||||
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)
|
||||
function solveMesonFields!(s::system, iterations=10)
|
||||
(Φ0, W0, B0, A0) = (s.Φ0, s.W0, s.B0, s.A0)
|
||||
(ρ_sp, ρ_vp, ρ_sn, ρ_vn) = (s.ρ_sp, s.ρ_vp, s.ρ_sn, s.ρ_vn)
|
||||
|
||||
|
|
@ -41,32 +55,10 @@ function solveMesonFields!(s::system, iterations=15, oscillation_control_paramet
|
|||
@. src_Φ0 = g2_s * ((κ_ss/ħc)/2 * (Φ0/ħc)^2 + (λ/6) * (Φ0/ħc)^3) - g2_s * (ρ_sp + ρ_sn)
|
||||
@. src_W0 = g2_v * ((ζ/6) * (W0/ħc)^3 + 2Λv * (2B0/ħc)^2 * (W0/ħc)) - g2_v * (ρ_vp + ρ_vn)
|
||||
@. src_B0 = (2Λv * g2_ρ * (W0/ħc)^2 * (2B0/ħc) - g2_ρ/2 * (ρ_vp - ρ_vn)) / 2
|
||||
|
||||
Φ0 .= solveKG(m_s, src_Φ0, s)
|
||||
|
||||
# W0 and B0 keep a fraction of their previous solutions to suppress oscillations
|
||||
W0 .*= (1 - oscillation_control_parameter)
|
||||
B0 .*= (1 - oscillation_control_parameter)
|
||||
|
||||
W0 .+= solveKG(m_ω, src_W0, s) .* oscillation_control_parameter
|
||||
B0 .+= solveKG(m_ρ, src_B0, s) .* oscillation_control_parameter
|
||||
W0 .= (solveKG(m_ω, src_W0, s) .+ W0) ./ 2 # to supress oscillation
|
||||
B0 .= (solveKG(m_ρ, src_B0, s) .+ B0) ./ 2 # to supress oscillation
|
||||
end
|
||||
|
||||
return (Φ0, W0, B0, A0)
|
||||
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)
|
||||
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_ωρ
|
||||
end
|
||||
|
||||
return simpsons_integrate(E_densities .* rs(s).^2, Δr(s); coefficient=4π * ħc)
|
||||
end
|
||||
|
|
|
|||
175
nucleons.jl
175
nucleons.jl
|
|
@ -1,129 +1,88 @@
|
|||
using LinearAlgebra, DifferentialEquations, Interpolations
|
||||
include("bisection.jl")
|
||||
include("common.jl")
|
||||
include("system.jl")
|
||||
using DifferentialEquations, Roots
|
||||
|
||||
const ħc = 197.33 # MeVfm
|
||||
const M_n = 939.0 # MeV/c2
|
||||
const M_p = 939.0 # MeV/c2
|
||||
|
||||
const r_reg = 1E-8 # fm # regulator for the centrifugal term
|
||||
|
||||
"The spherical Dirac equation that returns du=[dg, df] in-place where
|
||||
u=[g, f] are the reduced radial components evaluated at r,
|
||||
κ is the generalized angular momentum,
|
||||
p is true for proton and false for neutron,
|
||||
E in the energy in MeV,
|
||||
f1(r) = M-Φ0(r)-W0(r)-(p-0.5)*2B0(r)-p*A0(r) is a function of r in MeV (see optimized_dirac_potentials()),
|
||||
f2(r) = -M+Φ0(r)-W0(r)-(p-0.5)*2B0(r)-p*A0(r) is a function of r in MeV (see optimized_dirac_potentials()),
|
||||
Φ0, W0, B0, A0 are the mean-field potentials (couplings included) in MeV as functions of r in fm,
|
||||
r is the radius in fm.
|
||||
Reference: P. Giuliani, K. Godbey, E. Bonilla, F. Viens, and J. Piekarewicz, Frontiers in Physics 10, (2023)"
|
||||
function dirac!(du::Vector{Float64}, u::Vector{Float64}, (κ, E, f1, f2), r::Float64) # TODO: Static typing
|
||||
function dirac!(du, u, (κ, p, E, Φ0, W0, B0, A0), r)
|
||||
M = p ? M_p : M_n
|
||||
common1 = E - W0(r) - (p - 0.5) * 2B0(r) - p * A0(r)
|
||||
common2 = M - Φ0(r)
|
||||
(g, f) = u
|
||||
@inbounds du[1] = -(κ/(r + r_reg)) * g + (E + f1(r)) * f / ħc
|
||||
@inbounds du[2] = (κ/(r + r_reg)) * f - (E + f2(r)) * g / ħc
|
||||
du[1] = -(κ/(r + r_reg)) * g + (common1 + common2) * f / ħc
|
||||
du[2] = (κ/(r + r_reg)) * f - (common1 - common2) * g / ħc
|
||||
end
|
||||
|
||||
"Get the potentials f1 and f2 that goes into the Dirac equation, defined as
|
||||
f1(r) = M-Φ0(r)-W0(r)-(p-0.5)*2B0(r)-p*A0(r),
|
||||
f2(r) = -M+Φ0(r)-W0(r)-(p-0.5)*2B0(r)-p*A0(r)."
|
||||
function optimized_dirac_potentials(p::Bool, s::system)
|
||||
M = p ? M_p : M_n
|
||||
|
||||
f1s = zero_array(s)
|
||||
f2s = zero_array(s)
|
||||
|
||||
@. f1s = M - s.Φ0 - s.W0 - (p - 0.5) * 2 * s.B0 - p * s.A0
|
||||
@. f2s = -M + s.Φ0 - s.W0 - (p - 0.5) * 2 * s.B0 - p * s.A0
|
||||
|
||||
f1 = linear_interpolation(rs(s), f1s)
|
||||
f2 = linear_interpolation(rs(s), f2s)
|
||||
|
||||
return (f1, f2)
|
||||
end
|
||||
|
||||
"Approximate boundary condition for u(r)=[g(r), f(r)] at r -> ∞ where
|
||||
κ is the generalized angular momentum,
|
||||
p is true for proton and false for neutron,
|
||||
E is the energy in MeV,
|
||||
r is the radius in fm."
|
||||
function asymp_BC(κ::Int, p::Bool, E::Float64, r::Float64)
|
||||
M = p ? M_p : M_n
|
||||
g = 1.0
|
||||
f = ħc / (E + M) * (-√(M^2 - E^2) + κ/r) * g
|
||||
return [g, f]
|
||||
end
|
||||
|
||||
"Initial boundary condition for u(r)=[g(r), f(r)] at r=0"
|
||||
init_BC() = [1.0, 1.0] # TODO: Why not [0.0, 1.0]?
|
||||
|
||||
"Solve the Dirac equation and return the wave function u(r)=[g(r), f(r)] where
|
||||
divs is the number of mesh divisions so solution would be returned as a 2×(1+divs) matrix,
|
||||
shooting method divides the interval into two partitions at r_max/2, ensuring convergence at both r=0 and r=r_max,
|
||||
the other parameters are the same from dirac!(...)."
|
||||
function solveNucleonWf(κ, p::Bool, E, s::system; normalize=true, algo=Vern9())
|
||||
(f1, f2) = optimized_dirac_potentials(p, s)
|
||||
function solveNucleonWf(κ, p, E, s::system; shooting=true, normalize=true)
|
||||
(Φ0, W0, B0, A0) = fields_interp(s)
|
||||
|
||||
# partitioning
|
||||
mid_idx = s.divs ÷ 2
|
||||
r_mid = rs(s)[mid_idx]
|
||||
left_r = rs(s)[1:mid_idx]
|
||||
right_r = rs(s)[mid_idx:end]
|
||||
|
||||
# left partition
|
||||
prob = ODEProblem(dirac!, init_BC(), (0, r_mid))
|
||||
sol = solve(prob, algo, p=(κ, E, f1, f2), saveat=left_r)
|
||||
wf_left = hcat(sol.u...)
|
||||
|
||||
# right partition
|
||||
prob = ODEProblem(dirac!, asymp_BC(κ, p, E, s.r_max), (s.r_max, r_mid))
|
||||
sol = solve(prob, algo, p=(κ, E, f1, f2), saveat=right_r)
|
||||
wf_right = reverse(hcat(sol.u...); dims=2)
|
||||
|
||||
# join two segments
|
||||
u1 = wf_left[:, end]
|
||||
u2 = wf_right[:, 1]
|
||||
if norm(u2) < 1e-10
|
||||
@warn "Right partition too small to rescale, setting to zero"
|
||||
wf_right .= 0.0
|
||||
if shooting
|
||||
@assert s.divs % 2 == 0 "divs must be an even number when shooting=true"
|
||||
prob = ODEProblem(dirac!, [0, 1], (s.r_max, s.r_max / 2))
|
||||
sol = solve(prob, RK4(), p=(κ, p, E, Φ0, W0, B0, A0), saveat=Δr(s))
|
||||
wf_right = reverse(hcat(sol.u...); dims=2)
|
||||
next_r_max = s.r_max / 2 # for the next step
|
||||
else
|
||||
proj = only(u1' * u2) / norm(u2)^2
|
||||
wf_right .*= proj
|
||||
next_r_max = s.r_max
|
||||
end
|
||||
|
||||
prob = ODEProblem(dirac!, [0, 1], (0, next_r_max))
|
||||
sol = solve(prob, RK4(), p=(κ, p, E, Φ0, W0, B0, A0), saveat=Δr(s))
|
||||
wf = hcat(sol.u...)
|
||||
|
||||
if shooting # join two segments
|
||||
rescale_factor_g = wf[1, end] / wf_right[1, 1]
|
||||
rescale_factor_f = wf[2, end] / wf_right[2, 1]
|
||||
@assert isfinite(rescale_factor_g) && isfinite(rescale_factor_f) "Cannot rescale the right partition"
|
||||
isapprox(rescale_factor_g, rescale_factor_f; rtol=0.03) || @warn "Discontinuity between two partitions"
|
||||
wf_right_rescaled = wf_right .* ((rescale_factor_g + rescale_factor_f) / 2)
|
||||
wf = hcat(wf[:, 1:(end - 1)], wf_right_rescaled)
|
||||
end
|
||||
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)
|
||||
norm = sum(wf .* wf) * Δr(s) # integration by Reimann sum
|
||||
wf = wf ./ sqrt(norm)
|
||||
end
|
||||
|
||||
return wf
|
||||
end
|
||||
|
||||
"Returns a function that solves the Dirac equation in two partitions and returns
|
||||
the determinant of [g_left(r) g_right(r); f_left(r) f_right(r)],
|
||||
where is r is in fm."
|
||||
function determinantFunc(κ, p::Bool, s::system, r::Float64=s.r_max/2, algo=Vern9())
|
||||
(f1, f2) = optimized_dirac_potentials(p, s)
|
||||
prob_left = ODEProblem(dirac!, init_BC(), (0, r))
|
||||
function func(E)
|
||||
prob_right = ODEProblem(dirac!, asymp_BC(κ, p, E, s.r_max), (s.r_max, r))
|
||||
u_left = solve(prob_left, algo, p=(κ, E, f1, f2), saveat=[r])
|
||||
u_right = solve(prob_right, algo, p=(κ, E, f1, f2), saveat=[r])
|
||||
return u_left[1, 1] * u_right[2, 1] - u_right[1, 1] * u_left[2, 1]
|
||||
end
|
||||
"Returns a function that solves the Dirac equation and returns g(r=r_max) where
|
||||
r_max is the outer boundary in fm,
|
||||
the other parameters are the same from dirac!(...)."
|
||||
function boundaryValueFunc(κ, p, s::system; dtype=Float64, algo=Tsit5())
|
||||
prob = ODEProblem(dirac!, convert.(dtype, [0, 1]), (0, s.r_max))
|
||||
(Φ0, W0, B0, A0) = fields_interp(s)
|
||||
func(E) = solve(prob, algo, p=(κ, p, E, Φ0, W0, B0, A0), saveat=[s.r_max], save_idxs=[1])[1, 1]
|
||||
return func
|
||||
end
|
||||
|
||||
"Find all bound energies between E_min (=850.0) and E_max (=938.0) where
|
||||
"Find all bound energies between E_min (=0) and E_max (=mass) where
|
||||
tol_digits determines the precision for root finding and the threshold for identifying duplicates,
|
||||
the other parameters are the same from dirac!(...)."
|
||||
function findEs(κ, p::Bool, s::system, E_min=850.0, E_max=938.0; tol_digits=8)
|
||||
func = determinantFunc(κ, p, s)
|
||||
Es = find_all_zeros(func, E_min, E_max; partitions=20, tol=1/10^tol_digits)
|
||||
function findEs(κ, p, s::system, E_min=0, E_max=(p ? M_p : M_n), tol_digits=5)
|
||||
func = boundaryValueFunc(κ, p, s)
|
||||
Es = find_zeros(func, (E_min, E_max); xatol=1/10^tol_digits)
|
||||
return unique(E -> round(E; digits=tol_digits), Es)
|
||||
end
|
||||
|
||||
"Find all orbitals and return two lists containing κ values and corresponding energies for a single species where
|
||||
the other parameters are defined above"
|
||||
function findAllOrbitals(p::Bool, s::system, E_min=850.0, E_max=938.0)
|
||||
function findAllOrbitals(p, s::system, E_min=0, E_max=(p ? M_p : M_n))
|
||||
κs = Int[]
|
||||
Es = Float64[]
|
||||
# start from κ=-1 and go both up and down
|
||||
|
|
@ -138,8 +97,8 @@ function findAllOrbitals(p::Bool, s::system, E_min=850.0, E_max=938.0)
|
|||
return (κs, Es)
|
||||
end
|
||||
|
||||
"For a given list of κ values with corresponding energies, attempt to fill Z_or_N lowest lying orbitals and return the spectrum"
|
||||
function fillNucleons(Z_or_N::Int, κs, Es)::spectrum
|
||||
"For a given list of κ values with corresponding energies, attempt to fill Z_or_N lowest lying orbitals and return occupancy numbers"
|
||||
function fillNucleons(Z_or_N::Int, κs, Es)
|
||||
sort_i = sortperm(Es)
|
||||
|
||||
occ = zeros(Int, length(κs))
|
||||
|
|
@ -151,8 +110,8 @@ function fillNucleons(Z_or_N::Int, κs, Es)::spectrum
|
|||
Z_or_N -= occ[i]
|
||||
end
|
||||
|
||||
@assert Z_or_N == 0 "All orbitals could not be filled"
|
||||
return spectrum(κs, Es, occ)
|
||||
Z_or_N == 0 || @warn "All orbitals could not be filled"
|
||||
return occ
|
||||
end
|
||||
|
||||
"Total angular momentum j for a given κ value"
|
||||
|
|
@ -164,14 +123,11 @@ l_κ(κ::Int) = abs(κ) - (κ < 0) # since true = 1 and false = 0
|
|||
"Calculate scalar and vector densities of a nucleon species on [0,r_max] divided into (divs+1) points and returns them as vectors (ρ_s, ρ_v) where
|
||||
the arrays κs, Es, occs tabulate the energies and occupation numbers corresponding to each κ,
|
||||
the other parameters are defined above"
|
||||
function calculateNucleonDensity(p::Bool, s::system)::Tuple{Vector{Float64}, Vector{Float64}}
|
||||
spectrum = p ? s.p_spectrum : s.n_spectrum
|
||||
(κs, Es, occs) = (spectrum.κ, spectrum.E, spectrum.occ)
|
||||
|
||||
function calculateNucleonDensity(κs, Es, occs, p, s::system)
|
||||
ρr2 = zeros(2, s.divs + 1) # ρ×r² values
|
||||
|
||||
for (κ, E, occ) in zip(κs, Es, occs)
|
||||
wf = solveNucleonWf(κ, p, E, s; normalize=true)
|
||||
wf = solveNucleonWf(κ, p, E, s; shooting=true, normalize=true)
|
||||
wf2 = wf .* wf
|
||||
ρr2 += (occ / (4 * pi)) * wf2 # 2j+1 factor is accounted in the occupancy number
|
||||
end
|
||||
|
|
@ -186,25 +142,10 @@ function calculateNucleonDensity(p::Bool, s::system)::Tuple{Vector{Float64}, Vec
|
|||
return (ρ_s, ρ_v)
|
||||
end
|
||||
|
||||
"For a nucleon species, solve the Dirac equation and save the spectrum & densities in-place where
|
||||
"Solve the Dirac equation and calculate scalar and vector densities of a nucleon species where
|
||||
the other parameters are defined above"
|
||||
function solveNucleonDensity!(p::Bool, s::system, E_min=850.0, E_max=938.0)
|
||||
function solveNucleonDensity(p, s::system, E_min=800, E_max=939)
|
||||
κs, Es = findAllOrbitals(p, s, E_min, E_max)
|
||||
spec = fillNucleons(Z_or_N(s, p), κs, Es)
|
||||
if p
|
||||
s.p_spectrum = spec
|
||||
else
|
||||
s.n_spectrum = spec
|
||||
end
|
||||
(ρ_s, ρ_v) = calculateNucleonDensity(p, s)
|
||||
if p
|
||||
s.ρ_sp = ρ_s
|
||||
s.ρ_vp = ρ_v
|
||||
else
|
||||
s.ρ_sn = ρ_s
|
||||
s.ρ_vn = ρ_v
|
||||
end
|
||||
occs = fillNucleons(Z_or_N(s, p), κs, Es)
|
||||
return calculateNucleonDensity(κs, Es, occs, p, s)
|
||||
end
|
||||
|
||||
"Total energy of filled nucleons in the system"
|
||||
nucleon_E(s::system) = sum(s.p_spectrum.occ .* (s.p_spectrum.E .- M_p)) + sum(s.n_spectrum.occ .* (s.n_spectrum.E .- M_n))
|
||||
|
|
|
|||
110
system.jl
110
system.jl
|
|
@ -1,52 +1,13 @@
|
|||
"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}
|
||||
E::Vector{Float64}
|
||||
occ::Vector{Int}
|
||||
end
|
||||
|
||||
"Initializes an unfilled spectrum"
|
||||
unfilled_spectrum() = spectrum(Int[], Float64[], Int[])
|
||||
using Interpolations, PolyLog, Plots
|
||||
|
||||
"Defines a nuclear system containing relevant parameters and meson/nucleon densities"
|
||||
mutable struct system
|
||||
Z::Int
|
||||
N::Int
|
||||
|
||||
param::parameters
|
||||
r_max::Float64
|
||||
divs::Int
|
||||
|
||||
p_spectrum::spectrum
|
||||
n_spectrum::spectrum
|
||||
|
||||
Φ0::Vector{Float64}
|
||||
W0::Vector{Float64}
|
||||
B0::Vector{Float64}
|
||||
|
|
@ -58,10 +19,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, [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"
|
||||
|
|
@ -78,3 +39,66 @@ Z_or_N(s::system, p::Bool)::Int = p ? s.Z : s.N
|
|||
|
||||
"Create an empty array for the size of the mesh"
|
||||
zero_array(s::system) = zeros(1 + s.divs)
|
||||
|
||||
"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)))
|
||||
|
||||
"Create linear interpolations for the meson fields"
|
||||
function fields_interp(s::system)
|
||||
S_interp = linear_interpolation(rs(s), s.Φ0)
|
||||
V_interp = linear_interpolation(rs(s), s.W0)
|
||||
R_interp = linear_interpolation(rs(s), s.B0)
|
||||
A_interp = linear_interpolation(rs(s), s.A0)
|
||||
return (S_interp, V_interp, R_interp, A_interp)
|
||||
end
|
||||
|
||||
include("nucleons.jl")
|
||||
include("mesons.jl")
|
||||
|
||||
"Run the full program by self-consistent solution of nucleon and meson densities"
|
||||
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
|
||||
@. s.ρ_vp = s.Z * dens_guess
|
||||
@. s.ρ_sn = s.N * dens_guess
|
||||
@. s.ρ_vn = s.N * dens_guess
|
||||
end
|
||||
|
||||
if monitor_plot
|
||||
p = plot(legends=false, size=(1024, 768), layout=(2, 4), title=["ρₛₚ" "ρᵥₚ" "ρₛₙ" "ρᵥₙ" "Φ₀" "W₀" "B₀" "A₀"])
|
||||
end
|
||||
|
||||
E_total_previous = NaN
|
||||
|
||||
while true
|
||||
# mesons
|
||||
@time "Meson fields" solveMesonFields!(s, isnan(E_total_previous) ? 50 : 5)
|
||||
|
||||
# protons
|
||||
@time "Proton spectrum" (κs_p, Es_p) = findAllOrbitals(true, s)
|
||||
occs_p = fillNucleons(s.Z, κs_p, Es_p)
|
||||
@time "Proton densities" (s.ρ_sp, s.ρ_vp) = calculateNucleonDensity(κs_p, Es_p, occs_p, true, s)
|
||||
|
||||
# neutrons
|
||||
@time "Neutron spectrum" (κs_n, Es_n) = findAllOrbitals(false, s)
|
||||
occs_n = fillNucleons(s.N, κs_n, Es_n)
|
||||
@time "Neutron densities" (s.ρ_sn, s.ρ_vn) = calculateNucleonDensity(κs_n, Es_n, occs_n, false, s)
|
||||
|
||||
if monitor_plot
|
||||
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
|
||||
|
||||
# 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))")
|
||||
|
||||
# check convergence
|
||||
abs(E_total_previous - E_total) < 0.1 && break
|
||||
E_total_previous = E_total
|
||||
end
|
||||
end
|
||||
|
|
|
|||
|
|
@ -1,9 +1,9 @@
|
|||
using DelimitedFiles, Plots
|
||||
include("../nucleons.jl")
|
||||
using DelimitedFiles, Interpolations, Plots
|
||||
include("../system.jl")
|
||||
|
||||
# test data generated from Hartree.f
|
||||
# format: x S(x) V(x) R(x) A(x)
|
||||
test_data = readdlm("test/LinearFldsFSUGarnet.csv")
|
||||
test_data = readdlm("test/LinearFldsFSUGarnet.csv")
|
||||
xs = test_data[:, 1]
|
||||
Ss = test_data[:, 2]
|
||||
Vs = test_data[:, 3]
|
||||
|
|
@ -13,27 +13,30 @@ 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
|
||||
s.B0 = Rs
|
||||
s.A0 = As
|
||||
|
||||
κs, Es = findAllOrbitals(p, s)
|
||||
E_min = 860
|
||||
E_max = 939
|
||||
|
||||
κs, Es = findAllOrbitals(p, s, E_min, E_max)
|
||||
Ebinds = (p ? M_p : M_n) .- Es
|
||||
spec = fillNucleons(Z_or_N(s, p), κs, Es)
|
||||
occ = fillNucleons(Z_or_N(s, p), κs, Es)
|
||||
|
||||
# format: proton, kappa, filling, gnodes, fnodes, Ebind
|
||||
bench_data, _ = readdlm("test/LinearSpectrum.csv", ','; header=true)
|
||||
mask = findall(bench_data[:, 1] .== Int(p))
|
||||
bench_κs = bench_data[mask, 2]
|
||||
bench_Ebinds = bench_data[mask, 6]
|
||||
bench_filling = bench_data[mask, 3]
|
||||
bench_filling = bench_data[mask, 3]
|
||||
bench_occ = Int.(bench_filling .* (2 .* j_κ.(Int.(bench_κs)) .+ 1))
|
||||
|
||||
scatter(κs, Ebinds, label=(p ? "proton" : "neutron") * " spectrum")
|
||||
annotate!(κs .+ 0.3, Ebinds .+ 1, [(o, 9) for o in spec.occ])
|
||||
annotate!(κs .+ 0.3, Ebinds .+ 1, [(o, 9) for o in occ])
|
||||
|
||||
scatter!(bench_κs, bench_Ebinds, label="benchmark"; markershape=:cross, markercolor="red")
|
||||
annotate!(bench_κs .+ 0.3, bench_Ebinds .- 1, [(o, 9, :red) for o in bench_occ])
|
||||
|
|
|
|||
|
|
@ -1,14 +1,5 @@
|
|||
include("../NuclearRMF.jl")
|
||||
include("../system.jl")
|
||||
|
||||
Z = 82
|
||||
N = 126
|
||||
s = system(82, 126, 20.0, 400)
|
||||
|
||||
# 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)
|
||||
|
||||
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
|
||||
|
|
@ -1,9 +1,9 @@
|
|||
using DelimitedFiles, Plots
|
||||
include("../mesons.jl")
|
||||
using DelimitedFiles, Interpolations, Plots
|
||||
include("../system.jl")
|
||||
|
||||
# test data generated from Hartree.f
|
||||
# format: x S(x) V(x) R(x) A(x)
|
||||
bench_data = readdlm("test/Pb208FldsFSUGarnet.csv")
|
||||
bench_data = readdlm("test/Pb208FldsFSUGarnet.csv")
|
||||
xs_bench = bench_data[:, 1]
|
||||
Φ0_bench = bench_data[:, 2]
|
||||
W0_bench = bench_data[:, 3]
|
||||
|
|
@ -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]
|
||||
|
|
|
|||
|
|
@ -1,5 +1,5 @@
|
|||
using DelimitedFiles, Plots
|
||||
include("../nucleons.jl")
|
||||
include("../system.jl")
|
||||
|
||||
# test data generated from Hartree.f
|
||||
# format: x S(x) V(x) R(x) A(x)
|
||||
|
|
@ -21,18 +21,18 @@ s.W0 = Vs
|
|||
s.B0 = Rs
|
||||
s.A0 = As
|
||||
|
||||
E_min = 850.0
|
||||
E_max = 938.0
|
||||
E_min = 850
|
||||
E_max = 939
|
||||
|
||||
boundEs = findEs(κ, p, s, E_min, E_max)
|
||||
|
||||
binding_Es = (p ? M_p : M_n) .- boundEs
|
||||
println("binding energies = $binding_Es")
|
||||
|
||||
func = determinantFunc(κ, p, s)
|
||||
Es = collect(E_min:0.1:E_max)
|
||||
dets = [func(E)^2 for E in Es]
|
||||
func = boundaryValueFunc(κ, p, s)
|
||||
Es = collect(E_min:0.5:E_max)
|
||||
boundaryVals = [func(E)^2 for E in Es]
|
||||
|
||||
plot(Es, dets, yscale=:log10, label="determinant^2")
|
||||
plot(Es, boundaryVals, yscale=:log10, label="g(r_max)^2")
|
||||
vline!(boundEs, label="bound E")
|
||||
xlabel!("E (MeV)")
|
||||
|
|
@ -1,9 +1,9 @@
|
|||
using DelimitedFiles, Plots
|
||||
include("../nucleons.jl")
|
||||
include("../system.jl")
|
||||
|
||||
# test data generated from Hartree.f
|
||||
# format: x S(x) V(x) R(x) A(x)
|
||||
test_data = readdlm("test/Pb208FldsFSUGarnet.csv")
|
||||
test_data = readdlm("test/Pb208FldsFSUGarnet.csv")
|
||||
xs = test_data[:, 1]
|
||||
Ss = test_data[:, 2]
|
||||
Vs = test_data[:, 3]
|
||||
|
|
@ -14,22 +14,22 @@ 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
|
||||
s.B0 = Rs
|
||||
s.A0 = As
|
||||
|
||||
solveNucleonDensity!(true, s)
|
||||
(ρ_sp, ρ_vp) = solveNucleonDensity(true, s)
|
||||
|
||||
p_sp = plot(rs(s), s.ρ_sp, xlabel="r (fm)", label="ρₛₚ(r) calculated")
|
||||
p_vp = plot(rs(s), s.ρ_vp, xlabel="r (fm)", label="ρᵥₚ(r) calculated")
|
||||
p_sp = plot(rs(s), ρ_sp, xlabel="r (fm)", label="ρₛₚ(r) calculated")
|
||||
p_vp = plot(rs(s), ρ_vp, xlabel="r (fm)", label="ρᵥₚ(r) calculated")
|
||||
|
||||
solveNucleonDensity!(false, s)
|
||||
(ρ_sn, ρ_vn) = solveNucleonDensity(false, s)
|
||||
|
||||
p_sn = plot(rs(s), s.ρ_sn, xlabel="r (fm)", label="ρₛₙ(r) calculated")
|
||||
p_vn = plot(rs(s), s.ρ_vn, xlabel="r (fm)", label="ρᵥₙ(r) calculated")
|
||||
p_sn = plot(rs(s), ρ_sn, xlabel="r (fm)", label="ρₛₙ(r) calculated")
|
||||
p_vn = plot(rs(s), ρ_vn, xlabel="r (fm)", label="ρᵥₙ(r) calculated")
|
||||
|
||||
# benchmark data generated from Hartree.f
|
||||
# format: x Rhos(n) Rhov(n) Rhot(n) Rhos(p) Rhov(p) Rhot(p)
|
||||
|
|
|
|||
|
|
@ -1,5 +1,5 @@
|
|||
using DelimitedFiles, Plots
|
||||
include("../nucleons.jl")
|
||||
using DelimitedFiles, Interpolations, Plots
|
||||
include("../system.jl")
|
||||
|
||||
# test data generated from Hartree.f
|
||||
# format: x S(x) V(x) R(x) A(x)
|
||||
|
|
@ -15,16 +15,19 @@ 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
|
||||
s.B0 = Rs
|
||||
s.A0 = As
|
||||
|
||||
κs, Es = findAllOrbitals(p, s)
|
||||
E_min = 800
|
||||
E_max = 939
|
||||
|
||||
κs, Es = findAllOrbitals(p, s, E_min, E_max)
|
||||
Ebinds = (p ? M_p : M_n) .- Es
|
||||
spec = fillNucleons(Z_or_N(s, p), κs, Es)
|
||||
occ = fillNucleons(Z_or_N(s, p), κs, Es)
|
||||
|
||||
# format: proton, kappa, filling, gnodes, fnodes, Ebind
|
||||
bench_data, _ = readdlm("test/Pb208Spectrum.csv", ','; header=true)
|
||||
|
|
@ -35,7 +38,7 @@ bench_filling = bench_data[mask, 3]
|
|||
bench_occ = Int.(bench_filling .* (2 .* j_κ.(Int.(bench_κs)) .+ 1))
|
||||
|
||||
scatter(κs, Ebinds, label=(p ? "proton" : "neutron") * " spectrum")
|
||||
annotate!(κs .+ 0.3, Ebinds .+ 1, [(o, 9) for o in spec.occ])
|
||||
annotate!(κs .+ 0.3, Ebinds .+ 1, [(o, 9) for o in occ])
|
||||
|
||||
scatter!(bench_κs, bench_Ebinds, label="benchmark"; markershape=:cross, markercolor="red")
|
||||
annotate!(bench_κs .+ 0.3, bench_Ebinds .- 1, [(o, 9, :red) for o in bench_occ])
|
||||
|
|
|
|||
|
|
@ -1,5 +1,5 @@
|
|||
using DelimitedFiles, Plots
|
||||
include("../nucleons.jl")
|
||||
using DelimitedFiles, Interpolations, Plots
|
||||
include("../system.jl")
|
||||
|
||||
# test data generated from Hartree.f
|
||||
# format: x S(x) V(x) R(x) A(x)
|
||||
|
|
@ -21,7 +21,10 @@ s.W0 = Vs
|
|||
s.B0 = Rs
|
||||
s.A0 = As
|
||||
|
||||
groundE = findEs(κ, p, s) |> minimum
|
||||
E_min = 800
|
||||
E_max = 939
|
||||
|
||||
groundE = findEs(κ, p, s, E_min, E_max) |> minimum
|
||||
println("ground state E = $groundE")
|
||||
|
||||
wf = solveNucleonWf(κ, p, groundE, s)
|
||||
|
|
|
|||
Loading…
Reference in New Issue