Functions to compute Clebsch-Gordon basis for addition of two angular momentum states.  Support for computing up to three columns of such basis elements is supported.
This is used then applied to compute the (Subscript[l, 1 ])= 2 ⊕ (Subscript[l, 2 ])= 1 coefficients for phy1520 problem set 6, problem 2.

ClearAll[lower, lowerKet, lowerTwo, sumByList, lowerAll, formatKetTot, formatKetInd, formatKetList, normalPair, ketGreaterInd,ketGreater,sortKets, normalTriplet, dotOne, innerProd, ket, a, d, l, m, mmin, in, r1, r2, c, scale, k]

(*Compute the lowering coefficient √(l+m)(l-m+1)
input ket stored as a list {l,m}
 *)
lower[ket_, mmin_] := Module[{l,m},
{l,m} = ket ;
Boole[m ≥ mmin] 
Sqrt[(l+m)(l-m+1)] 
] 

(*|l,m> -> |l,m-1>
input ket stored as a list {l,m}
 *)
lowerKet[ket_, d_ :  1] := {ket[[1]], ket[[2]]-d} 

(*
|l1,m1 ; l2, m2> -> ((SubMinus[L])^1)|l1,m1 ; l2, m2> + ((SubMinus[L])^2)|l1,m1 ; l2, m2>

Input ket stored as a set of two kets and a numeric coefficent c:

{c, k1, k2}, where k1 = {l1, m1}, k2 = {l2,m2}

output is a list of the form:
{{c', k1', k2'}, {c'', k1'', k2''}}

List input param mmin are the minimum values of m1, m2 respectively.
*)
lowerTwo[in_, mmin_List] := Module[{r1, r2, c, k1, k2}, {c, {k1,k2}} = in ;
r1 = { c lower[k1, mmin[[1]]],{lowerKet[k1], k2}} ;
r2 ={ c lower[ k2, mmin[[2]]] ,{k1, lowerKet[k2]}} ;
{r1, r2}
]

(*Given a list
l = {{v_1,k_1}, {v_2, k_2}, ...} 
find and sum the v_i values for any k_j == k
*)
sumByList[l_,k_] := #[[1]] &/@ Select[l, #[[2]] == k &] // Total;

(*comparison function for two individual kets each specified by {l, m}*)
ketGreaterInd[k1_, k2_] := ( (k1[[1]] > k2[[1]]) ||
(k1[[1]] == k2[[1]] && k1[[2]] > k2[[2]]) ) ;

(*comparison function used to sort lists of {c, ket}, where ket = {{l1,m1}, {l2,m2}}.  This compares two kets specified by k1 = {c,ket}, k2 = {c', ket'} by comparing ket against. ket'.
 Only the basis portion (not the coefficient) is used in the comparision, since it is assumed that same-basekets have been merged when generated (i.e. as done in lowerAll). *)
ketGreater[k1_, k2_] := Module[{t1, t2, l1, l2, f1, f2},
t1 = k1 // Last ;
t2 = k2 // Last ;
f1 = t1 // First ;
f2 = t2 // First ;
l1 = t1 // Last ;
l2 = t2 // Last ;
ketGreaterInd[f1, f2] || ((f1 == f2) && ketGreaterInd[l1, l2])
] ;

(*sort a list of coefficient-ket pairs by their base kets*)
sortKets[a_] :=Sort[a, ketGreater[#1, #2]&] ;

(*
apply lowerTwo to all elements in a list (sum) of kets, weed out any zeros and merge (sum) the coefficients for any kets that are duplicated.
Each coefficient is divided by scale.
 *)
lowerAll[in_, mmin_List, scale_] := Module[{lowered, flattened, nonzero, kets, sumAndScale},
lowered = lowerTwo[#, mmin] & /@ in ;
flattened = Flatten[lowered, 1] ;
nonzero = Select[flattened, #[[1]] ≠ 0 & ] ;
kets = (#[[2]] &/@ nonzero) // DeleteDuplicates ;
sumAndScale ={sumByList[nonzero, #]/scale // Simplify, #} &/@ kets ;
sortKets[sumAndScale]
]




(*A display method for a total angular momentum ket {l,m}*)
formatKetTot[k_] := Row[{"\[VerticalSeparator]",k[[1]],",", k[[2]],"⟩"}] ;
(*A display method for an individual two particle ket {c, {{l1, m1},{l2,m2}}}
*)
formatKetInd[k_] := Module[{c,k1, k2},
c = k // First ;
{k1, k2} = k// Last ;Row[{c,"\[VerticalSeparator]",k1[[1]],",", k1[[2]], ";", k2[[1]],",", k2[[2]],"⟩"}] 
]
formatKetList[k_] := formatKetInd[#] &/@ k ;


(*Given a unit normalized ket composed of two weighted states, compute a state normal to that*)
normalPair[k_] := Module[{c1, k1, c2,k2},
{c1, k1} = First[k] ;
{c2, k2} = Last[k] ;
{{c2, k1}, {-c1, k2}}
]

(*Given two orthonomal kets with three components each, compute a state normal to that.  Assumes the respective basis lists are co-sorted (which is the case after lowerAll). *)
normalTriplet[k1_, k2_] := Module[{k, t, c1,c2, cr},
k = (#[[2]] &/@ k1) ;
t = (#[[2]] &/@ k2) ;
Assert[k == t] ;

c1 = (#[[1]] &/@ k1) ;
c2 = (#[[1]] &/@ k2) ;
cr = Cross[c1, c2] ;
{cr[[#]],k[[#]]} &/@ Range[3]
]

(* for two kets each with a nonzero projection onto state k,take the dot products of those projections (assuming real coefficients) *)
dotOne[k1_, k2_, k_] := Module[{c1, c2},
c1 = Select[k1, #[[2]]==k &] ;
c2 = Select[k2, #[[2]]==k &] ;
c1 = #[[1]] &/@ c1 ;
c2 = #[[1]] &/@ c2 ;
c1.c2
]

(*compute the inner product of two kets in the |l1,m1;l2,m2> basis. *)
innerProd[k1_, k2_] := Module[{t1, t2, k, dotlist},
t1 = (#[[2]] &/@ k1) ;
t2 = (#[[2]] &/@ k2) ;
k = Intersection[ t1, t2] ;
If [ k ≠ {},
Total[dotOne[k1, k2, #] &/@ k],
0 ]
]

Compute the Clebsch-Gordon coefficients for (Subscript[l, 1 ])= 2, (Subscript[l, 2])= 1.

ClearAll[l1, l2, lt, k1, kt1, k2, kt2, k3, kt3, n, minM]
l1 = 2 ;
l2 = 1 ;
lt = l1+l2 ;
n = 2 lt + 1 ;
minM = -{l1, l2} ;


k1 = Range[n]  ; (*dummy array to reserve storage*)
kt1 = Table[{lt,lt-i+1}, {i,1,n}] ;
k1[[1]] = {{1, {{l1,l1},{l2, l2}}}} ;
For[i=2,i <= n, i++,k1[[i]] = lowerAll[k1[[i-1]], minM + {1,1}, lower[kt1[[i-1]], -lt]] ] ;
(*lowerAll[k1[[n]], minM, 1]*)


k2 = Range[n-2] ; (*dummy array to reserve storage*)
kt2 = Table[{lt-1,lt-1-i+1}, {i,1,n-2}] ;
k2[[1]] = normalPair[ k1[[2]] ] ;
For[i=2,i <= n-2, i++,k2[[i]] = lowerAll[k2[[i-1]], minM, lower[kt2[[i-1]], -lt+1]] ] ;
(*lowerAll[k2[[n-2]], minM, 1]*)

k3 = Range[n-4] ; (*dummy array to reserve storage*)
kt3 = Table[{lt-2,lt-2-i+1}, {i,1,n-4}] ;
k3[[1]] = normalTriplet[ k1[[3]],k2[[2]] ] ;
For[i=2,i <= n-4, i++,k3[[i]] = lowerAll[k3[[i-1]], minM, lower[kt3[[i-1]], -lt+2]] ] ;
(*lowerAll[k3[[n-4]], minM, 1]*)

(* Note: Total was an easy way to display the results, but it looses the internal sort order of the basekets *)
Print["(Subscript[l, tot]) = 3: "]
For[ i = 1, i ≤ n, i++,Row[{formatKetTot[kt1[[i]]], " = ", formatKetList[ k1[[i]] ] // Total }] // Print]

Print["(Subscript[l, tot]) = 2: "]
For[ i = 1, i ≤ n-2, i++,Row[{formatKetTot[kt2[[i]]], " = ", formatKetList[ k2[[i]] ] // Total }] // Print]

Print["(Subscript[l, tot]) = 1: "]
For[ i = 1, i ≤ n-4, i++,Row[{formatKetTot[kt3[[i]]], " = ", formatKetList[ k3[[i]] ] // Total }] // Print]


(*check that all pairs are normal*)
On[Assert] ;
For[ i = 1, i ≤ n, i++,
For[ j = 1, j ≤ n-2, j++,
innerProd[k1[[i]], k2[[j]]] == 0 // Assert
] ;
For[ j = 1, j ≤ n-4, j++,
innerProd[k1[[i]], k3[[j]]] == 0 // Assert
]
]
For[ i = 1, i ≤ n-2, i++,
For[ j = 1, j ≤ n-4, j++,
innerProd[k2[[i]], k3[[j]]] == 0 // Assert
]
]

(* check that all results are unit kets *)
For[ i = 1, i ≤ n, i++,
innerProd[k1[[i]], k1[[i]]] == 1 // Assert
]
For[ i = 1, i ≤ n-2, i++,
innerProd[k2[[i]], k2[[i]]] == 1 // Assert
] ;
For[ i = 1, i ≤ n-4, i++,
innerProd[k3[[i]], k3[[i]]] == 1 // Assert
]

(*lower[{1,1}, -1]
lower[{1,0}, -1]
lower[{1/2,1/2}, -1/2]
lower[{3/2,3/2}, -3/2]
lower[{3/2,-1/2}, -3/2]
lower[{3/2,-3/2}, -3/2]

a = lowerTwo[{1, {{1,1},{1/2, 1/2}}}, {-1, -1/2}]
formatKetTot[{1,2}]
formatKetInd[{Sqrt[2], {{1,1}, {1/2, -1/2}}}]
formatKetInd[#] & /@ a

a = {{(Sqrt[(2/5)]),{{2,0},{1,-1}}},{2 (Sqrt[(2/15)]),{{2,1},{1,0}}},{(1/(Sqrt[15])),{{2,2},{1,-1}}}, {(Sqrt[(2/5)]),{{2,0},{2,1}}} , {(Sqrt[(2/5)]),{{2,0},{1,1}}}} 

ketGreater[ {(Sqrt[(2/5)]),{{2,0},{1,1}}},{(Sqrt[(2/5)]),{{2,0},{2,1}}}]
ketGreater[ {(Sqrt[(2/5)]),{{2,0},{2,1}}},{(Sqrt[(2/5)]),{{2,0},{1,-1}}} ]
ketGreaterInd[ {2,1}, {1,-1}]

On[Assert];
Off[Assert];



a = {{(2/(Sqrt[15])),{{2,2},{1,-1}}},{2 (Sqrt[(2/15)]),{{2,1},{1,0}}},{(Sqrt[(2/5)]),{{2,0},{1,1}}}} 
b = {{(1/(Sqrt[3])),{{2,2},{1,-1}}},{(1/(Sqrt[6])),{{3,1},{1,0}}},{-(1/(Sqrt[2])),{{3,0},{1,1}}}}

innerProd[a, b]*)
(*dotOne[a, b, {{2,1},{1,0}}]*)