Reorganize shooting method (not optional anymore)

nucleons.jl

@@ -56,27 +56,28 @@ 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))
+    # 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]
+
+    # 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)    
-        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...)
-
-    if shooting # join two segments
-        u1 = wf[:, end]
+    # 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"
@@ -85,8 +86,7 @@ function solveNucleonWf(κ, p::Bool, E, s::system; shooting=true, normalize=true
         proj = only(u1' * u2) / norm(u2)^2
         wf_right .*= proj
     end
-        wf = hcat(wf[:, 1:(end - 1)], wf_right)
-    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