NF = 3
NB = 0

opp1cr = NewOperator(NF,NB, {{ 0,1},{ 1,1},{ 2,1}})
opp1an = NewOperator(NF,NB, {{-0,1},{-1,1},{-2,1}})
opp2 = opp1cr * opp1an
opp3 = opp1an * opp2
opp4 = opp1cr * opp3

opp = 1 + opp1cr + opp1an + opp2 + opp3 + opp4

-- The operatore contains strings of creation and annihilation operators of length 0 to 4
-- l=0  1 term
-- l=1  6 terms (3 creation, 3 annihilation)
-- l=2  9 terms (3 creation * 3 annihilation)
-- l=3  9 terms (cr an an) note that an0 an1 = - an1 an0
-- l=4  9 terms 
-- the full operator is
print("================ opp =============")
print(opp)

-- We only keep terms that at least act once on orbital 1
print("================ PartialOperator  {1} include =============")
print(PartialOperator(opp,{1},"include"))

-- We only keep terms that at least act once on orbital 0 or at orbital 1
print("================ PartialOperator  {0,2} include =============")
print(PartialOperator(opp,{0,2},"include"))

-- We remove all terms that act on orbital 1
print("================ PartialOperator  {1} exclude =============")
print(PartialOperator(opp,{1},"exclude"))

-- We exclude all terms that act on orbital 0 and we exclude all terms that act on orbital 2
print("================ PartialOperator  {0,2} exclude =============")
print(PartialOperator(opp,{0,2},"exclude"))

-- We only keep the terms that keep the occupation of orbital 1 conserved
print("================ PartialOperator  {1} conserve =============")
print(PartialOperator(opp,{1},"conserve"))

-- We only keep the terms that keep the occupation of orbital 0 conserved and that keep
-- the occupation of orbital 2 conserved
print("================ PartialOperator  {0,2} conserve =============")
print(PartialOperator(opp,{0,2},"conserve"))

-- We only keep those terms that keep the sum of the occupation of orbital 0 and 2 conserved
print("================ PartialOperator  {{0,2}} conserve =============")
print(PartialOperator(opp,{{0,2}},"conserve"))