From 3a9f30e4f5d7fc970a3ce00c8f048cce65a7d000 Mon Sep 17 00:00:00 2001
From: Nuwan Yapa <yapanuwan@gmail.com>
Date: Wed, 12 Mar 2025 20:23:13 -0400
Subject: [PATCH] Gauss-Legendre added

---
 Project.toml  |  1 +
 nucleons.jl   |  4 ++--
 system.jl     | 18 ++++++++++++++++--
 test/Pb208.jl |  4 ++--
 4 files changed, 21 insertions(+), 6 deletions(-)

diff --git a/Project.toml b/Project.toml
index 3f55188..609fd15 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,5 +1,6 @@
 [deps]
 DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
+FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838"
 Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 PolyLog = "85e3b03c-9856-11eb-0374-4dc1f8670e7f"
diff --git a/nucleons.jl b/nucleons.jl
index 2d50817..90e5eec 100644
--- a/nucleons.jl
+++ b/nucleons.jl
@@ -33,8 +33,8 @@ function optimized_dirac_potentials(p::Bool, s::system)
     @. f1s =  M - s.Φ0 - s.W0 - (p - 0.5) * 2 * s.B0 - p * s.A0
     @. f2s = -M + s.Φ0 - s.W0 - (p - 0.5) * 2 * s.B0 - p * s.A0
 
-    f1 = linear_interpolation(s.r_mesh.r, f1s)
-    f2 = linear_interpolation(s.r_mesh.r, f2s)
+    f1 = linear_interpolation(s.r_mesh.r, f1s, extrapolation_bc = Flat())
+    f2 = linear_interpolation(s.r_mesh.r, f2s, extrapolation_bc = Flat())
 
     return (f1, f2)
 end
diff --git a/system.jl b/system.jl
index 3ccf6cb..f63056e 100644
--- a/system.jl
+++ b/system.jl
@@ -1,3 +1,5 @@
+using FastGaussQuadrature
+
 "Tabulates a nucleon spectrum (protons or neutrons) containing κ and occupancy"
 struct spectrum
     κ::Vector{Int}
@@ -22,6 +24,15 @@ function uniform_mesh(r_max::Float64, divs::Int)::mesh
     return mesh(r_max, r, w)
 end
 
+"Create a Gauss-Legendre mesh"
+function gauleg_mesh(r_max::Float64, n::Int)::mesh
+    a, b = (0.0, r_max)
+    nodes, weights = gausslegendre(n)
+    transformed_nodes = (b - a) / 2 .* nodes .+ (b + a) / 2
+    transformed_weights = (b - a) / 2 .* weights
+    return mesh(r_max, transformed_nodes, transformed_weights)
+end
+
 "Defines a nuclear system containing relevant parameters and meson/nucleon densities"
 mutable struct system
     Z::Int
@@ -42,11 +53,14 @@ mutable struct system
     ρ_sn::Vector{Float64}
     ρ_vn::Vector{Float64}
 
-    "Initialize an unsolved system"
-    system(Z, N, r_max, divs) = new(Z, N, uniform_mesh(r_max, divs), unfilled_spectrum(), unfilled_spectrum(), [zeros(1 + divs) for _ in 1:8]...)
+    "Initialize an unsolved system with a uniform r-mesh"
+    system(Z::Int, N::Int, r_max::Float64, divs::Int) = new(Z, N, uniform_mesh(r_max, divs), unfilled_spectrum(), unfilled_spectrum(), [zeros(1 + divs) for _ in 1:8]...)
 
     "Dummy struct to define the mesh"
     system(r_max, divs) = system(0, 0, r_max, divs)
+
+    "Initialize an unsolved system with a Gauss-Legendre r-mesh"
+    system(Z::Int, N::Int, mesh_size::Int, r_max::Float64) = new(Z, N, gauleg_mesh(r_max, mesh_size), unfilled_spectrum(), unfilled_spectrum(), [zeros(mesh_size) for _ in 1:8]...)
 end
 
 "Get mass number of nucleus"
diff --git a/test/Pb208.jl b/test/Pb208.jl
index cae50b9..b733adc 100644
--- a/test/Pb208.jl
+++ b/test/Pb208.jl
@@ -3,8 +3,8 @@ include("../NuclearRMF.jl")
 Z = 82
 N = 126
 r_max = 20.0
-divs = 400
+mesh_points = 400
 
-s = system(Z, N, r_max, divs)
+s = system(Z, N, mesh_points, r_max)
 
 solve_system!(s; monitor_plot=true)