From f77e5dfed8b57287108844db8ca39845a0998684 Mon Sep 17 00:00:00 2001
From: Nuwan Yapa <yapanuwan@gmail.com>
Date: Tue, 16 Jul 2024 18:29:35 -0400
Subject: [PATCH] Fixed unstability of wave functions

---
 dirac.jl         | 19 +++++++++++++------
 test/Pb208_wf.jl | 10 +++++-----
 2 files changed, 18 insertions(+), 11 deletions(-)

diff --git a/dirac.jl b/dirac.jl
index 2317cf0..b2761a3 100644
--- a/dirac.jl
+++ b/dirac.jl
@@ -26,18 +26,18 @@ end
 "Solve the Dirac equation and return the wave function u(r)=[g(r), f(r)] where
     divs is the number of mesh divisions if the solution should be discretized as a 2×(1+divs) matrix (keep divs=0 to obtain an interpolating function),
     the other parameters are the same from dirac!(...)."
-function solveWf(κ, p, E, Φ0, W0, B0, A0, r_max, divs::Int=0)
-    prob = ODEProblem(dirac!, [0, 1], (0, r_max))
-    sol = solve(prob, RK4(), p=(κ, p, E, Φ0, W0, B0, A0), saveat=(divs == 0 ? [] : r_max/divs))
+function solveWf(κ, p, E, Φ0, W0, B0, A0, r_max, divs::Int=0; dtype=BigFloat, algo=Feagin12())
+    prob = ODEProblem(dirac!, convert.(dtype, [0, 1]), (0, r_max))
+    sol = solve(prob, algo, p=(κ, p, E, Φ0, W0, B0, A0), saveat=(divs == 0 ? [] : r_max/divs))
     return divs == 0 ? sol : hcat(sol.u...)
 end
 
 "Solve the Dirac equation and return g(r=r_max) where
     r_max is the outer boundary in fm,
     the other parameters are the same from dirac!(...)."
-function boundaryValue(κ, p, E, Φ0, W0, B0, A0, r_max)
-    prob = ODEProblem(dirac!, [0, 1], (0, r_max))
-    sol = solve(prob, RK4(), p=(κ, p, E, Φ0, W0, B0, A0), saveat=[r_max], save_idxs=[1])
+function boundaryValue(κ, p, E, Φ0, W0, B0, A0, r_max; dtype=Float64, algo=RK4())
+    prob = ODEProblem(dirac!, convert.(dtype, [0, 1]), (0, r_max))
+    sol = solve(prob, algo, p=(κ, p, E, Φ0, W0, B0, A0), saveat=[r_max], save_idxs=[1])
     return sol[1, 1]
 end
 
@@ -48,6 +48,13 @@ function findEs(κ, p, Φ0, W0, B0, A0, r_max, E_min=0, E_max=(p ? M_p : M_n))
     return find_zeros(f, (E_min, E_max))
 end
 
+"Find more precise bound energies for a given list of Es where
+    the other parameters are the same from dirac!(...)."
+function refineEs(κ, p, Φ0, W0, B0, A0, r_max, Es)
+    f(E) = boundaryValue(κ, p, E, Φ0, W0, B0, A0, r_max; dtype=BigFloat, algo=Feagin12())
+    return [find_zero(f, E) for E in Es]
+end
+
 "Find all orbitals and return two lists containing κ values and corresponding energies for a single species where
     the other parameters are defined above"
 function findAllOrbitals(p, Φ0, W0, B0, A0, r_max, E_min=0, E_max=(p ? M_p : M_n))
diff --git a/test/Pb208_wf.jl b/test/Pb208_wf.jl
index e87c3c9..1d10236 100644
--- a/test/Pb208_wf.jl
+++ b/test/Pb208_wf.jl
@@ -21,14 +21,14 @@ r_max = maximum(xs)
 E_min = 880
 E_max = 939
 
-boundEs = findEs(κ, p, S_interp, V_interp, R_interp, A_interp, r_max, E_min, E_max)
-groundE = minimum(boundEs)
+approxE = findEs(κ, p, S_interp, V_interp, R_interp, A_interp, r_max, E_min, E_max) |> minimum
+groundE = refineEs(κ, p, S_interp, V_interp, R_interp, A_interp, r_max, [approxE])[1]
+
 println("ground state E = $groundE")
 
-plot_r_max = r_max * 0.75
 divs = 50
-wf = solveWf(κ, p, groundE, S_interp, V_interp, R_interp, A_interp, plot_r_max, divs)
-rs = collect(0: plot_r_max/divs : plot_r_max)
+wf = solveWf(κ, p, groundE, S_interp, V_interp, R_interp, A_interp, r_max, divs)
+rs = range(0, r_max, length=divs+1)
 gs = wf[1, :]
 fs = wf[2, :]