diff --git a/nucleons.jl b/nucleons.jl
index 52eb893..008f467 100644
--- a/nucleons.jl
+++ b/nucleons.jl
@@ -138,29 +138,99 @@ function findAllOrbitals(p::Bool, s::system, E_min=850.0, E_max=938.0)
     return (κs, Es)
 end
 
-"For a given list of κ values with corresponding energies, attempt to fill Z_or_N lowest lying orbitals and return the spectrum"
-function fillNucleons(Z_or_N::Int, κs, Es)::spectrum
-    sort_i = sortperm(Es)
-
-    occ = zeros(Int, length(κs))
-
-    for i in sort_i
-        if Z_or_N ≤ 0; break; end;
-        max_occ = 2 * j_κ(κs[i]) + 1
-        occ[i] = min(max_occ, Z_or_N)
-        Z_or_N -= occ[i]
-    end
-
-    @assert Z_or_N == 0 "All orbitals could not be filled"
-    return spectrum(κs, Es, occ)
-end
-
 "Total angular momentum j for a given κ value"
 j_κ(κ::Int) = abs(κ) - 1/2
 
 "Orbital angular momentum l for a given κ value"
 l_κ(κ::Int) = abs(κ) - (κ < 0) # since true = 1 and false = 0
 
+"Pair degeneracy Ω = (2j+1)/2 for a given κ (equals |κ|)"
+Ω_κ(κ::Int) = abs(κ)
+
+"Constant-G monopole pairing strength: G = 20/A MeV [Bohr-Mottelson Vol. 1]"
+pairing_G(s::system)::Float64 = 20.0 / A(s)
+
+"Energy window (MeV) around the Fermi level within which BCS pairing is active;
+    orbitals below (above) the window are treated as fully filled (empty) core states.
+    This cut-off regularizes the otherwise ill-defined constant-G gap equation."
+const E_pair_window = 5.0
+
+"Sharp (non-pairing) filling of the Z_or_N lowest orbitals — kept for diagnostics."
+function fillNucleons(Z_or_N::Int, κs, Es)::spectrum
+    occ = zeros(Float64, length(κs))
+    rem = Z_or_N
+    for i in sortperm(Es)
+        rem ≤ 0 && break
+        take = min(2 * Ω_κ(κs[i]), rem)
+        occ[i] = take
+        rem -= take
+    end
+    @assert rem == 0 "All orbitals could not be filled"
+    return spectrum(κs, Es, occ, 0.0, 0.0)
+end
+
+"Solve BCS equations for a constant-G monopole pairing force. Returns (λ, Δ, occ) where
+    occ[i] = 2 Ω_i v²_i sums to N_target, and the canonical-basis occupations are
+    v²_i = ½(1 - (ε_i − λ)/√((ε_i − λ)² + Δ²)) [Vretenar et al. 2005, Eq. 50].
+    Orbitals outside ±E_window of the Fermi level are frozen as fully filled / empty."
+function solveBCS(κs::Vector{Int}, Es::Vector{Float64}, G::Float64, N_target::Int;
+                  E_window::Float64=E_pair_window, tol::Float64=1e-10, max_iter::Int=200)
+    isempty(κs) && return (0.0, 0.0, Float64[])
+
+    Ω = Float64.(Ω_κ.(κs))
+    occ = zeros(Float64, length(κs))
+
+    # Sharp-filling estimate of the Fermi level
+    idx = sortperm(Es)
+    cum = cumsum(2 .* Ω[idx])
+    i_HOMO = findfirst(≥(N_target), cum)
+    i_HOMO === nothing && error("Not enough orbitals for N_target=$N_target particles")
+    λ_F = Es[idx[i_HOMO]]
+
+    # Core (frozen filled), active (paired), and above-window (empty)
+    core   = findall(e -> e < λ_F - E_window, Es)
+    active = findall(e -> abs(e - λ_F) ≤ E_window, Es)
+    occ[core] .= 2 .* Ω[core]
+    N_active = N_target - sum(occ[core])
+    max_active = 2 * sum(Ω[active])
+
+    # Closed / trivial case: no pairing to solve in the active window
+    if N_active ≤ 0 || N_active ≥ max_active
+        rem = N_active
+        for i in sort(active, by = i -> Es[i])
+            take = clamp(rem, 0.0, 2 * Ω[i])
+            occ[i] = take
+            rem -= take
+        end
+        return (λ_F, 0.0, occ)
+    end
+
+    Es_a, Ω_a = Es[active], Ω[active]
+    v²(λ, Δ)  = @. 0.5 * (1 - (Es_a - λ) / sqrt((Es_a - λ)^2 + Δ^2))
+    N_of(λ, Δ) = sum(2 .* Ω_a .* v²(λ, Δ))
+    gap_f(λ, Δ) = G/2 * sum(@. Ω_a / sqrt((Es_a - λ)^2 + Δ^2)) - 1
+
+    λ, Δ = λ_F, 1.0
+    λ_lo, λ_hi = minimum(Es_a) - 10.0, maximum(Es_a) + 10.0
+    for _ in 1:max_iter
+        λ_new = bisection(λ -> N_of(λ, Δ) - N_active, λ_lo, λ_hi; tol=tol)
+        Δ_new = gap_f(λ_new, 0.0) < 0 ? 0.0 :
+                bisection(Δ -> gap_f(λ_new, Δ), 0.0, 50.0; tol=tol)
+        converged = abs(λ_new - λ) < tol && abs(Δ_new - Δ) < tol
+        λ, Δ = λ_new, Δ_new
+        converged && break
+    end
+
+    occ[active] .= 2 .* Ω_a .* v²(λ, Δ)
+    return (λ, Δ, occ)
+end
+
+"Fill orbitals via BCS with constant-G monopole pairing and return the spectrum."
+function fillNucleonsBCS(N_target::Int, κs, Es, G::Float64)::spectrum
+    (λ, Δ, occ) = solveBCS(κs, Es, G, N_target)
+    return spectrum(κs, Es, occ, Δ, λ)
+end
+
 "Calculate scalar and vector densities of a nucleon species on [0,r_max] divided into (divs+1) points and returns them as vectors (ρ_s, ρ_v) where
     the arrays κs, Es, occs tabulate the energies and occupation numbers corresponding to each κ,
     the other parameters are defined above"
@@ -171,6 +241,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)
+        occ < 1e-10 && continue
         wf = solveNucleonWf(κ, p, E, s; normalize=true)
         wf2 = wf .* wf
         ρr2 += (occ / (4 * pi)) * wf2 # 2j+1 factor is accounted in the occupancy number
@@ -190,7 +261,7 @@ end
     the other parameters are defined above"
 function solveNucleonDensity!(p::Bool, s::system, E_min=850.0, E_max=938.0)
     κs, Es = findAllOrbitals(p, s, E_min, E_max)
-    spec = fillNucleons(Z_or_N(s, p), κs, Es)
+    spec = fillNucleonsBCS(Z_or_N(s, p), κs, Es, pairing_G(s))
     if p
         s.p_spectrum = spec
     else
@@ -206,5 +277,12 @@ function solveNucleonDensity!(p::Bool, s::system, E_min=850.0, E_max=938.0)
     end
 end
 
-"Total energy of filled nucleons in the system"
-nucleon_E(s::system) = sum(s.p_spectrum.occ .* (s.p_spectrum.E .- M_p)) + sum(s.n_spectrum.occ .* (s.n_spectrum.E .- M_n))
+"Total energy of filled nucleons in the system, including the BCS pairing energy
+    E_pair = −Δ²/G for each species [Vretenar et al. 2005, Eq. 34 with constant-G pp-force]."
+function nucleon_E(s::system)
+    G = pairing_G(s)
+    E_sp = sum(s.p_spectrum.occ .* (s.p_spectrum.E .- M_p)) +
+           sum(s.n_spectrum.occ .* (s.n_spectrum.E .- M_n))
+    E_pair = -(s.p_spectrum.Δ^2 + s.n_spectrum.Δ^2) / G
+    return E_sp + E_pair
+end
diff --git a/system.jl b/system.jl
index 0298438..83dda83 100644
--- a/system.jl
+++ b/system.jl
@@ -25,15 +25,19 @@ struct parameters
     end
 end
 
-"Tabulates a nucleon spectrum (protons or neutrons) containing κ and occupancy"
+"Tabulates a nucleon spectrum (protons or neutrons) containing κ, single-particle energies,
+    occupation numbers (2Ω·v² per κ, possibly fractional under BCS), the pairing gap Δ,
+    and the chemical potential λ."
 struct spectrum
     κ::Vector{Int}
     E::Vector{Float64}
-    occ::Vector{Int}
+    occ::Vector{Float64}
+    Δ::Float64
+    λ::Float64
 end
 
 "Initializes an unfilled spectrum"
-unfilled_spectrum() = spectrum(Int[], Float64[], Int[])
+unfilled_spectrum() = spectrum(Int[], Float64[], Float64[], 0.0, 0.0)
 
 "Defines a nuclear system containing relevant parameters and meson/nucleon densities"
 mutable struct system