"""Build your own shell model -- step 1: the valence space.

A shell-model calculation starts by throwing almost everything away. For a
nucleus like 18O we declare the 16O core FROZEN -- its 8 protons and 8 neutrons
never move -- and keep only the 2 leftover neutrons. Those two live in the
handful of orbitals just above the core: the sd shell, made of the orbitals
d5/2, s1/2, d3/2 (spectroscopic names: letter = orbital angular momentum
l = 2, 0, 2; subscript = total angular momentum j).

Each orbital with total angular momentum j holds 2j+1 magnetic substates,
m = -j ... +j. Listing every (orbital, m) combination for ONE species gives the
single-particle m-states -- the "seats" a valence neutron can occupy. That
list is the whole universe of this calculation.

The only other ingredient this step needs: a Slater determinant (the
antisymmetric many-body state built by filling seats, one particle per seat --
the Pauli principle as bookkeeping) is labeled by WHICH seats are filled.
Because the interaction conserves total Jz, we only need determinants whose
m-values add up to the Jz we ask for. Counting them for 18O is the first
real number the tutorial produces.
"""
from itertools import combinations

# label -> (name, l, 2j).  The numbering matches the interaction file coming
# in step 3 (usdb.int uses 1 = d3/2, 2 = d5/2, 3 = s1/2).
ORBITS = {1: ('d3/2', 2, 3), 2: ('d5/2', 2, 5), 3: ('s1/2', 0, 1)}

# every (orbital, m) seat for one species, in a fixed order
SEATS = []
for label, (name, l, j2) in ORBITS.items():
    for m2 in range(-j2, j2 + 1, 2):          # 2m runs -2j ... +2j in steps of 2
        SEATS.append((label, name, l, j2, m2))

print("the sd-shell single-particle m-states (one species):")
print(f"{'idx':>3} {'orbital':>8} {'l':>3} {'j':>5} {'m':>6}")
for i, (label, name, l, j2, m2) in enumerate(SEATS):
    print(f"{i:>3} {name:>8} {l:>3} {j2}/2 {m2:>4}/2")
print(f"total seats: {len(SEATS)}  (= (3+1) + (1+1) + (5+1) = 2j+1 summed)")

# 18O: two valence neutrons, total Jz = 0 -> pick 2 distinct seats whose
# m-values cancel. combinations() enforces the Pauli principle for free.
count = 0
for (s1, s2) in combinations(range(len(SEATS)), 2):
    if SEATS[s1][4] + SEATS[s2][4] == 0:
        count += 1
print(f"\n18O basis: 2 neutrons, Jz = 0  ->  {count} Slater determinants")
print("fourteen states. The whole 18O problem is a 14 x 14 matrix.")