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

Fixed unstability of wave functions

This commit is contained in:
Nuwan Yapa 2024-07-16 18:29:35 -04:00
parent 6538c30fa8
commit f77e5dfed8
2 changed files with 18 additions and 11 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

19
dirac.jl
View File

@ -26,18 +26,18 @@ end
"Solve the Dirac equation and return the wave function u(r)=[g(r), f(r)] where
divs is the number of mesh divisions if the solution should be discretized as a 2×(1+divs) matrix (keep divs=0 to obtain an interpolating function),
the other parameters are the same from dirac!(...)."
function solveWf(κ, p, E, Φ0, W0, B0, A0, r_max, divs::Int=0)
prob = ODEProblem(dirac!, [0, 1], (0, r_max))
sol = solve(prob, RK4(), p=(κ, p, E, Φ0, W0, B0, A0), saveat=(divs == 0 ? [] : r_max/divs))
function solveWf(κ, p, E, Φ0, W0, B0, A0, r_max, divs::Int=0; dtype=BigFloat, algo=Feagin12())
prob = ODEProblem(dirac!, convert.(dtype, [0, 1]), (0, r_max))
sol = solve(prob, algo, p=(κ, p, E, Φ0, W0, B0, A0), saveat=(divs == 0 ? [] : r_max/divs))
return divs == 0 ? sol : hcat(sol.u...)
end
"Solve the Dirac equation and return g(r=r_max) where
r_max is the outer boundary in fm,
the other parameters are the same from dirac!(...)."
function boundaryValue(κ, p, E, Φ0, W0, B0, A0, r_max)
prob = ODEProblem(dirac!, [0, 1], (0, r_max))
sol = solve(prob, RK4(), p=(κ, p, E, Φ0, W0, B0, A0), saveat=[r_max], save_idxs=[1])
function boundaryValue(κ, p, E, Φ0, W0, B0, A0, r_max; dtype=Float64, algo=RK4())
prob = ODEProblem(dirac!, convert.(dtype, [0, 1]), (0, r_max))
sol = solve(prob, algo, p=(κ, p, E, Φ0, W0, B0, A0), saveat=[r_max], save_idxs=[1])
return sol[1, 1]
end
@ -48,6 +48,13 @@ function findEs(κ, p, Φ0, W0, B0, A0, r_max, E_min=0, E_max=(p ? M_p : M_n))
return find_zeros(f, (E_min, E_max))
end
"Find more precise bound energies for a given list of Es where
the other parameters are the same from dirac!(...)."
function refineEs(κ, p, Φ0, W0, B0, A0, r_max, Es)
f(E) = boundaryValue(κ, p, E, Φ0, W0, B0, A0, r_max; dtype=BigFloat, algo=Feagin12())
return [find_zero(f, E) for E in 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, Φ0, W0, B0, A0, r_max, E_min=0, E_max=(p ? M_p : M_n))

10
test/Pb208_wf.jl
View File

@ -21,14 +21,14 @@ 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)
groundE = minimum(boundEs)
approxE = findEs(κ, p, S_interp, V_interp, R_interp, A_interp, r_max, E_min, E_max) |> minimum
groundE = refineEs(κ, p, S_interp, V_interp, R_interp, A_interp, r_max, [approxE])[1]
println("ground state E = $groundE")
plot_r_max = r_max * 0.75
divs = 50
wf = solveWf(κ, p, groundE, S_interp, V_interp, R_interp, A_interp, plot_r_max, divs)
rs = collect(0: plot_r_max/divs : plot_r_max)
wf = solveWf(κ, p, groundE, S_interp, V_interp, R_interp, A_interp, r_max, divs)
rs = range(0, r_max, length=divs+1)
gs = wf[1, :]
fs = wf[2, :]