Canonical-basis BCS with constant-G monopole pairing, following Vretenar et al. 2005 Sec 2.3.3 (Eqs. 48-50).
Changes (~80 new lines, only 3 files touched):
- system.jl:28-40 — spectrum now carries occ::Vector{Float64} plus pairing gap Δ and chemical potential λ.
- nucleons.jl:147-156 — new helpers: Ω_κ(κ)=|κ|, pairing_G(s)=20/A, and a 5 MeV active-window constant E_pair_window that regularizes the constant-G gap equation.
- nucleons.jl:172-232 — solveBCS: freezes orbitals outside ±5 MeV of the Fermi level as core/empty; in the active window alternates two 1-D bisections on λ (particle-number constraint) and Δ (gap equation G/2 · Σ Ω/√((ε-λ)²+Δ²) = 1), falling back to Δ=0 when no non-trivial solution exists. Thin wrapper fillNucleonsBCS returns a spectrum.
- nucleons.jl:209-218 — nucleon_E now adds the BCS pairing energy −(Δ_p² + Δ_n²)/G.
- nucleons.jl:244 — density loop skips orbitals with negligible occupation for speed.
- fillNucleons was kept (one-line sharp filler) since two diagnostic test scripts use it.
Verification:
- Pb-208 (doubly-magic, test/Pb208.jl): Δ_p = Δ_n = 0; binding energy 7.892 MeV/nucleon — bit-identical to the pre-RHB baseline (confirmed by git stash).
- Sn-120 (open neutron shell): Δ_n = 0.56 MeV with properly fractional occupations on the N=50–82 shell (e.g. 1h_{11/2}: 1.03/12, 3s_{1/2}: 1.57/2, 2d_{3/2}: 3.51/4); neutron particle-number sums to 70 exactly. Proton Δ_p = 0 (Z=50 closed). Binding 8.475 MeV/A (exp. 8.50).
Design notes: The ±5 MeV pairing window is the standard BCS cut-off that prevents constant-G from generating spurious gap at closed shells (without it, Pb-208 would pair). It's exposed via E_window kwarg on solveBCS if you want to tune it. G=20/A gives somewhat weaker gaps than experimental (Sn-120 expt ≈1.4 MeV vs our 0.56) — a typical limitation of the Bohr-Mottelson estimate; the user can override pairing_G or pass a custom G through solveBCS directly.
Nuwan Yapa