From 0af47ef1613790711d989fddd9ce72552104941a Mon Sep 17 00:00:00 2001
From: Nuwan Yapa <yapanuwan@gmail.com>
Date: Wed, 16 Apr 2025 16:08:32 -0400
Subject: [PATCH] Reorganize shooting method (not optional anymore)

---
 nucleons.jl | 54 ++++++++++++++++++++++++++---------------------------
 1 file changed, 27 insertions(+), 27 deletions(-)

diff --git a/nucleons.jl b/nucleons.jl
index 3c3f81d..393637c 100644
--- a/nucleons.jl
+++ b/nucleons.jl
@@ -56,37 +56,37 @@ 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
     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,
     the other parameters are the same from dirac!(...)."
-function solveNucleonWf(κ, p::Bool, E, s::system; shooting=true, normalize=true, algo=Vern9())
+function solveNucleonWf(κ, p::Bool, E, s::system; normalize=true, algo=Vern9())
     (f1, f2) = optimized_dirac_potentials(p, s)
 
-    if shooting
-        @assert s.divs % 2 == 0 "divs must be an even number when shooting=true"
-        prob = ODEProblem(dirac!, asymp_BC(κ, p, E, s.r_max), (s.r_max, s.r_max / 2))
-        sol = solve(prob, algo, p=(κ, E, f1, f2), saveat=Δr(s))
-        wf_right = reverse(hcat(sol.u...); dims=2)
-        next_r_max = s.r_max / 2 # for the next step
-    else
-        next_r_max = s.r_max
-    end
-    
-    prob = ODEProblem(dirac!, init_BC(), (0, next_r_max))
-    sol = solve(prob, algo, p=(κ, E, f1, f2), saveat=Δr(s))
-    wf = hcat(sol.u...)
+    # partitioning
+    mid_idx = s.divs ÷ 2
+    r_mid = rs(s)[mid_idx]
+    left_r = rs(s)[1:mid_idx]
+    right_r = rs(s)[mid_idx:end]
 
-    if shooting # join two segments
-        u1 = wf[:, end]
-        u2 = wf_right[:, 1]
-        if norm(u2) < 1e-10
-            @warn "Right partition too small to rescale, setting to zero"
-            wf_right .= 0.0
-        else
-            proj = only(u1' * u2) / norm(u2)^2
-            wf_right .*= proj
-        end
-        wf = hcat(wf[:, 1:(end - 1)], wf_right)
+    # left partition
+    prob = ODEProblem(dirac!, init_BC(), (0, r_mid))
+    sol = solve(prob, algo, p=(κ, E, f1, f2), saveat=left_r)
+    wf_left = hcat(sol.u...)
+
+    # right partition
+    prob = ODEProblem(dirac!, asymp_BC(κ, p, E, s.r_max), (s.r_max, r_mid))
+    sol = solve(prob, algo, p=(κ, E, f1, f2), saveat=right_r)
+    wf_right = reverse(hcat(sol.u...); dims=2)    
+    
+    # join two segments
+    u1 = wf_left[:, end]
+    u2 = wf_right[:, 1]
+    if norm(u2) < 1e-10
+        @warn "Right partition too small to rescale, setting to zero"
+        wf_right .= 0.0
+    else
+        proj = only(u1' * u2) / norm(u2)^2
+        wf_right .*= proj
     end
+    wf = hcat(wf_left[:, 1:(end - 1)], wf_right)
 
     if normalize
         wf ./= norm(wf) * sqrt(Δr(s)) # integration by Reimann sum
@@ -169,7 +169,7 @@ function calculateNucleonDensity(p::Bool, s::system)::Tuple{Vector{Float64}, Vec
     ρr2 = zeros(2, s.divs + 1) # ρ×r² values
     
     for (κ, E, occ) in zip(κs, Es, occs)
-        wf = solveNucleonWf(κ, p, E, s; shooting=true, normalize=true)
+        wf = solveNucleonWf(κ, p, E, s; normalize=true)
         wf2 = wf .* wf
         ρr2 += (occ / (4 * pi)) * wf2 # 2j+1 factor is accounted in the occupancy number
     end