(* ::Package:: *)

(* ::Title:: *)
(*Subroutine for 9-J symbol*)


nj[{j1_, j2_, j12_}, {j3_, j4_, j34_}, {j13_, j24_, j_}] := 
  Simplify[Sum[(-1)^(2*g)*(2*g + 1)*sixj[{j1, j34}, {j2, j}, {j12, g}]*
     sixj[{j3, j2}, {j4, g}, {j34, j24}]*sixj[{j13, g}, {j24, j1}, {j, j3}], 
    {g, Max[Abs[j1 - j], Abs[j2 - j34], Abs[j3 - j24]], 
     Min[j1 + j, j2 + j34, j3 + j24], 1 - Mod[2*(j1 + j), 2]/2}]]
 
sixj[{j1_, g1_}, {j2_, g2_}, {j3_, g3_}] := 
  Simplify[Sum[f3[g1, g2, g3, j1, j2, j3, m1, m2, m3], {m1, -g1, g1}, 
    {m2, -g2, g2}, {m3, -g3, g3}]]
 
f3[g1_, g2_, g3_, j1_, j2_, j3_, m1_, m2_, m3_] := 
  (-1)^(g1 + g2 + g3 + j1 - j2)*(-1)^(2*m1 + m3)*
   f1[g1, g2, g3, j3, j1, m1, m2, m3]*f2[g1, g3, j1, j2, j3, m1, m2, m3]*
   Sqrt[1/((2*j1 + 1)*(2*j2 + 1))]
 
f1[g1_, g2_, j3_, j1_, m1_, m2_, m3_] := 
  Cg[{g1, m1}, {g2, -m2}, {j3, m1 - m2}]*
   Cg[{g2, m2}, {g3, -m3}, {j1, m2 - m3}]
 
f1[g1_, g2_, g3_, j3_, j1_, m1_, m2_, m3_] := 
  Cg[{g1, m1}, {g2, -m2}, {j3, j3}]*Cg[{g2, m2}, {g3, -m3}, {j1, m2 - m3}]
 
Cg[{j1_, m1_}, {j2_, m2_}, {j_, m_}] := 
  If[m != m1 + m2 || Abs[m1] > j1 || Abs[m2] > j2 || Abs[m] > j || 
    j > j1 + j2 || j1 > j + j2 || j2 > j + j1, 0, 
   Block[{t}, t = C1[j, j1, j2, m1, m2]; 
     Sign[t]*Sqrt[t^2*C2[j, j1, j2]*C3[j, j1, j2, m, m1, m2]]]]
 
C1[j_, j1_, j2_, m1_, m2_] := 
  Block[{k}, Sum[(-1)^k/
     (k!*(j1 + j2 - j - k)!*(j1 - m1 - k)!*(j2 + m2 - k)!*(j - j2 + m1 + k)!*
       (j - j1 - m2 + k)!), {k, Max[0, j2 - j - m1, j1 - j + m2], 
     Min[j1 + j2 - j, j1 - m1, j2 + m2]}]]
 
C2[j_, j1_, j2_] := 
  ((-j + j1 + j2)!*(j - j1 + j2)!*(j + j1 - j2)!*(1 + 2*j))/(1 + j + j1 + j2)!
 
C3[j_, j1_, j2_, m_, m1_, m2_] := 
  (j - m)!*(j + m)!*(j1 - m1)!*(j1 + m1)!*(j2 - m2)!*(j2 + m2)!
 
f2[g1_, g3_, j1_, j2_, j3_, m1_, m2_, m3_] := 
  Cg[{g3, m3}, {g1, -m1}, {j2, m3 - m1}]*
   Cg[{j1, m3 - m2}, {j2, m1 - m3}, {j3, j3}]
\.1a
