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

Implemented determinant and removed boundary value

This commit is contained in:
Nuwan Yapa 2025-02-19 12:49:43 -05:00
parent b8e82644a8
commit df8ed2fe12
2 changed files with 16 additions and 11 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

19
nucleons.jl
View File

@ -78,13 +78,18 @@ function solveNucleonWf(κ, p::Bool, E, s::system; shooting=true, normalize=true
return wf
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::Bool, s::system; dtype=Float64, algo=Tsit5())
prob = ODEProblem(dirac!, convert.(dtype, [0, 1]), (0, s.r_max))
"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=Tsit5())
(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]
prob_left = ODEProblem(dirac!, [0.0, 1.0], (0, r))
prob_right = ODEProblem(dirac!, [0.0, 1.0], (s.r_max, r))
function func(E)
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
@ -92,7 +97,7 @@ end
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=5)
func = boundaryValueFunc(κ, p, s)
func = determinantFunc(κ, 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)
end

8
test/Pb208_nucleon_basic.jl → test/Pb208_shooting.jl
View File

@ -29,10 +29,10 @@ boundEs = findEs(κ, p, s, E_min, E_max)
binding_Es = (p ? M_p : M_n) .- boundEs
println("binding energies = $binding_Es")
func = boundaryValueFunc(κ, p, s)
Es = collect(E_min:0.5:E_max)
boundaryVals = [func(E)^2 for E in Es]
func = determinantFunc(κ, p, s)
Es = collect(E_min:0.1:E_max)
dets = [func(E)^2 for E in Es]
plot(Es, boundaryVals, yscale=:log10, label="g(r_max)^2")
plot(Es, dets, yscale=:log10, label="determinant^2")
vline!(boundEs, label="bound E")
xlabel!("E (MeV)")