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

Regular instead of r_min

This commit is contained in:
Nuwan Yapa 2024-06-21 10:58:51 -04:00
parent d23792b977
commit 6deee7e47c
2 changed files with 10 additions and 9 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

15
dirac.jl
View File

@ -4,6 +4,8 @@ using DifferentialEquations, Roots
M_n = 939.5654133 # Neutron mass in MeV/c2
M_p = 938.2720813 # Proton mass in MeV/c2
const r_reg = 1E-6 # regulator for the centrifugal term in fm
"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,
@ -17,23 +19,22 @@ function dirac!(du, u, (κ, p, E, Φ0, W0, B0, A0), r)
common1 = E - W0(r) - (p - 0.5) * B0(r) - p * A0(r)
common2 = M - Φ0(r)
(g, f) = u
du[1] = -(κ/r) * g + (common1 + common2) * f / ħc
du[2] = (κ/r) * f - (common1 - common2) * g / ħc
du[1] = -(κ/(r + r_reg)) * g + (common1 + common2) * f / ħc
du[2] = (κ/(r + r_reg)) * f - (common1 - common2) * g / ħc
end
"Solve the Dirac equation and return g(r=r_max) for given scalar and vector potentials where
r_max is the outer boundary in fm,
r_min (=r_max/1000) is inside boundary in fm which cannot be 0 due to the centrifugal term,
the other parameters are the same from dirac!(...)."
function boundaryValue(κ, p, E, Φ0, W0, B0, A0, r_max, r_min=r_max/1000)
prob = ODEProblem(dirac!, [0, 1], (r_min, r_max))
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))
return sol(r_max)[1]
end
"Find all bound energies between E_min (=0) and E_max (=mass) where
the other parameters are the same from dirac!(...)."
function findEs(κ, p, Φ0, W0, B0, A0, r_max, r_min=r_max/1000, E_min=0, E_max=(p ? M_p : M_n))
f(E) = boundaryValue(κ, p, E, Φ0, W0, B0, A0, r_max, r_min)
function findEs(κ, p, Φ0, W0, B0, A0, r_max, E_min=0, E_max=(p ? M_p : M_n))
f(E) = boundaryValue(κ, p, E, Φ0, W0, B0, A0, r_max)
return find_zeros(f, (E_min, E_max))
end

4
test/Pb208.jl
View File

@ -18,10 +18,10 @@ A_interp = linear_interpolation(xs, As)
κ = -1
p = true
r_max = maximum(xs)
E_min = 850
E_min = 880
E_max = 939
boundEs = findEs(κ, p, S_interp, V_interp, R_interp, A_interp, r_max, r_max/1000, E_min, E_max)
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)