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

Changed initial boundary condition

This commit is contained in:
Nuwan Yapa 2025-04-08 12:29:01 -04:00
parent 268ce12b6c
commit 897fdf9214
1 changed files with 5 additions and 2 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

7
nucleons.jl
View File

@ -51,6 +51,9 @@ function asymp_BC(κ::Int, p::Bool, E::Float64, r::Float64)
return [g, f] return [g, f]
end 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, shooting method divides the interval into two partitions at r_max/2, ensuring convergence at both r=0 and r=r_max,
@ -68,7 +71,7 @@ function solveNucleonWf(κ, p::Bool, E, s::system; shooting=true, normalize=true
next_r_max = s.r_max next_r_max = s.r_max
end end
prob = ODEProblem(dirac!, [0, 1], (0, next_r_max)) prob = ODEProblem(dirac!, init_BC(), (0, next_r_max))
sol = solve(prob, algo, p=(κ, E, f1, f2), saveat=Δr(s)) sol = solve(prob, algo, p=(κ, E, f1, f2), saveat=Δr(s))
wf = hcat(sol.u...) wf = hcat(sol.u...)
@ -97,7 +100,7 @@ end
where is r is in fm." where is r is in fm."
function determinantFunc(κ, p::Bool, s::system, r::Float64=s.r_max/2, algo=Vern9()) function determinantFunc(κ, p::Bool, s::system, r::Float64=s.r_max/2, algo=Vern9())
(f1, f2) = optimized_dirac_potentials(p, s) (f1, f2) = optimized_dirac_potentials(p, s)
prob_left = ODEProblem(dirac!, [0.0, 1.0], (0, r)) prob_left = ODEProblem(dirac!, init_BC(), (0, r))
function func(E) function func(E)
prob_right = ODEProblem(dirac!, asymp_BC(κ, p, E, s.r_max), (s.r_max, r)) 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_left = solve(prob_left, algo, p=(κ, E, f1, f2), saveat=[r])