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

Some static typing

This commit is contained in:
Nuwan Yapa 2025-02-18 19:18:00 -05:00
parent 5fced75f41
commit b8e82644a8
1 changed files with 8 additions and 8 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

16
nucleons.jl
View File

@ -16,7 +16,7 @@ const r_reg = 1E-8 # fm # regulator for the centrifugal term
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.
Reference: P. Giuliani, K. Godbey, E. Bonilla, F. Viens, and J. Piekarewicz, Frontiers in Physics 10, (2023)"
function dirac!(du, u, (κ, E, f1, f2), r)
function dirac!(du::Vector{Float64}, u::Vector{Float64}, (κ, E, f1, f2), r::Float64) # TODO: Static typing
(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
@ -25,7 +25,7 @@ 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, s::system)
function optimized_dirac_potentials(p::Bool, s::system)
M = p ? M_p : M_n
f1s = zero_array(s)
@ -44,7 +44,7 @@ end
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, E, s::system; shooting=true, normalize=true, algo=Tsit5())
function solveNucleonWf(κ, p::Bool, E, s::system; shooting=true, normalize=true, algo=Tsit5())
(f1, f2) = optimized_dirac_potentials(p, s)
if shooting
@ -81,7 +81,7 @@ 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())
function boundaryValueFunc(κ, p::Bool, s::system; dtype=Float64, algo=Tsit5())
prob = ODEProblem(dirac!, convert.(dtype, [0, 1]), (0, s.r_max))
(f1, f2) = optimized_dirac_potentials(p, s)
func(E) = solve(prob, algo, p=(κ, E, f1, f2), saveat=[s.r_max], save_idxs=[1])[1, 1]
@ -91,7 +91,7 @@ end
"Find all bound energies between E_min (=850.0) and E_max (=938.0) 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, s::system, E_min=850.0, E_max=938.0, tol_digits=5)
function findEs(κ, p::Bool, s::system, E_min=850.0, E_max=938.0, tol_digits=5)
func = boundaryValueFunc(κ, p, s)
Es = find_all_zeros(func, E_min, E_max; partitions=20, tol=1/10^tol_digits)
return unique(E -> round(E; digits=tol_digits), Es)
@ -99,7 +99,7 @@ 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, s::system, E_min=850.0, E_max=938.0)
function findAllOrbitals(p::Bool, s::system, E_min=850.0, E_max=938.0)
κs = Int[]
Es = Float64[]
# start from κ=-1 and go both up and down
@ -140,7 +140,7 @@ 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(κs, Es, occs, p, s::system)
function calculateNucleonDensity(κs, Es, occs, p::Bool, s::system)
ρr2 = zeros(2, s.divs + 1) # ρ×r² values
for (κ, E, occ) in zip(κs, Es, occs)
@ -161,7 +161,7 @@ end
"Solve the Dirac equation and calculate scalar and vector densities of a nucleon species where
the other parameters are defined above"
function solveNucleonDensity(p, s::system, E_min=850.0, E_max=938.0)
function solveNucleonDensity(p::Bool, s::system, E_min=850.0, E_max=938.0)
κs, Es = findAllOrbitals(p, s, E_min, E_max)
occs = fillNucleons(Z_or_N(s, p), κs, Es)
return calculateNucleonDensity(κs, Es, occs, p, s)