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

Fix silly mistake

This commit is contained in:
Nuwan Yapa 2025-02-04 16:56:48 -05:00
parent 83f5039eca
commit 9f1ed7cc6e
1 changed files with 5 additions and 5 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

10
nucleons.jl
View File

@ -28,18 +28,18 @@ end
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,
the other parameters are the same from dirac!(...)." the other parameters are the same from dirac!(...)."
function solveNucleonWf(κ, p, E, Φ0, W0, B0, A0, r_max, divs; shooting=true, normalize=true) function solveNucleonWf(κ, p, E, Φ0, W0, B0, A0, r_max, divs; shooting=true, normalize=true)
saveat = r_max / divs Δr = r_max / divs
if shooting if shooting
@assert divs % 2 == 0 "divs must be an even number when shooting=true" @assert divs % 2 == 0 "divs must be an even number when shooting=true"
prob = ODEProblem(dirac!, [0, 1], (r_max, r_max / 2)) prob = ODEProblem(dirac!, [0, 1], (r_max, r_max / 2))
sol = solve(prob, RK4(), p=(κ, p, E, Φ0, W0, B0, A0), saveat=saveat) sol = solve(prob, RK4(), p=(κ, p, E, Φ0, W0, B0, A0), saveat=Δr)
wf_right = reverse(hcat(sol.u...); dims=2) wf_right = reverse(hcat(sol.u...); dims=2)
r_max = r_max / 2 # for the next step r_max = r_max / 2 # for the next step
end end
prob = ODEProblem(dirac!, [0, 1], (0, r_max)) prob = ODEProblem(dirac!, [0, 1], (0, r_max))
sol = solve(prob, RK4(), p=(κ, p, E, Φ0, W0, B0, A0), saveat=saveat) sol = solve(prob, RK4(), p=(κ, p, E, Φ0, W0, B0, A0), saveat=Δr)
wf = hcat(sol.u...) wf = hcat(sol.u...)
if shooting # join two segments if shooting # join two segments
@ -52,8 +52,8 @@ function solveNucleonWf(κ, p, E, Φ0, W0, B0, A0, r_max, divs; shooting=true, n
end end
if normalize if normalize
norm = sum(wf .* wf) * r_max / divs # integration by Reimann sum norm = sum(wf .* wf) * Δr # integration by Reimann sum
wf = wf ./ sqrt(norm * 2) # WHY FACTOR OF 2? wf = wf ./ sqrt(norm)
end end
return wf return wf