From 9f1ed7cc6e016cc553ec80c070a6e8f389433e57 Mon Sep 17 00:00:00 2001 From: Nuwan Yapa Date: Tue, 4 Feb 2025 16:56:48 -0500 Subject: [PATCH] Fix silly mistake --- nucleons.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/nucleons.jl b/nucleons.jl index 18d791d..1acdf96 100644 --- a/nucleons.jl +++ b/nucleons.jl @@ -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, the other parameters are the same from dirac!(...)." 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 @assert divs % 2 == 0 "divs must be an even number when shooting=true" 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) r_max = r_max / 2 # for the next step end 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...) if shooting # join two segments @@ -52,8 +52,8 @@ function solveNucleonWf(κ, p, E, Φ0, W0, B0, A0, r_max, divs; shooting=true, n end if normalize - norm = sum(wf .* wf) * r_max / divs # integration by Reimann sum - wf = wf ./ sqrt(norm * 2) # WHY FACTOR OF 2? + norm = sum(wf .* wf) * Δr # integration by Reimann sum + wf = wf ./ sqrt(norm) end return wf