diff --git a/Examples/HydrogenAnion.jl b/Examples/HydrogenAnion.jl
index cce5095..84b59d4 100644
--- a/Examples/HydrogenAnion.jl
+++ b/Examples/HydrogenAnion.jl
@@ -6,9 +6,9 @@ w_list = [ [1, -1, 0], [1, 0, -1], [0, 1, -1] ]
 
 K = [0 0 0; 0 1/2 0; 0 0 1/2]
 J, U = jacobi_transform(masses)
-K_transformedformed = J * K * J'
+K_transformed = J * K * J'
 w_transformed = [U' * w for w in w_list]
 Theortical_value = -0.527751016523
 
-p, Energy, bases = run_simulation(50, :quasirandom, w_transformedformedformed, K_transformedformed)
+p, Energy, bases = run_simulation(50, :quasirandom, w_transformed, K_transformed)
 
diff --git a/Examples/PionNucleon.jl b/Examples/PionNucleon.jl
index 80d12ea..b383f5c 100644
--- a/Examples/PionNucleon.jl
+++ b/Examples/PionNucleon.jl
@@ -1,4 +1,5 @@
 using Plots
+using .FewBodyPhysics
 
 b = 3.9
 S = 41.5
diff --git a/src/FewBodyPhysics.jl b/src/FewBodyPhysics.jl
index 9511fd0..7e3a15c 100644
--- a/src/FewBodyPhysics.jl
+++ b/src/FewBodyPhysics.jl
@@ -1,13 +1,14 @@
 module FewBodyPhysics
 
+export m_pi, m_pi0, m_p, m_n, μ, ħc, m_bare
+export jacobi_transform, generate_A_matrix, transform_list, shift_vectors, generate_weight_vector, transform_coordinates, inverse_transform_coordinates
+export S_elements, S_wave, S_energy, P_elements, pion_nucleon, ComputeEigenSystem, GetMinimumEnergy, OptimizeGlobalParameters
+export corput, halton, run_simulation, run_simulation_nuclear
+
+include("constants.jl")
 include("coordinates.jl")
 include("matrix_elements.jl")
 include("sampling.jl")
-include("constants.jl")
 
-export m_pi, m_pi0, m_p, m_n, μ, ħc, m_bare
-export jacobi_transform, generate_A_matrix, transform_list, shift_vectors, generate_weight_vector, transform_coordinates, inverse_transform_coordinates
-export S_elements, S_wave, S_energy, P_elements, pion_nucleon,compute_eigen_system, get_minimum_energy, optimize_global_parameters
-export corput, halton, run_simulation, run_simulation_nuclear
 
 end
diff --git a/src/constants.jl b/src/constants.jl
index ae3ada8..403dab3 100644
--- a/src/constants.jl
+++ b/src/constants.jl
@@ -1,7 +1,5 @@
 export m_pi, m_pi0, m_p, m_n, μ, ħc, m_bare
 
-# Physical Constants
-
 const m_p   = 938.27       # Proton mass in MeV/c^2
 const m_n   = 939.57       # Neutron mass in MeV/c^2
 const m_pi0 = 134.98       # Neutral pion mass (π⁰) in MeV/c^2
diff --git a/src/matrix_elements.jl b/src/matrix_elements.jl
index 577ce5f..1ad2586 100644
--- a/src/matrix_elements.jl
+++ b/src/matrix_elements.jl
@@ -143,9 +143,8 @@ function P_elements(a::Vector{Float64}, b::Vector{Float64}, A::PositiveDefiniteS
     @assert isposdef(D) "Matrix D must be positive definite."
     R = inv(D)
     M0, trace = S_elements(A.matrix, B.matrix, K)
-    M1 = 0.5 * (a' * R * b) * M0  # Overlap perturbation term
+    M1 = 0.5 * (a' * R * b) * M0 
 
-    # Kinetic energy perturbation term
     kinetic = 6 * trace * M1
     kinetic += (a' * K * b) * M0
     kinetic -= (a' * K * A.matrix * R * b) * M0
@@ -280,8 +279,8 @@ Perform global optimization over a given parameter space to find optimal paramet
 function OptimizeGlobalParameters(ngauss, dim, bmax, masses, params)
     E_list = []
     gaussians = []
-    coords = []
     eigenvectors = []
+    coords = []
 
     global E0S = 0.0
     masses_min = copy(masses)
diff --git a/src/sampling.jl b/src/sampling.jl
index 8f5a4cd..7ab45e6 100644
--- a/src/sampling.jl
+++ b/src/sampling.jl
@@ -2,7 +2,6 @@ using Plots
 
 export corput, halton, run_simulation, run_simulation_nuclear
 
-# Generate the nth element of the van der Corput sequence in base b
 function corput(n, b=3)
     q, bk = 0.0, 1 / b
     while n > 0
@@ -108,7 +107,6 @@ function run_simulation(num_gauss::Int, method::Symbol, w_transformed::Vector{Ve
         error("Invalid method provided!")
     end
 
-    # Plot results
     println("Best convergent numerical value: ", E_list[end])
     p = plot(gaussians, E_list, marker=:circle, label="Numerical result", linewidth=2)
     title!(p, "S-wave Convergence")
@@ -122,9 +120,6 @@ function run_simulation(num_gauss::Int, method::Symbol, w_transformed::Vector{Ve
     return p, E_list[end], bases
 end
 
-
-
-
 """
 Run a nuclear simulation and print the final energy.
 """
diff --git a/test/runtests.jl b/test/runtests.jl
index 93358c8..7b92ece 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,69 +1,71 @@
-using FewBodyPhysics
 using Test
 using LinearAlgebra
+include("../src/FewBodyPhysics.jl")
+
+using .FewBodyPhysics
 
 @testset "FewBodyPhysics.jl" begin
-    @testset "Coordinate Transformations" begin
-        masses = [1.0, 2.0, 3.0]
-        J, U = JacobiTransform(masses)
-        @test size(J) == (2, 3)
-        @test size(U) == (3, 2)
-    
-        bij = [0.5, 0.8]
-        w_list = [generate_weight_vector(3, 1, 2), generate_weight_vector(3, 2, 3)]
-        A = generate_A_matrix(bij, w_list)
-        @test size(A) == (3, 3)
-    
-        α = [1.0, 2.0, 3.0]
-        transformed_list = transform_list(α)
-        @test length(transformed_list) == 3
-        @test transformed_list[1] == [1.0]
-    
-        a = rand(3, 2)
-        b = rand(3, 2)
-        sum_val = shift_vectors(a, b)
-        @test isa(sum_val, Float64)
-    
-        w = generate_weight_vector(3, 1, 2)
-        @test w == [1, -1, 0]
-    
-        r = rand(3)
-        x = transform_coordinates(J, r)
-        r_back = inverse_transform_coordinates(U, x)
-        @test norm(r - r_back) < 1e-10
-    end
 
-    @testset "corput tests" begin
+    @testset "Corput Sequence Tests" begin
+        # Basic tests for base 2
         @test corput(1, 2) ≈ 0.5
         @test corput(2, 2) ≈ 0.25
         @test corput(3, 2) ≈ 0.75
+        @test corput(4, 2) ≈ 0.125
+        @test corput(5, 2) ≈ 0.625
+    
+        # Additional tests for base 3
+        @test corput(1, 3) ≈ 1/3
+        @test corput(2, 3) ≈ 2/3
+        @test corput(3, 3) ≈ 1/9
+        @test corput(4, 3) ≈ 4/9
+        @test corput(5, 3) ≈ 7/9
+    
+        # Additional tests for base 5
+        @test corput(1, 5) ≈ 0.2
+        @test corput(2, 5) ≈ 0.04
+        @test corput(3, 5) ≈ 0.6
+        @test corput(4, 5) ≈ 0.08
+        @test corput(5, 5) ≈ 0.16
+        @test corput(10, 5) ≈ 0.32
+    
     end
     
-    @testset "halton tests" begin
-        @test halton(1, 2) ≈ [0.5, 0.3333333333333333]
-        @test halton(2, 2) ≈ [0.25, 0.6666666666666666]
+    @testset "Halton Sequence Tests" begin
+        # Basic sequence tests for dimension 2
+        @test halton(1, 2) ≈ [0.5, 1/3]
+        @test halton(2, 2) ≈ [0.25, 2/3]
+        @test halton(3, 2) ≈ [0.75, 1/9]
+        @test halton(4, 2) ≈ [0.125, 4/9]
+        @test halton(5, 2) ≈ [0.625, 7/9]
+    
+        # Additional tests for higher dimensions
+        @test halton(1, 3) ≈ [0.5, 1/3, 0.2]
+        @test halton(2, 3) ≈ [0.25, 2/3, 0.4]
+        @test halton(3, 3) ≈ [0.75, 1/9, 0.6]
+        @test halton(4, 3) ≈ [0.125, 4/9, 0.8]
+        
+        
+        # Edge cases: checking very high `n`
+        @test halton(1000, 2) ≈ [corput(1000, 2), corput(1000, 3)]
+        @test halton(5000, 3) ≈ [corput(5000, 2), corput(5000, 3), corput(5000, 5)]
+    
+        # Edge case: requesting only 1 dimension
+        @test halton(1, 1) ≈ [0.5]
+        @test halton(10, 1) ≈ [corput(10, 2)]
+    
+        # Edge case: invalid dimension exceeding base list
         @test_throws AssertionError halton(1, 100)
-    end
     
-    @testset "run_simulation tests" begin
-        w_transformed = [1, 2, 3]
-        K_transformed = [1, 2, 3]
-        p, E_list, bases = run_simulation(15, :quasirandom, w_transformed, K_transformed, false)
-        @test E_list isa Float64
-        @test bases isa Array
-        @test_throws ErrorException run_simulation(15, :invalidmethod, w_transformed, K_transformed, false)
-    end
+        # Edge case: `n = 0` should return zero for any dimension
+        @test halton(0, 2) ≈ [0.0, 0.0]
+        @test halton(0, 5) ≈ [0.0, 0.0, 0.0, 0.0, 0.0]
     
-    @testset "run_simulation_nuclear tests" begin
-        masses = [1, 2, 3]
-        params = [1, 2, 3]
-        E_list, gaussians, eigenvectors, coords, masses = run_simulation_nuclear(2, 2, 5, masses, params)
-        @test E_list isa Array
-        @test gaussians isa Array
-        @test eigenvectors isa Array
-        @test coords isa Array
-        @test masses isa Array
+        # Random large n tests for robustness
+        @test halton(123, 4) ≈ [corput(123, 2), corput(123, 3), corput(123, 5), corput(123, 7)]
+        @test halton(999, 3) ≈ [corput(999, 2), corput(999, 3), corput(999, 5)]
     end
-
+    
     
 end
+