From a143edb49630f52adfcb7b44bae365970c07579b Mon Sep 17 00:00:00 2001
From: Nuwan Yapa <yapanuwan@gmail.com>
Date: Mon, 7 Apr 2025 19:24:42 -0400
Subject: [PATCH] Elegant method to attach two partitions

---
 nucleons.jl | 21 ++++++++++++---------
 1 file changed, 12 insertions(+), 9 deletions(-)

diff --git a/nucleons.jl b/nucleons.jl
index c70a9be..d465af4 100644
--- a/nucleons.jl
+++ b/nucleons.jl
@@ -1,4 +1,4 @@
-using DifferentialEquations, Interpolations
+using LinearAlgebra, DifferentialEquations, Interpolations
 include("bisection.jl")
 include("common.jl")
 include("system.jl")
@@ -73,17 +73,20 @@ function solveNucleonWf(κ, p::Bool, E, s::system; shooting=true, normalize=true
     wf = hcat(sol.u...)
 
     if shooting # join two segments
-        rescale_factor_g = wf[1, end] / wf_right[1, 1]
-        rescale_factor_f = wf[2, end] / wf_right[2, 1]
-        @assert isfinite(rescale_factor_g) && isfinite(rescale_factor_f) "Cannot rescale the right partition"
-        isapprox(rescale_factor_g, rescale_factor_f; rtol=0.03) || @warn "Discontinuity between two partitions"
-        wf_right_rescaled = wf_right .* ((rescale_factor_g + rescale_factor_f) / 2)
-        wf = hcat(wf[:, 1:(end - 1)], wf_right_rescaled)
+        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)
     end
 
     if normalize
-        norm = sum(wf .* wf) * Δr(s) # integration by Reimann sum
-        wf = wf ./ sqrt(norm)
+        wf ./= norm(wf) * Δr(s) # integration by Reimann sum
     end
 
     return wf