Compare commits
No commits in common. "main" and "mesons" 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)
|
|
||||||
|
|
@ -2,4 +2,4 @@
|
||||||
DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
|
DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
|
||||||
Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
|
Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
|
||||||
Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
|
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
|
|
||||||
86
mesons.jl
86
mesons.jl
|
|
@ -1,5 +1,23 @@
|
||||||
include("common.jl")
|
using DifferentialEquations
|
||||||
include("system.jl")
|
|
||||||
|
const ħc = 197.327 # 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 κ = 3.260178893447 # MeV
|
||||||
|
const λ = -0.003551486718 # dimensionless # LambdaSS
|
||||||
|
const ζ = 0.023499504053 # dimensionless # LambdaVV
|
||||||
|
const Λv = 0.043376933644 # dimensionless # LambdaVR
|
||||||
|
|
||||||
|
const r_reg = 1E-9 # fm # regulator for Green's functions
|
||||||
|
|
||||||
"Green's function for Klein-Gordon equation in natural units"
|
"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))
|
greensFunctionKG(m, r, rp) = -1 / (m * (r + r_reg) * (rp + r_reg)) * sinh(m * min(r, rp)) * exp(-m * max(r, rp))
|
||||||
|
|
@ -8,65 +26,43 @@ greensFunctionKG(m, r, rp) = -1 / (m * (r + r_reg) * (rp + r_reg)) * sinh(m * mi
|
||||||
greensFunctionP(r, rp) = -1 / (max(r, rp) + r_reg)
|
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
|
"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."
|
m is the mass of the meson in MeV/c2,
|
||||||
function solveKG(m, source, s::system)
|
r_max is the r-cutoff in fm."
|
||||||
|
function solveKG(m, source, r_max)
|
||||||
@assert s.divs % 2 == 0 "Number of mesh divisions must be even for Simpson's rule"
|
Δr = r_max / (length(source) - 1)
|
||||||
simpsons_weights = (Δr(s)/3) .* [1; repeat([2,4], s.divs ÷ 2)[2:end]; 1]
|
rs = range(0, r_max; length=length(source))
|
||||||
int_measure = ħc .* simpsons_weights .* rs(s) .^ 2
|
int_measure = ħc .* Δr .* rs .^ 2
|
||||||
|
|
||||||
greensFunction = m == 0 ? greensFunctionP : (r, rp) -> greensFunctionKG(m / ħc, r, rp)
|
greensFunction = m == 0 ? greensFunctionP : (r, rp) -> greensFunctionKG(m / ħc, r, rp)
|
||||||
greenMat = greensFunction.(rs(s), transpose(rs(s)))
|
greenMat = greensFunction.(rs, transpose(rs))
|
||||||
|
|
||||||
return greenMat * (int_measure .* source)
|
return greenMat * (int_measure .* source)
|
||||||
end
|
end
|
||||||
|
|
||||||
"Iteratively solve meson equations and return the wave functions u(r)=[Φ0(r), W0(r), B0(r), A0(r)] where
|
"Iteratively solve meson equations and return the wave functions u(r)=[Φ0(r), W0(r), B0(r), A0(r)] where
|
||||||
divs is the number of mesh divisions so the arrays are of length (1+divs),
|
divs is the number of mesh divisions so the arrays are of length (1+divs),
|
||||||
|
ρ_sp, ρ_vp are the scalar and vector proton densities as arrays,
|
||||||
|
ρ_sn, ρ_vn are the scalar and vector neutron densities as arrays,
|
||||||
|
Φ0, W0, B0, A0 are the mean-field potentials (couplings included) in MeV as as arrays,
|
||||||
r is the radius in fm,
|
r is the radius in fm,
|
||||||
the inital solutions are read from s and the final solutions are saved in-place.
|
An initial guess initial_sol=(Φ0, W0, B0, A0) can be provided to speed up convergence (permuting!).
|
||||||
Reference: P. Giuliani, K. Godbey, E. Bonilla, F. Viens, and J. Piekarewicz, Frontiers in Physics 10, (2023)"
|
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)
|
function solveMesonWfs(ρ_sp, ρ_vp, ρ_sn, ρ_vn, r_max, divs, iterations=3; initial_sol=(zeros(1 + divs) for _ in 1:4))
|
||||||
(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) = initial_sol
|
||||||
(Φ0, W0, B0, A0) = (s.Φ0, s.W0, s.B0, s.A0)
|
(src_Φ0, src_W0, src_B0, src_A0) = (zeros(1 + divs) for _ in 1:4)
|
||||||
(ρ_sp, ρ_vp, ρ_sn, ρ_vn) = (s.ρ_sp, s.ρ_vp, s.ρ_sn, s.ρ_vn)
|
|
||||||
|
|
||||||
(src_Φ0, src_W0, src_B0, src_A0) = (zero_array(s) for _ in 1:4)
|
|
||||||
|
|
||||||
# A0 doesn't need iterations
|
# A0 doesn't need iterations
|
||||||
@. src_A0 = -g2_γ * ρ_vp
|
@. src_A0 = -g2_γ * ρ_vp
|
||||||
A0 .= solveKG(m_γ, src_A0, s)
|
A0 .= solveKG(m_γ, src_A0, r_max)
|
||||||
|
|
||||||
for _ in 1:iterations
|
for _ in 1:iterations
|
||||||
@. src_Φ0 = g2_s * ((κ_ss/ħc)/2 * (Φ0/ħc)^2 + (λ/6) * (Φ0/ħc)^3) - g2_s * (ρ_sp + ρ_sn)
|
@. src_Φ0 = g2_s * ((κ/ħ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_W0 = g2_v * ((ζ/6) * (W0/ħc)^3 + 2Λv * (B0/ħc)^2 * (W0/ħc)) - g2_v * (ρ_vp + ρ_vn)
|
||||||
@. src_B0 = (2Λv * g2_ρ * (W0/ħc)^2 * (2B0/ħc) - g2_ρ/2 * (ρ_vp - ρ_vn)) / 2
|
@. src_B0 = 2Λv * g2_ρ * (W0/ħc)^2 * (B0/ħc) - g2_ρ/4 * (ρ_vp - ρ_vn)
|
||||||
|
Φ0 .= solveKG(m_s, src_Φ0, r_max)
|
||||||
Φ0 .= solveKG(m_s, src_Φ0, s)
|
W0 .= solveKG(m_ω, src_W0, r_max)
|
||||||
|
B0 .= solveKG(m_ρ, src_B0, r_max)
|
||||||
# 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
|
|
||||||
end
|
end
|
||||||
|
|
||||||
return (Φ0, W0, B0, A0)
|
return (Φ0, W0, B0, A0)
|
||||||
end
|
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
|
|
||||||
|
|
|
||||||
198
nucleons.jl
198
nucleons.jl
|
|
@ -1,135 +1,89 @@
|
||||||
using LinearAlgebra, DifferentialEquations, Interpolations
|
using DifferentialEquations, Roots
|
||||||
include("bisection.jl")
|
|
||||||
include("common.jl")
|
|
||||||
include("system.jl")
|
|
||||||
|
|
||||||
const M_n = 939.0 # MeV/c2
|
const ħc = 197.327 # MeVfm
|
||||||
const M_p = 939.0 # MeV/c2
|
const M_n = 939.5654133 # MeV/c2
|
||||||
|
const M_p = 938.2720813 # MeV/c2
|
||||||
|
|
||||||
|
const r_reg = 1E-6 # fm # regulator for the centrifugal term
|
||||||
|
|
||||||
"The spherical Dirac equation that returns du=[dg, df] in-place where
|
"The spherical Dirac equation that returns du=[dg, df] in-place where
|
||||||
u=[g, f] are the reduced radial components evaluated at r,
|
u=[g, f] are the reduced radial components evaluated at r,
|
||||||
κ is the generalized angular momentum,
|
κ is the generalized angular momentum,
|
||||||
p is true for proton and false for neutron,
|
p is true for proton and false for neutron,
|
||||||
E in the energy in MeV,
|
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()),
|
Φ0, W0, B0, A0 are the mean-field potentials (couplings included) in MeV as functions of r in fm,
|
||||||
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()),
|
|
||||||
r is the radius 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)"
|
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) * B0(r) - p * A0(r)
|
||||||
|
common2 = M - Φ0(r)
|
||||||
(g, f) = u
|
(g, f) = u
|
||||||
@inbounds du[1] = -(κ/(r + r_reg)) * g + (E + f1(r)) * f / ħc
|
du[1] = -(κ/(r + r_reg)) * g + (common1 + common2) * f / ħc
|
||||||
@inbounds du[2] = (κ/(r + r_reg)) * f - (E + f2(r)) * g / ħc
|
du[2] = (κ/(r + r_reg)) * f - (common1 - common2) * g / ħc
|
||||||
end
|
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
|
"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,
|
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!(...)."
|
the other parameters are the same from dirac!(...)."
|
||||||
function solveNucleonWf(κ, p::Bool, E, s::system; normalize=true, algo=Vern9())
|
function solveNucleonWf(κ, p, E, Φ0, W0, B0, A0, r_max, divs; shooting=true, normalize=true)
|
||||||
(f1, f2) = optimized_dirac_potentials(p, s)
|
saveat = r_max / divs
|
||||||
|
|
||||||
# partitioning
|
if shooting
|
||||||
mid_idx = s.divs ÷ 2
|
@assert divs % 2 == 0 "divs must be an even number when shooting=true"
|
||||||
r_mid = rs(s)[mid_idx]
|
prob = ODEProblem(dirac!, [0, 1], (r_max, r_max / 2))
|
||||||
left_r = rs(s)[1:mid_idx]
|
sol = solve(prob, RK4(), p=(κ, p, E, Φ0, W0, B0, A0), saveat=saveat)
|
||||||
right_r = rs(s)[mid_idx:end]
|
wf_right = reverse(hcat(sol.u...); dims=2)
|
||||||
|
r_max = r_max / 2 # for the next step
|
||||||
# left partition
|
end
|
||||||
prob = ODEProblem(dirac!, init_BC(), (0, r_mid))
|
|
||||||
sol = solve(prob, algo, p=(κ, E, f1, f2), saveat=left_r)
|
prob = ODEProblem(dirac!, [0, 1], (0, r_max))
|
||||||
wf_left = hcat(sol.u...)
|
sol = solve(prob, RK4(), p=(κ, p, E, Φ0, W0, B0, A0), saveat=saveat)
|
||||||
|
wf = hcat(sol.u...)
|
||||||
# right partition
|
|
||||||
prob = ODEProblem(dirac!, asymp_BC(κ, p, E, s.r_max), (s.r_max, r_mid))
|
if shooting # join two segments
|
||||||
sol = solve(prob, algo, p=(κ, E, f1, f2), saveat=right_r)
|
rescale_factor_g = wf[1, end] / wf_right[1, 1]
|
||||||
wf_right = reverse(hcat(sol.u...); dims=2)
|
rescale_factor_f = wf[2, end] / wf_right[2, 1]
|
||||||
|
@assert isfinite(rescale_factor_g) && isfinite(rescale_factor_f) "Cannot rescale the right partition"
|
||||||
# join two segments
|
isapprox(rescale_factor_g, rescale_factor_f; rtol=0.03) || @warn "Discontinuity between two partitions"
|
||||||
u1 = wf_left[:, end]
|
wf_right_rescaled = wf_right .* ((rescale_factor_g + rescale_factor_f) / 2)
|
||||||
u2 = wf_right[:, 1]
|
wf = hcat(wf[:, 1:(end - 1)], wf_right_rescaled)
|
||||||
if norm(u2) < 1e-10
|
|
||||||
@warn "Right partition too small to rescale, setting to zero"
|
|
||||||
wf_right .= 0.0
|
|
||||||
else
|
|
||||||
proj = only(u1' * u2) / norm(u2)^2
|
|
||||||
wf_right .*= proj
|
|
||||||
end
|
end
|
||||||
wf = hcat(wf_left[:, 1:(end - 1)], wf_right)
|
|
||||||
|
|
||||||
if normalize
|
if normalize
|
||||||
g2_int = simpsons_integrate(wf[1, :] .^ 2, Δr(s))
|
norm = sum(wf .* wf) * r_max / divs # integration by Reimann sum
|
||||||
f2_int = simpsons_integrate(wf[2, :] .^ 2, Δr(s))
|
wf = wf ./ sqrt(norm)
|
||||||
wf ./= sqrt(g2_int + f2_int)
|
|
||||||
end
|
end
|
||||||
|
|
||||||
return wf
|
return wf
|
||||||
end
|
end
|
||||||
|
|
||||||
"Returns a function that solves the Dirac equation in two partitions and returns
|
"Solve the Dirac equation and return g(r=r_max) where
|
||||||
the determinant of [g_left(r) g_right(r); f_left(r) f_right(r)],
|
r_max is the outer boundary in fm,
|
||||||
where is r is in fm."
|
the other parameters are the same from dirac!(...)."
|
||||||
function determinantFunc(κ, p::Bool, s::system, r::Float64=s.r_max/2, algo=Vern9())
|
function boundaryValue(κ, p, E, Φ0, W0, B0, A0, r_max; dtype=Float64, algo=RK4())
|
||||||
(f1, f2) = optimized_dirac_potentials(p, s)
|
prob = ODEProblem(dirac!, convert.(dtype, [0, 1]), (0, r_max))
|
||||||
prob_left = ODEProblem(dirac!, init_BC(), (0, r))
|
sol = solve(prob, algo, p=(κ, p, E, Φ0, W0, B0, A0), saveat=[r_max], save_idxs=[1])
|
||||||
function func(E)
|
return sol[1, 1]
|
||||||
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
|
|
||||||
return func
|
|
||||||
end
|
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!(...)."
|
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)
|
function findEs(κ, p, Φ0, W0, B0, A0, r_max, E_min=0, E_max=(p ? M_p : M_n))
|
||||||
func = determinantFunc(κ, p, s)
|
f(E) = boundaryValue(κ, p, E, Φ0, W0, B0, A0, r_max)
|
||||||
Es = find_all_zeros(func, E_min, E_max; partitions=20, tol=1/10^tol_digits)
|
return find_zeros(f, (E_min, E_max))
|
||||||
return unique(E -> round(E; digits=tol_digits), Es)
|
|
||||||
end
|
end
|
||||||
|
|
||||||
"Find all orbitals and return two lists containing κ values and corresponding energies for a single species where
|
"Find all orbitals and return two lists containing κ values and corresponding energies for a single species where
|
||||||
the other parameters are defined above"
|
the other parameters are defined above"
|
||||||
function findAllOrbitals(p::Bool, s::system, E_min=850.0, E_max=938.0)
|
function findAllOrbitals(p, Φ0, W0, B0, A0, r_max, E_min=0, E_max=(p ? M_p : M_n))
|
||||||
κs = Int[]
|
κs = Int[]
|
||||||
Es = Float64[]
|
Es = Float64[]
|
||||||
# start from κ=-1 and go both up and down
|
# start from κ=-1 and go both up and down
|
||||||
for direction in [-1, 1]
|
for direction in [-1, 1]
|
||||||
for κ in direction * (1:100) # cutoff is 100
|
for κ in direction * (1:100) # cutoff is 100
|
||||||
new_Es = findEs(κ, p, s, E_min, E_max)
|
new_Es = findEs(κ, p, Φ0, W0, B0, A0, r_max, E_min, E_max)
|
||||||
if isempty(new_Es); break; end
|
if isempty(new_Es); break; end
|
||||||
append!(Es, new_Es)
|
append!(Es, new_Es)
|
||||||
append!(κs, fill(κ, length(new_Es)))
|
append!(κs, fill(κ, length(new_Es)))
|
||||||
|
|
@ -138,21 +92,21 @@ function findAllOrbitals(p::Bool, s::system, E_min=850.0, E_max=938.0)
|
||||||
return (κs, Es)
|
return (κs, Es)
|
||||||
end
|
end
|
||||||
|
|
||||||
"For a given list of κ values with corresponding energies, attempt to fill Z_or_N lowest lying orbitals and return the spectrum"
|
"For a given list of κ values with corresponding energies, attempt to fill N lowest lying orbitals and return occupancy numbers"
|
||||||
function fillNucleons(Z_or_N::Int, κs, Es)::spectrum
|
function fillNucleons(N::Int, κs, Es)
|
||||||
sort_i = sortperm(Es)
|
sort_i = sortperm(Es)
|
||||||
|
|
||||||
occ = zeros(Int, length(κs))
|
occ = zeros(Int, length(κs))
|
||||||
|
|
||||||
for i in sort_i
|
for i in sort_i
|
||||||
if Z_or_N ≤ 0; break; end;
|
if N ≤ 0; break; end;
|
||||||
max_occ = 2 * j_κ(κs[i]) + 1
|
max_occ = 2 * j_κ(κs[i]) + 1
|
||||||
occ[i] = min(max_occ, Z_or_N)
|
occ[i] = min(max_occ, N)
|
||||||
Z_or_N -= occ[i]
|
N -= occ[i]
|
||||||
end
|
end
|
||||||
|
|
||||||
@assert Z_or_N == 0 "All orbitals could not be filled"
|
N == 0 || @warn "All orbitals could not be filled"
|
||||||
return spectrum(κs, Es, occ)
|
return occ
|
||||||
end
|
end
|
||||||
|
|
||||||
"Total angular momentum j for a given κ value"
|
"Total angular momentum j for a given κ value"
|
||||||
|
|
@ -162,21 +116,20 @@ j_κ(κ::Int) = abs(κ) - 1/2
|
||||||
l_κ(κ::Int) = abs(κ) - (κ < 0) # since true = 1 and false = 0
|
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
|
"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"
|
the other parameters are defined above"
|
||||||
function calculateNucleonDensity(p::Bool, s::system)::Tuple{Vector{Float64}, Vector{Float64}}
|
function calculateNucleonDensity(N, p, Φ0, W0, B0, A0, r_max, divs, E_min=0, E_max=(p ? M_p : M_n))
|
||||||
spectrum = p ? s.p_spectrum : s.n_spectrum
|
κs, Es = findAllOrbitals(p, Φ0, W0, B0, A0, r_max, E_min, E_max)
|
||||||
(κs, Es, occs) = (spectrum.κ, spectrum.E, spectrum.occ)
|
occs = fillNucleons(N, κs, Es)
|
||||||
|
|
||||||
ρr2 = zeros(2, s.divs + 1) # ρ×r² values
|
ρr2 = zeros(2, divs + 1) # ρ×r² values
|
||||||
|
|
||||||
for (κ, E, occ) in zip(κs, Es, occs)
|
for (κ, E, occ) in zip(κs, Es, occs)
|
||||||
wf = solveNucleonWf(κ, p, E, s; normalize=true)
|
wf = solveNucleonWf(κ, p, E, Φ0, W0, B0, A0, r_max, divs; shooting=true, normalize=true)
|
||||||
wf2 = wf .* wf
|
wf2 = wf .* wf
|
||||||
ρr2 += (occ / (4 * pi)) * wf2 # 2j+1 factor is accounted in the occupancy number
|
ρr2 += (occ / (4 * pi)) * wf2 # 2j+1 factor is accounted in the occupancy number
|
||||||
end
|
end
|
||||||
|
|
||||||
r2s = rs(s).^2
|
r2s = (collect ∘ range)(0, r_max, length=divs+1).^2
|
||||||
ρ = ρr2 ./ transpose(r2s)
|
ρ = ρr2 ./ transpose(r2s)
|
||||||
ρ[:, 1] .= ρ[:, 2] # dirty fix for NaN at r=0
|
ρ[:, 1] .= ρ[:, 2] # dirty fix for NaN at r=0
|
||||||
|
|
||||||
|
|
@ -184,27 +137,4 @@ function calculateNucleonDensity(p::Bool, s::system)::Tuple{Vector{Float64}, Vec
|
||||||
ρ_v = ρ[1, :] + ρ[2, :] # g^2 + f^2
|
ρ_v = ρ[1, :] + ρ[2, :] # g^2 + f^2
|
||||||
|
|
||||||
return (ρ_s, ρ_v)
|
return (ρ_s, ρ_v)
|
||||||
end
|
end
|
||||||
|
|
||||||
"For a nucleon species, solve the Dirac equation and save the spectrum & densities in-place where
|
|
||||||
the other parameters are defined above"
|
|
||||||
function solveNucleonDensity!(p::Bool, s::system, E_min=850.0, E_max=938.0)
|
|
||||||
κ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
|
|
||||||
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))
|
|
||||||
80
system.jl
80
system.jl
|
|
@ -1,80 +0,0 @@
|
||||||
"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[])
|
|
||||||
|
|
||||||
"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}
|
|
||||||
A0::Vector{Float64}
|
|
||||||
|
|
||||||
ρ_sp::Vector{Float64}
|
|
||||||
ρ_vp::Vector{Float64}
|
|
||||||
ρ_sn::Vector{Float64}
|
|
||||||
ρ_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]...)
|
|
||||||
|
|
||||||
"Dummy struct to define the mesh (for testing)"
|
|
||||||
system(r_max, divs) = system(0, 0, parameters(), r_max, divs)
|
|
||||||
end
|
|
||||||
|
|
||||||
"Get mass number of nucleus"
|
|
||||||
A(s::system)::Int = s.Z + s.N
|
|
||||||
|
|
||||||
"Get r values in the mesh"
|
|
||||||
rs(s::system)::StepRangeLen = range(0, s.r_max, length=s.divs+1)
|
|
||||||
|
|
||||||
"Get Δr value for the mesh"
|
|
||||||
Δr(s::system)::Float64 = s.r_max / s.divs
|
|
||||||
|
|
||||||
"Get the number of protons or neutrons in the system"
|
|
||||||
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)
|
|
||||||
|
|
@ -1,401 +0,0 @@
|
||||||
0.000 0.44010E+03 0.36339E+03 -0.36327E+00 0.20168E-12
|
|
||||||
0.050 0.44022E+03 0.36351E+03 -0.36349E+00 0.56072E+01
|
|
||||||
0.100 0.44057E+03 0.36388E+03 -0.36414E+00 0.56056E+01
|
|
||||||
0.150 0.44115E+03 0.36447E+03 -0.36522E+00 0.56029E+01
|
|
||||||
0.200 0.44195E+03 0.36528E+03 -0.36669E+00 0.55991E+01
|
|
||||||
0.250 0.44294E+03 0.36630E+03 -0.36852E+00 0.55942E+01
|
|
||||||
0.300 0.44411E+03 0.36751E+03 -0.37067E+00 0.55882E+01
|
|
||||||
0.350 0.44544E+03 0.36888E+03 -0.37309E+00 0.55811E+01
|
|
||||||
0.400 0.44689E+03 0.37039E+03 -0.37571E+00 0.55729E+01
|
|
||||||
0.450 0.44843E+03 0.37201E+03 -0.37846E+00 0.55635E+01
|
|
||||||
0.500 0.45004E+03 0.37370E+03 -0.38129E+00 0.55530E+01
|
|
||||||
0.550 0.45167E+03 0.37543E+03 -0.38412E+00 0.55412E+01
|
|
||||||
0.600 0.45328E+03 0.37716E+03 -0.38686E+00 0.55283E+01
|
|
||||||
0.650 0.45484E+03 0.37886E+03 -0.38946E+00 0.55142E+01
|
|
||||||
0.700 0.45630E+03 0.38049E+03 -0.39182E+00 0.54988E+01
|
|
||||||
0.750 0.45763E+03 0.38202E+03 -0.39387E+00 0.54822E+01
|
|
||||||
0.800 0.45878E+03 0.38339E+03 -0.39555E+00 0.54643E+01
|
|
||||||
0.850 0.45971E+03 0.38458E+03 -0.39678E+00 0.54452E+01
|
|
||||||
0.900 0.46038E+03 0.38555E+03 -0.39750E+00 0.54247E+01
|
|
||||||
0.950 0.46075E+03 0.38626E+03 -0.39764E+00 0.54030E+01
|
|
||||||
1.000 0.46079E+03 0.38668E+03 -0.39716E+00 0.53799E+01
|
|
||||||
1.050 0.46046E+03 0.38678E+03 -0.39600E+00 0.53556E+01
|
|
||||||
1.100 0.45973E+03 0.38652E+03 -0.39413E+00 0.53299E+01
|
|
||||||
1.150 0.45857E+03 0.38587E+03 -0.39149E+00 0.53029E+01
|
|
||||||
1.200 0.45694E+03 0.38482E+03 -0.38808E+00 0.52745E+01
|
|
||||||
1.250 0.45484E+03 0.38334E+03 -0.38385E+00 0.52449E+01
|
|
||||||
1.300 0.45224E+03 0.38140E+03 -0.37881E+00 0.52140E+01
|
|
||||||
1.350 0.44911E+03 0.37899E+03 -0.37293E+00 0.51817E+01
|
|
||||||
1.400 0.44545E+03 0.37610E+03 -0.36622E+00 0.51482E+01
|
|
||||||
1.450 0.44125E+03 0.37272E+03 -0.35869E+00 0.51135E+01
|
|
||||||
1.500 0.43649E+03 0.36884E+03 -0.35034E+00 0.50775E+01
|
|
||||||
1.550 0.43119E+03 0.36445E+03 -0.34121E+00 0.50403E+01
|
|
||||||
1.600 0.42534E+03 0.35956E+03 -0.33130E+00 0.50020E+01
|
|
||||||
1.650 0.41895E+03 0.35417E+03 -0.32067E+00 0.49625E+01
|
|
||||||
1.700 0.41203E+03 0.34829E+03 -0.30935E+00 0.49219E+01
|
|
||||||
1.750 0.40458E+03 0.34192E+03 -0.29738E+00 0.48803E+01
|
|
||||||
1.800 0.39664E+03 0.33509E+03 -0.28482E+00 0.48377E+01
|
|
||||||
1.850 0.38822E+03 0.32782E+03 -0.27173E+00 0.47942E+01
|
|
||||||
1.900 0.37935E+03 0.32012E+03 -0.25817E+00 0.47497E+01
|
|
||||||
1.950 0.37005E+03 0.31202E+03 -0.24420E+00 0.47045E+01
|
|
||||||
2.000 0.36036E+03 0.30355E+03 -0.22991E+00 0.46584E+01
|
|
||||||
2.050 0.35031E+03 0.29474E+03 -0.21535E+00 0.46117E+01
|
|
||||||
2.100 0.33994E+03 0.28562E+03 -0.20062E+00 0.45643E+01
|
|
||||||
2.150 0.32928E+03 0.27624E+03 -0.18578E+00 0.45163E+01
|
|
||||||
2.200 0.31838E+03 0.26663E+03 -0.17092E+00 0.44678E+01
|
|
||||||
2.250 0.30728E+03 0.25683E+03 -0.15612E+00 0.44188E+01
|
|
||||||
2.300 0.29603E+03 0.24688E+03 -0.14145E+00 0.43695E+01
|
|
||||||
2.350 0.28465E+03 0.23682E+03 -0.12700E+00 0.43199E+01
|
|
||||||
2.400 0.27321E+03 0.22671E+03 -0.11283E+00 0.42701E+01
|
|
||||||
2.450 0.26175E+03 0.21656E+03 -0.99026E-01 0.42201E+01
|
|
||||||
2.500 0.25030E+03 0.20644E+03 -0.85640E-01 0.41700E+01
|
|
||||||
2.550 0.23890E+03 0.19638E+03 -0.72739E-01 0.41199E+01
|
|
||||||
2.600 0.22760E+03 0.18642E+03 -0.60376E-01 0.40698E+01
|
|
||||||
2.650 0.21644E+03 0.17659E+03 -0.48602E-01 0.40199E+01
|
|
||||||
2.700 0.20545E+03 0.16693E+03 -0.37457E-01 0.39701E+01
|
|
||||||
2.750 0.19467E+03 0.15747E+03 -0.26975E-01 0.39205E+01
|
|
||||||
2.800 0.18412E+03 0.14824E+03 -0.17184E-01 0.38712E+01
|
|
||||||
2.850 0.17383E+03 0.13927E+03 -0.81038E-02 0.38223E+01
|
|
||||||
2.900 0.16383E+03 0.13058E+03 0.25458E-03 0.37737E+01
|
|
||||||
2.950 0.15413E+03 0.12220E+03 0.78861E-02 0.37256E+01
|
|
||||||
3.000 0.14477E+03 0.11412E+03 0.14793E-01 0.36779E+01
|
|
||||||
3.050 0.13575E+03 0.10638E+03 0.20984E-01 0.36308E+01
|
|
||||||
3.100 0.12708E+03 0.98976E+02 0.26475E-01 0.35842E+01
|
|
||||||
3.150 0.11878E+03 0.91920E+02 0.31284E-01 0.35382E+01
|
|
||||||
3.200 0.11084E+03 0.85215E+02 0.35439E-01 0.34928E+01
|
|
||||||
3.250 0.10328E+03 0.78861E+02 0.38968E-01 0.34480E+01
|
|
||||||
3.300 0.96095E+02 0.72858E+02 0.41905E-01 0.34039E+01
|
|
||||||
3.350 0.89280E+02 0.67203E+02 0.44284E-01 0.33604E+01
|
|
||||||
3.400 0.82832E+02 0.61888E+02 0.46144E-01 0.33177E+01
|
|
||||||
3.450 0.76746E+02 0.56906E+02 0.47523E-01 0.32756E+01
|
|
||||||
3.500 0.71014E+02 0.52249E+02 0.48462E-01 0.32342E+01
|
|
||||||
3.550 0.65628E+02 0.47905E+02 0.49000E-01 0.31936E+01
|
|
||||||
3.600 0.60577E+02 0.43863E+02 0.49178E-01 0.31537E+01
|
|
||||||
3.650 0.55849E+02 0.40110E+02 0.49033E-01 0.31145E+01
|
|
||||||
3.700 0.51432E+02 0.36632E+02 0.48604E-01 0.30760E+01
|
|
||||||
3.750 0.47313E+02 0.33417E+02 0.47927E-01 0.30382E+01
|
|
||||||
3.800 0.43480E+02 0.30450E+02 0.47036E-01 0.30012E+01
|
|
||||||
3.850 0.39917E+02 0.27717E+02 0.45962E-01 0.29649E+01
|
|
||||||
3.900 0.36612E+02 0.25205E+02 0.44737E-01 0.29292E+01
|
|
||||||
3.950 0.33550E+02 0.22898E+02 0.43388E-01 0.28943E+01
|
|
||||||
4.000 0.30718E+02 0.20784E+02 0.41941E-01 0.28600E+01
|
|
||||||
4.050 0.28101E+02 0.18850E+02 0.40418E-01 0.28265E+01
|
|
||||||
4.100 0.25688E+02 0.17082E+02 0.38840E-01 0.27936E+01
|
|
||||||
4.150 0.23465E+02 0.15469E+02 0.37227E-01 0.27613E+01
|
|
||||||
4.200 0.21419E+02 0.13999E+02 0.35594E-01 0.27297E+01
|
|
||||||
4.250 0.19538E+02 0.12661E+02 0.33956E-01 0.26987E+01
|
|
||||||
4.300 0.17811E+02 0.11444E+02 0.32325E-01 0.26684E+01
|
|
||||||
4.350 0.16228E+02 0.10339E+02 0.30713E-01 0.26386E+01
|
|
||||||
4.400 0.14777E+02 0.93363E+01 0.29128E-01 0.26095E+01
|
|
||||||
4.450 0.13448E+02 0.84272E+01 0.27579E-01 0.25809E+01
|
|
||||||
4.500 0.12233E+02 0.76036E+01 0.26070E-01 0.25529E+01
|
|
||||||
4.550 0.11123E+02 0.68582E+01 0.24608E-01 0.25255E+01
|
|
||||||
4.600 0.10109E+02 0.61840E+01 0.23195E-01 0.24986E+01
|
|
||||||
4.650 0.91839E+01 0.55746E+01 0.21836E-01 0.24722E+01
|
|
||||||
4.700 0.83402E+01 0.50241E+01 0.20531E-01 0.24463E+01
|
|
||||||
4.750 0.75714E+01 0.45271E+01 0.19283E-01 0.24210E+01
|
|
||||||
4.800 0.68712E+01 0.40786E+01 0.18092E-01 0.23961E+01
|
|
||||||
4.850 0.62339E+01 0.36741E+01 0.16958E-01 0.23718E+01
|
|
||||||
4.900 0.56542E+01 0.33094E+01 0.15881E-01 0.23479E+01
|
|
||||||
4.950 0.51270E+01 0.29807E+01 0.14859E-01 0.23244E+01
|
|
||||||
5.000 0.46480E+01 0.26845E+01 0.13892E-01 0.23014E+01
|
|
||||||
5.050 0.42128E+01 0.24177E+01 0.12978E-01 0.22788E+01
|
|
||||||
5.100 0.38176E+01 0.21774E+01 0.12116E-01 0.22567E+01
|
|
||||||
5.150 0.34589E+01 0.19611E+01 0.11304E-01 0.22350E+01
|
|
||||||
5.200 0.31334E+01 0.17664E+01 0.10540E-01 0.22136E+01
|
|
||||||
5.250 0.28382E+01 0.15911E+01 0.98225E-02 0.21927E+01
|
|
||||||
5.300 0.25704E+01 0.14333E+01 0.91489E-02 0.21721E+01
|
|
||||||
5.350 0.23277E+01 0.12914E+01 0.85173E-02 0.21520E+01
|
|
||||||
5.400 0.21077E+01 0.11636E+01 0.79258E-02 0.21321E+01
|
|
||||||
5.450 0.19084E+01 0.10486E+01 0.73723E-02 0.21127E+01
|
|
||||||
5.500 0.17278E+01 0.94514E+00 0.68547E-02 0.20936E+01
|
|
||||||
5.550 0.15642E+01 0.85202E+00 0.63712E-02 0.20748E+01
|
|
||||||
5.600 0.14161E+01 0.76821E+00 0.59199E-02 0.20563E+01
|
|
||||||
5.650 0.12820E+01 0.69277E+00 0.54988E-02 0.20382E+01
|
|
||||||
5.700 0.11605E+01 0.62487E+00 0.51061E-02 0.20204E+01
|
|
||||||
5.750 0.10506E+01 0.56374E+00 0.47403E-02 0.20029E+01
|
|
||||||
5.800 0.95110E+00 0.50870E+00 0.43996E-02 0.19856E+01
|
|
||||||
5.850 0.86104E+00 0.45914E+00 0.40825E-02 0.19687E+01
|
|
||||||
5.900 0.77954E+00 0.41451E+00 0.37874E-02 0.19521E+01
|
|
||||||
5.950 0.70578E+00 0.37430E+00 0.35130E-02 0.19357E+01
|
|
||||||
6.000 0.63903E+00 0.33807E+00 0.32579E-02 0.19196E+01
|
|
||||||
6.050 0.57863E+00 0.30543E+00 0.30208E-02 0.19038E+01
|
|
||||||
6.100 0.52398E+00 0.27601E+00 0.28006E-02 0.18882E+01
|
|
||||||
6.150 0.47452E+00 0.24948E+00 0.25960E-02 0.18729E+01
|
|
||||||
6.200 0.42977E+00 0.22556E+00 0.24062E-02 0.18578E+01
|
|
||||||
6.250 0.38927E+00 0.20399E+00 0.22299E-02 0.18429E+01
|
|
||||||
6.300 0.35263E+00 0.18454E+00 0.20663E-02 0.18283E+01
|
|
||||||
6.350 0.31947E+00 0.16698E+00 0.19146E-02 0.18140E+01
|
|
||||||
6.400 0.28946E+00 0.15113E+00 0.17739E-02 0.17998E+01
|
|
||||||
6.450 0.26229E+00 0.13682E+00 0.16434E-02 0.17859E+01
|
|
||||||
6.500 0.23771E+00 0.12391E+00 0.15223E-02 0.17721E+01
|
|
||||||
6.550 0.21546E+00 0.11224E+00 0.14102E-02 0.17586E+01
|
|
||||||
6.600 0.19532E+00 0.10170E+00 0.13062E-02 0.17453E+01
|
|
||||||
6.650 0.17708E+00 0.92168E-01 0.12098E-02 0.17322E+01
|
|
||||||
6.700 0.16057E+00 0.83556E-01 0.11205E-02 0.17193E+01
|
|
||||||
6.750 0.14562E+00 0.75769E-01 0.10377E-02 0.17066E+01
|
|
||||||
6.800 0.13208E+00 0.68727E-01 0.96107E-03 0.16940E+01
|
|
||||||
6.850 0.11982E+00 0.62355E-01 0.89005E-03 0.16817E+01
|
|
||||||
6.900 0.10872E+00 0.56589E-01 0.82426E-03 0.16695E+01
|
|
||||||
6.950 0.98657E-01 0.51371E-01 0.76332E-03 0.16575E+01
|
|
||||||
7.000 0.89542E-01 0.46645E-01 0.70688E-03 0.16456E+01
|
|
||||||
7.050 0.81283E-01 0.42366E-01 0.65461E-03 0.16340E+01
|
|
||||||
7.100 0.73799E-01 0.38489E-01 0.60620E-03 0.16225E+01
|
|
||||||
7.150 0.67015E-01 0.34976E-01 0.56137E-03 0.16111E+01
|
|
||||||
7.200 0.60865E-01 0.31792E-01 0.51986E-03 0.15999E+01
|
|
||||||
7.250 0.55290E-01 0.28905E-01 0.48142E-03 0.15889E+01
|
|
||||||
7.300 0.50234E-01 0.26286E-01 0.44583E-03 0.15780E+01
|
|
||||||
7.350 0.45649E-01 0.23911E-01 0.41287E-03 0.15673E+01
|
|
||||||
7.400 0.41490E-01 0.21756E-01 0.38236E-03 0.15567E+01
|
|
||||||
7.450 0.37717E-01 0.19801E-01 0.35410E-03 0.15463E+01
|
|
||||||
7.500 0.34293E-01 0.18025E-01 0.32794E-03 0.15360E+01
|
|
||||||
7.550 0.31186E-01 0.16412E-01 0.30372E-03 0.15258E+01
|
|
||||||
7.600 0.28366E-01 0.14948E-01 0.28129E-03 0.15158E+01
|
|
||||||
7.650 0.25806E-01 0.13617E-01 0.26053E-03 0.15058E+01
|
|
||||||
7.700 0.23481E-01 0.12408E-01 0.24130E-03 0.14961E+01
|
|
||||||
7.750 0.21370E-01 0.11308E-01 0.22350E-03 0.14864E+01
|
|
||||||
7.800 0.19452E-01 0.10309E-01 0.20702E-03 0.14769E+01
|
|
||||||
7.850 0.17710E-01 0.93999E-02 0.19176E-03 0.14675E+01
|
|
||||||
7.900 0.16127E-01 0.85729E-02 0.17763E-03 0.14582E+01
|
|
||||||
7.950 0.14689E-01 0.78205E-02 0.16454E-03 0.14490E+01
|
|
||||||
8.000 0.13381E-01 0.71356E-02 0.15243E-03 0.14400E+01
|
|
||||||
8.050 0.12192E-01 0.65122E-02 0.14121E-03 0.14310E+01
|
|
||||||
8.100 0.11111E-01 0.59445E-02 0.13082E-03 0.14222E+01
|
|
||||||
8.150 0.10128E-01 0.54275E-02 0.12120E-03 0.14135E+01
|
|
||||||
8.200 0.92331E-02 0.49564E-02 0.11229E-03 0.14049E+01
|
|
||||||
8.250 0.84193E-02 0.45272E-02 0.10403E-03 0.13963E+01
|
|
||||||
8.300 0.76787E-02 0.41360E-02 0.96394E-04 0.13879E+01
|
|
||||||
8.350 0.70045E-02 0.37794E-02 0.89317E-04 0.13796E+01
|
|
||||||
8.400 0.63908E-02 0.34542E-02 0.82763E-04 0.13714E+01
|
|
||||||
8.450 0.58320E-02 0.31577E-02 0.76692E-04 0.13633E+01
|
|
||||||
8.500 0.53230E-02 0.28871E-02 0.71070E-04 0.13553E+01
|
|
||||||
8.550 0.48594E-02 0.26403E-02 0.65862E-04 0.13474E+01
|
|
||||||
8.600 0.44371E-02 0.24150E-02 0.61038E-04 0.13395E+01
|
|
||||||
8.650 0.40522E-02 0.22094E-02 0.56569E-04 0.13318E+01
|
|
||||||
8.700 0.37013E-02 0.20217E-02 0.52430E-04 0.13241E+01
|
|
||||||
8.750 0.33815E-02 0.18502E-02 0.48595E-04 0.13166E+01
|
|
||||||
8.800 0.30899E-02 0.16936E-02 0.45043E-04 0.13091E+01
|
|
||||||
8.850 0.28240E-02 0.15505E-02 0.41751E-04 0.13017E+01
|
|
||||||
8.900 0.25814E-02 0.14198E-02 0.38702E-04 0.12944E+01
|
|
||||||
8.950 0.23601E-02 0.13004E-02 0.35877E-04 0.12871E+01
|
|
||||||
9.000 0.21582E-02 0.11911E-02 0.33259E-04 0.12800E+01
|
|
||||||
9.050 0.19739E-02 0.10913E-02 0.30834E-04 0.12729E+01
|
|
||||||
9.100 0.18057E-02 0.99999E-03 0.28586E-04 0.12659E+01
|
|
||||||
9.150 0.16521E-02 0.91648E-03 0.26504E-04 0.12590E+01
|
|
||||||
9.200 0.15118E-02 0.84008E-03 0.24574E-04 0.12522E+01
|
|
||||||
9.250 0.13837E-02 0.77018E-03 0.22785E-04 0.12454E+01
|
|
||||||
9.300 0.12667E-02 0.70621E-03 0.21127E-04 0.12387E+01
|
|
||||||
9.350 0.11597E-02 0.64766E-03 0.19591E-04 0.12321E+01
|
|
||||||
9.400 0.10620E-02 0.59406E-03 0.18167E-04 0.12255E+01
|
|
||||||
9.450 0.97267E-03 0.54498E-03 0.16847E-04 0.12190E+01
|
|
||||||
9.500 0.89101E-03 0.50003E-03 0.15624E-04 0.12126E+01
|
|
||||||
9.550 0.81635E-03 0.45886E-03 0.14490E-04 0.12063E+01
|
|
||||||
9.600 0.74806E-03 0.42114E-03 0.13439E-04 0.12000E+01
|
|
||||||
9.650 0.68561E-03 0.38659E-03 0.12465E-04 0.11938E+01
|
|
||||||
9.700 0.62846E-03 0.35492E-03 0.11561E-04 0.11876E+01
|
|
||||||
9.750 0.57618E-03 0.32589E-03 0.10724E-04 0.11815E+01
|
|
||||||
9.800 0.52833E-03 0.29928E-03 0.99473E-05 0.11755E+01
|
|
||||||
9.850 0.48453E-03 0.27488E-03 0.92274E-05 0.11695E+01
|
|
||||||
9.900 0.44443E-03 0.25251E-03 0.85599E-05 0.11636E+01
|
|
||||||
9.950 0.40772E-03 0.23199E-03 0.79410E-05 0.11578E+01
|
|
||||||
10.000 0.37409E-03 0.21317E-03 0.73671E-05 0.11520E+01
|
|
||||||
10.050 0.34330E-03 0.19590E-03 0.68349E-05 0.11463E+01
|
|
||||||
10.100 0.31508E-03 0.18006E-03 0.63414E-05 0.11406E+01
|
|
||||||
10.150 0.28923E-03 0.16551E-03 0.58837E-05 0.11350E+01
|
|
||||||
10.200 0.26554E-03 0.15217E-03 0.54592E-05 0.11294E+01
|
|
||||||
10.250 0.24382E-03 0.13992E-03 0.50655E-05 0.11239E+01
|
|
||||||
10.300 0.22392E-03 0.12867E-03 0.47004E-05 0.11184E+01
|
|
||||||
10.350 0.20567E-03 0.11834E-03 0.43618E-05 0.11130E+01
|
|
||||||
10.400 0.18893E-03 0.10885E-03 0.40476E-05 0.11077E+01
|
|
||||||
10.450 0.17358E-03 0.10014E-03 0.37563E-05 0.11024E+01
|
|
||||||
10.500 0.15950E-03 0.92134E-04 0.34860E-05 0.10971E+01
|
|
||||||
10.550 0.14658E-03 0.84780E-04 0.32353E-05 0.10919E+01
|
|
||||||
10.600 0.13473E-03 0.78022E-04 0.30027E-05 0.10868E+01
|
|
||||||
10.650 0.12385E-03 0.71812E-04 0.27869E-05 0.10817E+01
|
|
||||||
10.700 0.11387E-03 0.66103E-04 0.25867E-05 0.10766E+01
|
|
||||||
10.750 0.10470E-03 0.60856E-04 0.24009E-05 0.10716E+01
|
|
||||||
10.800 0.96287E-04 0.56031E-04 0.22286E-05 0.10667E+01
|
|
||||||
10.850 0.88560E-04 0.51595E-04 0.20687E-05 0.10617E+01
|
|
||||||
10.900 0.81463E-04 0.47516E-04 0.19204E-05 0.10569E+01
|
|
||||||
10.950 0.74945E-04 0.43764E-04 0.17827E-05 0.10520E+01
|
|
||||||
11.000 0.68957E-04 0.40313E-04 0.16549E-05 0.10473E+01
|
|
||||||
11.050 0.63456E-04 0.37138E-04 0.15364E-05 0.10425E+01
|
|
||||||
11.100 0.58400E-04 0.34217E-04 0.14264E-05 0.10378E+01
|
|
||||||
11.150 0.53754E-04 0.31528E-04 0.13243E-05 0.10332E+01
|
|
||||||
11.200 0.49484E-04 0.29055E-04 0.12295E-05 0.10286E+01
|
|
||||||
11.250 0.45558E-04 0.26778E-04 0.11416E-05 0.10240E+01
|
|
||||||
11.300 0.41949E-04 0.24682E-04 0.10600E-05 0.10195E+01
|
|
||||||
11.350 0.38630E-04 0.22752E-04 0.98421E-06 0.10150E+01
|
|
||||||
11.400 0.35578E-04 0.20976E-04 0.91390E-06 0.10105E+01
|
|
||||||
11.450 0.32771E-04 0.19340E-04 0.84863E-06 0.10061E+01
|
|
||||||
11.500 0.30188E-04 0.17833E-04 0.78804E-06 0.10017E+01
|
|
||||||
11.550 0.27813E-04 0.16446E-04 0.73181E-06 0.99739E+00
|
|
||||||
11.600 0.25627E-04 0.15168E-04 0.67960E-06 0.99309E+00
|
|
||||||
11.650 0.23615E-04 0.13990E-04 0.63114E-06 0.98883E+00
|
|
||||||
11.700 0.21764E-04 0.12906E-04 0.58614E-06 0.98461E+00
|
|
||||||
11.750 0.20060E-04 0.11906E-04 0.54437E-06 0.98042E+00
|
|
||||||
11.800 0.18491E-04 0.10985E-04 0.50559E-06 0.97626E+00
|
|
||||||
11.850 0.17047E-04 0.10136E-04 0.46959E-06 0.97214E+00
|
|
||||||
11.900 0.15718E-04 0.93535E-05 0.43616E-06 0.96806E+00
|
|
||||||
11.950 0.14493E-04 0.86322E-05 0.40512E-06 0.96401E+00
|
|
||||||
12.000 0.13365E-04 0.79673E-05 0.37630E-06 0.95999E+00
|
|
||||||
12.050 0.12326E-04 0.73542E-05 0.34954E-06 0.95601E+00
|
|
||||||
12.100 0.11370E-04 0.67890E-05 0.32469E-06 0.95206E+00
|
|
||||||
12.150 0.10488E-04 0.62677E-05 0.30161E-06 0.94814E+00
|
|
||||||
12.200 0.96757E-05 0.57869E-05 0.28018E-06 0.94425E+00
|
|
||||||
12.250 0.89272E-05 0.53435E-05 0.26028E-06 0.94040E+00
|
|
||||||
12.300 0.82374E-05 0.49344E-05 0.24180E-06 0.93658E+00
|
|
||||||
12.350 0.76016E-05 0.45571E-05 0.22464E-06 0.93278E+00
|
|
||||||
12.400 0.70156E-05 0.42089E-05 0.20870E-06 0.92902E+00
|
|
||||||
12.450 0.64753E-05 0.38877E-05 0.19390E-06 0.92529E+00
|
|
||||||
12.500 0.59771E-05 0.35913E-05 0.18015E-06 0.92159E+00
|
|
||||||
12.550 0.55178E-05 0.33178E-05 0.16738E-06 0.91792E+00
|
|
||||||
12.600 0.50943E-05 0.30653E-05 0.15552E-06 0.91428E+00
|
|
||||||
12.650 0.47037E-05 0.28323E-05 0.14450E-06 0.91066E+00
|
|
||||||
12.700 0.43434E-05 0.26171E-05 0.13427E-06 0.90708E+00
|
|
||||||
12.750 0.40110E-05 0.24186E-05 0.12476E-06 0.90352E+00
|
|
||||||
12.800 0.37044E-05 0.22352E-05 0.11593E-06 0.89999E+00
|
|
||||||
12.850 0.34215E-05 0.20659E-05 0.10773E-06 0.89649E+00
|
|
||||||
12.900 0.31605E-05 0.19096E-05 0.10011E-06 0.89301E+00
|
|
||||||
12.950 0.29197E-05 0.17652E-05 0.93026E-07 0.88957E+00
|
|
||||||
13.000 0.26974E-05 0.16319E-05 0.86449E-07 0.88615E+00
|
|
||||||
13.050 0.24923E-05 0.15087E-05 0.80339E-07 0.88275E+00
|
|
||||||
13.100 0.23029E-05 0.13950E-05 0.74663E-07 0.87938E+00
|
|
||||||
13.150 0.21281E-05 0.12899E-05 0.69389E-07 0.87604E+00
|
|
||||||
13.200 0.19667E-05 0.11928E-05 0.64489E-07 0.87272E+00
|
|
||||||
13.250 0.18177E-05 0.11031E-05 0.59936E-07 0.86943E+00
|
|
||||||
13.300 0.16801E-05 0.10202E-05 0.55706E-07 0.86616E+00
|
|
||||||
13.350 0.15531E-05 0.94364E-06 0.51776E-07 0.86291E+00
|
|
||||||
13.400 0.14358E-05 0.87286E-06 0.48124E-07 0.85969E+00
|
|
||||||
13.450 0.13274E-05 0.80744E-06 0.44730E-07 0.85650E+00
|
|
||||||
13.500 0.12273E-05 0.74698E-06 0.41577E-07 0.85333E+00
|
|
||||||
13.550 0.11348E-05 0.69109E-06 0.38646E-07 0.85018E+00
|
|
||||||
13.600 0.10494E-05 0.63942E-06 0.35923E-07 0.84705E+00
|
|
||||||
13.650 0.97049E-06 0.59165E-06 0.33393E-07 0.84395E+00
|
|
||||||
13.700 0.89757E-06 0.54749E-06 0.31041E-07 0.84087E+00
|
|
||||||
13.750 0.83018E-06 0.50666E-06 0.28856E-07 0.83781E+00
|
|
||||||
13.800 0.76791E-06 0.46890E-06 0.26825E-07 0.83477E+00
|
|
||||||
13.850 0.71036E-06 0.43399E-06 0.24937E-07 0.83176E+00
|
|
||||||
13.900 0.65717E-06 0.40170E-06 0.23183E-07 0.82877E+00
|
|
||||||
13.950 0.60800E-06 0.37183E-06 0.21552E-07 0.82580E+00
|
|
||||||
14.000 0.56255E-06 0.34421E-06 0.20037E-07 0.82285E+00
|
|
||||||
14.050 0.52053E-06 0.31866E-06 0.18629E-07 0.81992E+00
|
|
||||||
14.100 0.48169E-06 0.29502E-06 0.17319E-07 0.81701E+00
|
|
||||||
14.150 0.44577E-06 0.27315E-06 0.16103E-07 0.81413E+00
|
|
||||||
14.200 0.41255E-06 0.25292E-06 0.14971E-07 0.81126E+00
|
|
||||||
14.250 0.38184E-06 0.23420E-06 0.13920E-07 0.80841E+00
|
|
||||||
14.300 0.35343E-06 0.21688E-06 0.12943E-07 0.80559E+00
|
|
||||||
14.350 0.32716E-06 0.20085E-06 0.12034E-07 0.80278E+00
|
|
||||||
14.400 0.30286E-06 0.18602E-06 0.11190E-07 0.79999E+00
|
|
||||||
14.450 0.28039E-06 0.17229E-06 0.10405E-07 0.79722E+00
|
|
||||||
14.500 0.25959E-06 0.15958E-06 0.96751E-08 0.79448E+00
|
|
||||||
14.550 0.24036E-06 0.14782E-06 0.89967E-08 0.79175E+00
|
|
||||||
14.600 0.22256E-06 0.13694E-06 0.83660E-08 0.78903E+00
|
|
||||||
14.650 0.20609E-06 0.12686E-06 0.77797E-08 0.78634E+00
|
|
||||||
14.700 0.19086E-06 0.11753E-06 0.72346E-08 0.78367E+00
|
|
||||||
14.750 0.17676E-06 0.10889E-06 0.67278E-08 0.78101E+00
|
|
||||||
14.800 0.16371E-06 0.10090E-06 0.62566E-08 0.77837E+00
|
|
||||||
14.850 0.15163E-06 0.93491E-07 0.58185E-08 0.77575E+00
|
|
||||||
14.900 0.14045E-06 0.86634E-07 0.54112E-08 0.77315E+00
|
|
||||||
14.950 0.13011E-06 0.80285E-07 0.50325E-08 0.77056E+00
|
|
||||||
15.000 0.12053E-06 0.74404E-07 0.46804E-08 0.76799E+00
|
|
||||||
15.050 0.11166E-06 0.68958E-07 0.43529E-08 0.76544E+00
|
|
||||||
15.100 0.10345E-06 0.63914E-07 0.40485E-08 0.76291E+00
|
|
||||||
15.150 0.95854E-07 0.59242E-07 0.37654E-08 0.76039E+00
|
|
||||||
15.200 0.88817E-07 0.54914E-07 0.35022E-08 0.75789E+00
|
|
||||||
15.250 0.82302E-07 0.50904E-07 0.32574E-08 0.75540E+00
|
|
||||||
15.300 0.76269E-07 0.47190E-07 0.30298E-08 0.75293E+00
|
|
||||||
15.350 0.70681E-07 0.43749E-07 0.28182E-08 0.75048E+00
|
|
||||||
15.400 0.65507E-07 0.40561E-07 0.26214E-08 0.74804E+00
|
|
||||||
15.450 0.60714E-07 0.37607E-07 0.24383E-08 0.74562E+00
|
|
||||||
15.500 0.56275E-07 0.34870E-07 0.22681E-08 0.74322E+00
|
|
||||||
15.550 0.52164E-07 0.32334E-07 0.21098E-08 0.74083E+00
|
|
||||||
15.600 0.48355E-07 0.29984E-07 0.19626E-08 0.73845E+00
|
|
||||||
15.650 0.44826E-07 0.27805E-07 0.18257E-08 0.73610E+00
|
|
||||||
15.700 0.41557E-07 0.25786E-07 0.16983E-08 0.73375E+00
|
|
||||||
15.750 0.38529E-07 0.23915E-07 0.15799E-08 0.73142E+00
|
|
||||||
15.800 0.35723E-07 0.22181E-07 0.14698E-08 0.72911E+00
|
|
||||||
15.850 0.33123E-07 0.20573E-07 0.13673E-08 0.72681E+00
|
|
||||||
15.900 0.30713E-07 0.19083E-07 0.12720E-08 0.72452E+00
|
|
||||||
15.950 0.28481E-07 0.17702E-07 0.11834E-08 0.72225E+00
|
|
||||||
16.000 0.26412E-07 0.16421E-07 0.11010E-08 0.71999E+00
|
|
||||||
16.050 0.24494E-07 0.15234E-07 0.10243E-08 0.71775E+00
|
|
||||||
16.100 0.22717E-07 0.14133E-07 0.95299E-09 0.71552E+00
|
|
||||||
16.150 0.21069E-07 0.13112E-07 0.88665E-09 0.71331E+00
|
|
||||||
16.200 0.19542E-07 0.12166E-07 0.82495E-09 0.71110E+00
|
|
||||||
16.250 0.18127E-07 0.11288E-07 0.76755E-09 0.70892E+00
|
|
||||||
16.300 0.16815E-07 0.10474E-07 0.71416E-09 0.70674E+00
|
|
||||||
16.350 0.15598E-07 0.97192E-08 0.66449E-09 0.70458E+00
|
|
||||||
16.400 0.14470E-07 0.90191E-08 0.61829E-09 0.70243E+00
|
|
||||||
16.450 0.13425E-07 0.83699E-08 0.57531E-09 0.70030E+00
|
|
||||||
16.500 0.12455E-07 0.77677E-08 0.53532E-09 0.69818E+00
|
|
||||||
16.550 0.11556E-07 0.72092E-08 0.49812E-09 0.69607E+00
|
|
||||||
16.600 0.10723E-07 0.66911E-08 0.46352E-09 0.69397E+00
|
|
||||||
16.650 0.99497E-08 0.62105E-08 0.43133E-09 0.69189E+00
|
|
||||||
16.700 0.92327E-08 0.57646E-08 0.40138E-09 0.68981E+00
|
|
||||||
16.750 0.85679E-08 0.53510E-08 0.37351E-09 0.68775E+00
|
|
||||||
16.800 0.79512E-08 0.49673E-08 0.34759E-09 0.68571E+00
|
|
||||||
16.850 0.73793E-08 0.46112E-08 0.32347E-09 0.68367E+00
|
|
||||||
16.900 0.68487E-08 0.42809E-08 0.30103E-09 0.68165E+00
|
|
||||||
16.950 0.63566E-08 0.39744E-08 0.28015E-09 0.67964E+00
|
|
||||||
17.000 0.59001E-08 0.36900E-08 0.26072E-09 0.67764E+00
|
|
||||||
17.050 0.54767E-08 0.34261E-08 0.24265E-09 0.67565E+00
|
|
||||||
17.100 0.50838E-08 0.31811E-08 0.22583E-09 0.67368E+00
|
|
||||||
17.150 0.47193E-08 0.29538E-08 0.21018E-09 0.67171E+00
|
|
||||||
17.200 0.43811E-08 0.27429E-08 0.19562E-09 0.66976E+00
|
|
||||||
17.250 0.40674E-08 0.25471E-08 0.18207E-09 0.66782E+00
|
|
||||||
17.300 0.37762E-08 0.23654E-08 0.16946E-09 0.66589E+00
|
|
||||||
17.350 0.35061E-08 0.21968E-08 0.15773E-09 0.66397E+00
|
|
||||||
17.400 0.32554E-08 0.20402E-08 0.14682E-09 0.66206E+00
|
|
||||||
17.450 0.30228E-08 0.18949E-08 0.13666E-09 0.66017E+00
|
|
||||||
17.500 0.28069E-08 0.17600E-08 0.12720E-09 0.65828E+00
|
|
||||||
17.550 0.26065E-08 0.16348E-08 0.11841E-09 0.65640E+00
|
|
||||||
17.600 0.24205E-08 0.15185E-08 0.11022E-09 0.65454E+00
|
|
||||||
17.650 0.22479E-08 0.14106E-08 0.10260E-09 0.65269E+00
|
|
||||||
17.700 0.20877E-08 0.13104E-08 0.95514E-10 0.65084E+00
|
|
||||||
17.750 0.19390E-08 0.12174E-08 0.88916E-10 0.64901E+00
|
|
||||||
17.800 0.18009E-08 0.11310E-08 0.82775E-10 0.64718E+00
|
|
||||||
17.850 0.16728E-08 0.10508E-08 0.77060E-10 0.64537E+00
|
|
||||||
17.900 0.15538E-08 0.97628E-09 0.71742E-10 0.64357E+00
|
|
||||||
17.950 0.14434E-08 0.90711E-09 0.66792E-10 0.64178E+00
|
|
||||||
18.000 0.13408E-08 0.84288E-09 0.62185E-10 0.63999E+00
|
|
||||||
18.050 0.12456E-08 0.78323E-09 0.57897E-10 0.63822E+00
|
|
||||||
18.100 0.11572E-08 0.72783E-09 0.53906E-10 0.63646E+00
|
|
||||||
18.150 0.10751E-08 0.67638E-09 0.50191E-10 0.63470E+00
|
|
||||||
18.200 0.99888E-09 0.62860E-09 0.46734E-10 0.63296E+00
|
|
||||||
18.250 0.92809E-09 0.58421E-09 0.43516E-10 0.63123E+00
|
|
||||||
18.300 0.86235E-09 0.54298E-09 0.40520E-10 0.62950E+00
|
|
||||||
18.350 0.80129E-09 0.50468E-09 0.37732E-10 0.62779E+00
|
|
||||||
18.400 0.74459E-09 0.46911E-09 0.35136E-10 0.62608E+00
|
|
||||||
18.450 0.69191E-09 0.43606E-09 0.32720E-10 0.62438E+00
|
|
||||||
18.500 0.64298E-09 0.40536E-09 0.30471E-10 0.62270E+00
|
|
||||||
18.550 0.59753E-09 0.37683E-09 0.28377E-10 0.62102E+00
|
|
||||||
18.600 0.55530E-09 0.35032E-09 0.26428E-10 0.61935E+00
|
|
||||||
18.650 0.51606E-09 0.32570E-09 0.24613E-10 0.61769E+00
|
|
||||||
18.700 0.47960E-09 0.30281E-09 0.22923E-10 0.61604E+00
|
|
||||||
18.750 0.44571E-09 0.28154E-09 0.21349E-10 0.61439E+00
|
|
||||||
18.800 0.41422E-09 0.26177E-09 0.19884E-10 0.61276E+00
|
|
||||||
18.850 0.38493E-09 0.24340E-09 0.18519E-10 0.61113E+00
|
|
||||||
18.900 0.35771E-09 0.22631E-09 0.17248E-10 0.60952E+00
|
|
||||||
18.950 0.33239E-09 0.21043E-09 0.16063E-10 0.60791E+00
|
|
||||||
19.000 0.30883E-09 0.19565E-09 0.14960E-10 0.60631E+00
|
|
||||||
19.050 0.28692E-09 0.18190E-09 0.13931E-10 0.60472E+00
|
|
||||||
19.100 0.26651E-09 0.16911E-09 0.12972E-10 0.60314E+00
|
|
||||||
19.150 0.24751E-09 0.15719E-09 0.12076E-10 0.60156E+00
|
|
||||||
19.200 0.22980E-09 0.14609E-09 0.11241E-10 0.59999E+00
|
|
||||||
19.250 0.21329E-09 0.13574E-09 0.10460E-10 0.59844E+00
|
|
||||||
19.300 0.19788E-09 0.12607E-09 0.97296E-11 0.59689E+00
|
|
||||||
19.350 0.18348E-09 0.11704E-09 0.90455E-11 0.59534E+00
|
|
||||||
19.400 0.17002E-09 0.10858E-09 0.84036E-11 0.59381E+00
|
|
||||||
19.450 0.15742E-09 0.10064E-09 0.78000E-11 0.59228E+00
|
|
||||||
19.500 0.14560E-09 0.93171E-10 0.72308E-11 0.59076E+00
|
|
||||||
19.550 0.13448E-09 0.86115E-10 0.66919E-11 0.58925E+00
|
|
||||||
19.600 0.12402E-09 0.79421E-10 0.61796E-11 0.58775E+00
|
|
||||||
19.650 0.11413E-09 0.73033E-10 0.56896E-11 0.58625E+00
|
|
||||||
19.700 0.10475E-09 0.66895E-10 0.52178E-11 0.58477E+00
|
|
||||||
19.750 0.95834E-10 0.60946E-10 0.47595E-11 0.58329E+00
|
|
||||||
19.800 0.87307E-10 0.55120E-10 0.43098E-11 0.58181E+00
|
|
||||||
19.850 0.79114E-10 0.49344E-10 0.38632E-11 0.58035E+00
|
|
||||||
19.900 0.71193E-10 0.43537E-10 0.34137E-11 0.57889E+00
|
|
||||||
19.950 0.63484E-10 0.37606E-10 0.29542E-11 0.57744E+00
|
|
||||||
20.000 0.55923E-10 0.31447E-10 0.24767E-11 0.57599E+00
|
|
||||||
|
|
|
@ -1,7 +0,0 @@
|
||||||
proton,kappa,filling,gnodes,fnodes,Ebind
|
|
||||||
1,-1,1,0,0,36.78504
|
|
||||||
1,-2,1,0,0,16.42374
|
|
||||||
1,1,1,0,1,8.61761
|
|
||||||
0,-1,1,0,0,40.99365
|
|
||||||
0,-2,1,0,0,20.30183
|
|
||||||
0,1,1,0,1,12.37064
|
|
||||||
|
|
|
@ -1,43 +0,0 @@
|
||||||
using DelimitedFiles, Plots
|
|
||||||
include("../nucleons.jl")
|
|
||||||
|
|
||||||
# test data generated from Hartree.f
|
|
||||||
# format: x S(x) V(x) R(x) A(x)
|
|
||||||
test_data = readdlm("test/LinearFldsFSUGarnet.csv")
|
|
||||||
xs = test_data[:, 1]
|
|
||||||
Ss = test_data[:, 2]
|
|
||||||
Vs = test_data[:, 3]
|
|
||||||
Rs = test_data[:, 4]
|
|
||||||
As = test_data[:, 5]
|
|
||||||
|
|
||||||
p = false
|
|
||||||
r_max = maximum(xs)
|
|
||||||
divs = length(xs) - 1
|
|
||||||
s = system(8, 8, parameters(), r_max, divs)
|
|
||||||
|
|
||||||
s.Φ0 = Ss
|
|
||||||
s.W0 = Vs
|
|
||||||
s.B0 = Rs
|
|
||||||
s.A0 = As
|
|
||||||
|
|
||||||
κs, Es = findAllOrbitals(p, s)
|
|
||||||
Ebinds = (p ? M_p : M_n) .- Es
|
|
||||||
spec = 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_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])
|
|
||||||
|
|
||||||
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])
|
|
||||||
|
|
||||||
xlabel!("κ")
|
|
||||||
ylabel!("Binding energy (MeV)")
|
|
||||||
yflip!()
|
|
||||||
|
|
@ -1,14 +0,0 @@
|
||||||
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)
|
|
||||||
|
|
||||||
solve_system!(s; live_plots=false)
|
|
||||||
|
|
@ -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,39 +0,0 @@
|
||||||
proton,kappa,filling,gnodes,fnodes,Ebind
|
|
||||||
0,-1,1,0,0,63.80958
|
|
||||||
0,-2,1,0,0,57.09607
|
|
||||||
0,1,1,0,1,56.51421
|
|
||||||
0,-3,1,0,0,48.96921
|
|
||||||
0,2,1,0,1,47.53564
|
|
||||||
0,-1,1,1,1,44.05962
|
|
||||||
0,-4,1,0,0,39.82046
|
|
||||||
0,3,1,0,1,37.11666
|
|
||||||
0,-2,1,1,1,32.52852
|
|
||||||
0,1,1,1,2,31.35759
|
|
||||||
0,-5,1,0,0,30.03271
|
|
||||||
0,4,1,0,1,25.74641
|
|
||||||
0,-3,1,1,1,21.54695
|
|
||||||
0,-6,1,0,0,19.94124
|
|
||||||
0,2,1,1,2,19.70747
|
|
||||||
0,-1,1,2,2,19.05405
|
|
||||||
0,5,1,0,1,13.96937
|
|
||||||
0,-4,1,1,1,11.33989
|
|
||||||
0,-7,1,0,0,9.84013
|
|
||||||
0,3,1,1,2,9.04831
|
|
||||||
0,-2,1,2,2,8.57808
|
|
||||||
0,1,1,2,3,7.67567
|
|
||||||
1,-1,1,0,0,50.48167
|
|
||||||
1,-2,1,0,0,44.63156
|
|
||||||
1,1,1,0,1,43.88757
|
|
||||||
1,-3,1,0,0,37.28703
|
|
||||||
1,2,1,0,1,35.553
|
|
||||||
1,-1,1,1,1,31.53505
|
|
||||||
1,-4,1,0,0,28.92962
|
|
||||||
1,3,1,0,1,25.80862
|
|
||||||
1,-2,1,1,1,21.04461
|
|
||||||
1,-5,1,0,0,19.94885
|
|
||||||
1,1,1,1,2,19.8522
|
|
||||||
1,4,1,0,1,15.19352
|
|
||||||
1,-3,1,1,1,11.01997
|
|
||||||
1,-6,1,0,0,10.64656
|
|
||||||
1,2,1,1,2,9.22511
|
|
||||||
1,-1,1,2,2,8.31272
|
|
||||||
|
|
|
@ -1,9 +1,9 @@
|
||||||
using DelimitedFiles, Plots
|
using DelimitedFiles, Interpolations, Plots
|
||||||
include("../mesons.jl")
|
include("../mesons.jl")
|
||||||
|
|
||||||
# test data generated from Hartree.f
|
# test data generated from Hartree.f
|
||||||
# format: x S(x) V(x) R(x) A(x)
|
# 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]
|
xs_bench = bench_data[:, 1]
|
||||||
Φ0_bench = bench_data[:, 2]
|
Φ0_bench = bench_data[:, 2]
|
||||||
W0_bench = bench_data[:, 3]
|
W0_bench = bench_data[:, 3]
|
||||||
|
|
@ -16,20 +16,15 @@ plot(xs_bench, hcat(Φ0_bench, W0_bench, B0_bench, A0_bench), layout=4, label=["
|
||||||
# format: x Rhos(n) Rhov(n) Rhot(n) Rhos(p) Rhov(p) Rhot(p)
|
# format: x Rhos(n) Rhov(n) Rhot(n) Rhos(p) Rhov(p) Rhot(p)
|
||||||
test_data = readdlm("test/Pb208DensFSUGarnet.csv")
|
test_data = readdlm("test/Pb208DensFSUGarnet.csv")
|
||||||
xs = test_data[:, 1]
|
xs = test_data[:, 1]
|
||||||
|
ρ_sn = test_data[:, 2]
|
||||||
|
ρ_vn = test_data[:, 3]
|
||||||
|
ρ_sp = test_data[:, 5]
|
||||||
|
ρ_vp = test_data[:, 6]
|
||||||
|
|
||||||
N_p = 82
|
|
||||||
N_n = 126
|
|
||||||
param = parameters(read("test/Pb208ParamsFSUGarnet.txt", String))
|
|
||||||
r_max = maximum(xs)
|
r_max = maximum(xs)
|
||||||
divs = length(xs) - 1
|
divs = length(xs) - 1
|
||||||
s = system(N_p, N_n, param, r_max, divs)
|
|
||||||
|
|
||||||
s.ρ_sn = test_data[:, 2]
|
(Φ0, W0, B0, A0) = solveMesonWfs(ρ_sp, ρ_vp, ρ_sn, ρ_vn, r_max, divs, 1; initial_sol=(Φ0_bench, W0_bench, B0_bench, A0_bench))
|
||||||
s.ρ_vn = test_data[:, 3]
|
|
||||||
s.ρ_sp = test_data[:, 5]
|
|
||||||
s.ρ_vp = test_data[:, 6]
|
|
||||||
|
|
||||||
(Φ0, W0, B0, A0) = solveMesonFields!(s, 200) # many iterations needed without an initial solution
|
|
||||||
|
|
||||||
plot!(xs, hcat(Φ0, W0, B0, A0), layout=4, label=["Φ0" "W0" "B0" "A0"])
|
plot!(xs, hcat(Φ0, W0, B0, A0), layout=4, label=["Φ0" "W0" "B0" "A0"])
|
||||||
xlabel!("r (fm)")
|
xlabel!("r (fm)")
|
||||||
|
|
@ -0,0 +1,32 @@
|
||||||
|
using DelimitedFiles, Interpolations, Plots
|
||||||
|
include("../nucleons.jl")
|
||||||
|
|
||||||
|
# test data generated from Hartree.f
|
||||||
|
# format: x S(x) V(x) R(x) A(x)
|
||||||
|
test_data = readdlm("test/Pb208FldsFSUGarnet.csv")
|
||||||
|
xs = test_data[:, 1]
|
||||||
|
Ss = test_data[:, 2]
|
||||||
|
Vs = test_data[:, 3]
|
||||||
|
Rs = test_data[:, 4]
|
||||||
|
As = test_data[:, 5]
|
||||||
|
|
||||||
|
S_interp = linear_interpolation(xs, Ss)
|
||||||
|
V_interp = linear_interpolation(xs, Vs)
|
||||||
|
R_interp = linear_interpolation(xs, Rs)
|
||||||
|
A_interp = linear_interpolation(xs, As)
|
||||||
|
|
||||||
|
κ = -1
|
||||||
|
p = true
|
||||||
|
r_max = maximum(xs)
|
||||||
|
E_min = 880
|
||||||
|
E_max = 939
|
||||||
|
|
||||||
|
boundEs = findEs(κ, p, S_interp, V_interp, R_interp, A_interp, r_max, E_min, E_max)
|
||||||
|
println("bound E = $boundEs")
|
||||||
|
|
||||||
|
Es = collect(E_min:0.5:E_max)
|
||||||
|
boundaryVals = [boundaryValue(κ, p, E, S_interp, V_interp, R_interp, A_interp, r_max)^2 for E in Es]
|
||||||
|
|
||||||
|
plot(Es, boundaryVals, yscale=:log10, label="g(r_max)^2")
|
||||||
|
vline!(boundEs, label="bound E")
|
||||||
|
xlabel!("E (MeV)")
|
||||||
|
|
@ -1,35 +1,38 @@
|
||||||
using DelimitedFiles, Plots
|
using DelimitedFiles, Interpolations, Plots
|
||||||
include("../nucleons.jl")
|
include("../nucleons.jl")
|
||||||
|
|
||||||
# test data generated from Hartree.f
|
# test data generated from Hartree.f
|
||||||
# format: x S(x) V(x) R(x) A(x)
|
# 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]
|
xs = test_data[:, 1]
|
||||||
Ss = test_data[:, 2]
|
Ss = test_data[:, 2]
|
||||||
Vs = test_data[:, 3]
|
Vs = test_data[:, 3]
|
||||||
Rs = test_data[:, 4]
|
Rs = test_data[:, 4]
|
||||||
As = test_data[:, 5]
|
As = test_data[:, 5]
|
||||||
|
|
||||||
|
S_interp = linear_interpolation(xs, Ss)
|
||||||
|
V_interp = linear_interpolation(xs, Vs)
|
||||||
|
R_interp = linear_interpolation(xs, Rs)
|
||||||
|
A_interp = linear_interpolation(xs, As)
|
||||||
|
|
||||||
N_p = 82
|
N_p = 82
|
||||||
N_n = 126
|
N_n = 126
|
||||||
r_max = maximum(xs)
|
r_max = maximum(xs)
|
||||||
divs = length(xs) - 1
|
E_min = 880
|
||||||
s = system(N_p, N_n, parameters(), r_max, divs)
|
E_max = 939
|
||||||
|
divs = 400
|
||||||
|
|
||||||
s.Φ0 = Ss
|
rs = range(0, r_max, length=divs+1)
|
||||||
s.W0 = Vs
|
|
||||||
s.B0 = Rs
|
|
||||||
s.A0 = As
|
|
||||||
|
|
||||||
solveNucleonDensity!(true, s)
|
(ρ_sp, ρ_vp) = calculateNucleonDensity(N_p, true, S_interp, V_interp, R_interp, A_interp, r_max, divs, E_min, E_max)
|
||||||
|
|
||||||
p_sp = plot(rs(s), s.ρ_sp, xlabel="r (fm)", label="ρₛₚ(r) calculated")
|
p_sp = plot(rs, ρ_sp, xlabel="r (fm)", label="ρₛₚ(r) calculated")
|
||||||
p_vp = plot(rs(s), s.ρ_vp, xlabel="r (fm)", label="ρᵥₚ(r) calculated")
|
p_vp = plot(rs, ρ_vp, xlabel="r (fm)", label="ρᵥₚ(r) calculated")
|
||||||
|
|
||||||
solveNucleonDensity!(false, s)
|
(ρ_sn, ρ_vn) = calculateNucleonDensity(N_n, false, S_interp, V_interp, R_interp, A_interp, r_max, divs, E_min, E_max)
|
||||||
|
|
||||||
p_sn = plot(rs(s), s.ρ_sn, xlabel="r (fm)", label="ρₛₙ(r) calculated")
|
p_sn = plot(rs, ρ_sn, xlabel="r (fm)", label="ρₛₙ(r) calculated")
|
||||||
p_vn = plot(rs(s), s.ρ_vn, xlabel="r (fm)", label="ρᵥₙ(r) calculated")
|
p_vn = plot(rs, ρ_vn, xlabel="r (fm)", label="ρᵥₙ(r) calculated")
|
||||||
|
|
||||||
# benchmark data generated from Hartree.f
|
# benchmark data generated from Hartree.f
|
||||||
# format: x Rhos(n) Rhov(n) Rhot(n) Rhos(p) Rhov(p) Rhot(p)
|
# format: x Rhos(n) Rhov(n) Rhot(n) Rhos(p) Rhov(p) Rhot(p)
|
||||||
|
|
|
||||||
|
|
@ -1,4 +1,4 @@
|
||||||
using DelimitedFiles, Plots
|
using DelimitedFiles, Interpolations, Plots
|
||||||
include("../nucleons.jl")
|
include("../nucleons.jl")
|
||||||
|
|
||||||
# test data generated from Hartree.f
|
# test data generated from Hartree.f
|
||||||
|
|
@ -10,36 +10,21 @@ Vs = test_data[:, 3]
|
||||||
Rs = test_data[:, 4]
|
Rs = test_data[:, 4]
|
||||||
As = test_data[:, 5]
|
As = test_data[:, 5]
|
||||||
|
|
||||||
|
S_interp = linear_interpolation(xs, Ss)
|
||||||
|
V_interp = linear_interpolation(xs, Vs)
|
||||||
|
R_interp = linear_interpolation(xs, Rs)
|
||||||
|
A_interp = linear_interpolation(xs, As)
|
||||||
|
|
||||||
p = true
|
p = true
|
||||||
N_p = 82
|
N = 82
|
||||||
N_n = 126
|
|
||||||
r_max = maximum(xs)
|
r_max = maximum(xs)
|
||||||
divs = length(xs) - 1
|
E_min = 880
|
||||||
s = system(N_p, N_n, parameters(), r_max, divs)
|
E_max = 939
|
||||||
|
|
||||||
s.Φ0 = Ss
|
κs, Es = findAllOrbitals(p, S_interp, V_interp, R_interp, A_interp, r_max, E_min, E_max)
|
||||||
s.W0 = Vs
|
occ = fillNucleons(N, κs, Es)
|
||||||
s.B0 = Rs
|
|
||||||
s.A0 = As
|
|
||||||
|
|
||||||
κs, Es = findAllOrbitals(p, s)
|
|
||||||
Ebinds = (p ? M_p : M_n) .- Es
|
|
||||||
spec = fillNucleons(Z_or_N(s, p), κs, Es)
|
|
||||||
|
|
||||||
# format: proton, kappa, filling, gnodes, fnodes, Ebind
|
|
||||||
bench_data, _ = readdlm("test/Pb208Spectrum.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_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])
|
|
||||||
|
|
||||||
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])
|
|
||||||
|
|
||||||
|
scatter(κs, Es, label=(p ? "proton" : "neutron") * " spectrum")
|
||||||
|
annotate!(κs .+ 0.3, Es, occ)
|
||||||
xlabel!("κ")
|
xlabel!("κ")
|
||||||
ylabel!("Binding energy (MeV)")
|
ylabel!("E (MeV)")
|
||||||
yflip!()
|
|
||||||
|
|
@ -1,4 +1,4 @@
|
||||||
using DelimitedFiles, Plots
|
using DelimitedFiles, Interpolations, Plots
|
||||||
include("../nucleons.jl")
|
include("../nucleons.jl")
|
||||||
|
|
||||||
# test data generated from Hartree.f
|
# test data generated from Hartree.f
|
||||||
|
|
@ -10,24 +10,26 @@ Vs = test_data[:, 3]
|
||||||
Rs = test_data[:, 4]
|
Rs = test_data[:, 4]
|
||||||
As = test_data[:, 5]
|
As = test_data[:, 5]
|
||||||
|
|
||||||
|
S_interp = linear_interpolation(xs, Ss)
|
||||||
|
V_interp = linear_interpolation(xs, Vs)
|
||||||
|
R_interp = linear_interpolation(xs, Rs)
|
||||||
|
A_interp = linear_interpolation(xs, As)
|
||||||
|
|
||||||
κ = -1
|
κ = -1
|
||||||
p = true
|
p = true
|
||||||
r_max = maximum(xs)
|
r_max = maximum(xs)
|
||||||
divs = length(xs) - 1
|
E_min = 880
|
||||||
s = system(r_max, divs)
|
E_max = 939
|
||||||
|
|
||||||
s.Φ0 = Ss
|
groundE = findEs(κ, p, S_interp, V_interp, R_interp, A_interp, r_max, E_min, E_max) |> minimum
|
||||||
s.W0 = Vs
|
|
||||||
s.B0 = Rs
|
|
||||||
s.A0 = As
|
|
||||||
|
|
||||||
groundE = findEs(κ, p, s) |> minimum
|
|
||||||
println("ground state E = $groundE")
|
println("ground state E = $groundE")
|
||||||
|
|
||||||
wf = solveNucleonWf(κ, p, groundE, s)
|
divs = 400
|
||||||
|
wf = solveNucleonWf(κ, p, groundE, S_interp, V_interp, R_interp, A_interp, r_max, divs)
|
||||||
|
rs = range(0, r_max, length=divs+1)
|
||||||
gs = wf[1, :]
|
gs = wf[1, :]
|
||||||
fs = wf[2, :]
|
fs = wf[2, :]
|
||||||
|
|
||||||
plot(rs(s), gs, label="g(r)")
|
plot(rs, gs, label="g(r)")
|
||||||
plot!(rs(s), fs, label="f(r)")
|
plot!(rs, fs, label="f(r)")
|
||||||
xlabel!("r (fm)")
|
xlabel!("r (fm)")
|
||||||
|
|
|
||||||
|
|
@ -1,38 +0,0 @@
|
||||||
using DelimitedFiles, Plots
|
|
||||||
include("../nucleons.jl")
|
|
||||||
|
|
||||||
# test data generated from Hartree.f
|
|
||||||
# format: x S(x) V(x) R(x) A(x)
|
|
||||||
test_data = readdlm("test/Pb208FldsFSUGarnet.csv")
|
|
||||||
xs = test_data[:, 1]
|
|
||||||
Ss = test_data[:, 2]
|
|
||||||
Vs = test_data[:, 3]
|
|
||||||
Rs = test_data[:, 4]
|
|
||||||
As = test_data[:, 5]
|
|
||||||
|
|
||||||
κ = -1
|
|
||||||
p = false
|
|
||||||
r_max = maximum(xs)
|
|
||||||
divs = length(xs) - 1
|
|
||||||
s = system(r_max, divs)
|
|
||||||
|
|
||||||
s.Φ0 = Ss
|
|
||||||
s.W0 = Vs
|
|
||||||
s.B0 = Rs
|
|
||||||
s.A0 = As
|
|
||||||
|
|
||||||
E_min = 850.0
|
|
||||||
E_max = 938.0
|
|
||||||
|
|
||||||
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]
|
|
||||||
|
|
||||||
plot(Es, dets, yscale=:log10, label="determinant^2")
|
|
||||||
vline!(boundEs, label="bound E")
|
|
||||||
xlabel!("E (MeV)")
|
|
||||||
Loading…
Reference in New Issue