(* :Title: Clebsch-Gordan Coefficients *)
 
(* :Author: Paul C. Abbott  *)
 
 
(* :Summary:
  This package computes the Clebsch Gordan, 3-j coefficients (Wigner),
  and 6-j coefficients (Wigner) for both numeric and symbolic arguments.
*)
 
(* :Context: ClebschGordan` *)
 
 
(* :Package Version: 1.1 *)
 
(* :Copyright: Copyright 1990 Wolfram Research, Inc. *)
 
(* :History:
  Version 1.0 by Stephen Wolfram (Wolfram Research), May 1987.
  Extended to Symbolic arguments by Paul C. Abbott, February 1990.
  Modified for new behaviour of Hypergeometric functions
  	by Paul C. Abbott, November 1991.
  Modified to rationalize real input, February 1992.
*)
 
(* :Keywords:
  Clebsch-Gordan coefficients, Wigner coefficients, 3-j symbols,
  6-j symbols, Racah coefficients, spin statistics, vector addition
  coefficients, spherical harmonics, angular momentum
*)
 
(* :Source:
  Edmonds, A. R.: "Angular momenta in quantum mechanics",
  Princeton University Press, 1974
*)
 
(* :Source:
  Rao, K. R. and Chiu, C. B.: Journal of Physics, A,
  Vol. 22 pp3779-88, 1989
*)
 
(* :Mathematica Version: 2.0 *)
 
Begin["System`"]
 
Unprotect[ ClebschGordan, ThreeJSymbol, SixJSymbol]

Map[ Clear, {ClebschGordan, ThreeJSymbol, SixJSymbol}]

ClebschGordan::usage=
"ClebschGordan[{j1,m1},{j2,m2},{j3,m3}] evaluates the Clebsch-Gordan
coefficient through its relation to the 3-j symbol.  The ji must
satisfy triangularity conditions and m1+m2 = m3.";
 
ThreeJSymbol::usage=
"ThreeJSymbol[{j1,m1},{j2,m2},{j3,m3}] evaluates the 3-j symbol.  The
ji may be symbolic and must satisfy triangularity conditions. Also note
that m1+m2+m3 = 0.";
 
SixJSymbol::usage=
"The SixJSymbol[{a,b,e},{d,e,f}] returns the value of the 6-j symbol.";
 
Begin["`Private`"];
 
(*  Triangularity conditions *)

Triangular[s1_,s2_,s3_]:=
Abs[s1-s2] <= s3 <= s1 + s2 &&
   Abs[s1-s3] <= s2 <= s1 + s3 &&
   Abs[s2-s3] <= s1 <= s2 + s3
 
(*  Tolerance for rationalization *)

tol = 10^-5;
 
(*  Physical conditions *)

Physical[ji_,mi_]:= 
Which[
	NumberQ[mi] && NumberQ[ji], 
		ji >= Abs[mi] && IntegerQ[ji-mi], 
	NumberQ[ji-mi], 
		IntegerQ[ji-mi] && ji-mi >= 0,
	NumberQ[ji+mi], 
		IntegerQ[ji+mi] && ji+mi >= 0, 
	True, 
		True
]

(* Factorial and Gamma reductions *)
 
reduce =
{Factorial[a_] :> Gamma[a+1],
   Gamma[k_Integer+a_] :> Gamma[a] Pochhammer[a,k],
   Gamma[k_Rational+a_] :> Gamma[a+1/2] Pochhammer[a+1/2,k-1/2]}
 
(* ClebschGordan definition *)
 
ClebschGordan::tri= "`1` is not triangular."
ClebschGordan::phy= "`1` is not physical."
 
ClebschGordan[{j1_,m1_},{j2_,m2_},{j3_,m3_}] :=
	(-1)^Rationalize[j1-j2+m3,tol]*
		Sqrt[Rationalize[2j3+1,tol]] *
			ThreeJSymbol[{j1,m1},{j2,m2},{j3,-m3}]
 
(* ThreeJSymbol definition *)
   
ThreeJSymbol[{j1_,m1_},{j2_,m2_},{j3_,m3_}]:=
Block[{jvec = Rationalize[{j1,j2,j3},tol], 
		mvec = Rationalize[{m1,m2,m3},tol],
   jtot,rmat,flat,prod,int,pos,
   p,q,r,a,b,c,d,e,d1,e1,sig,gam,ratio,
   tri, phy, threeJ},
   tri = Triangular @ (Sequence @@ jvec);
   phy = (Expand[Plus @@ mvec] == 0 &&
	  And @@ 
		Apply[Physical, Transpose[{jvec,mvec}], 1]);
   threeJ=HoldForm[ThreeJSymbol[{j1,m1},{j2,m2},{j3,m3}]];
   jtot = Plus @@ jvec;
   If[!tri, Message[ClebschGordan::tri, threeJ]; Return[0] ];
   If[!phy, Message[ClebschGordan::phy, threeJ]; Return[0] ];
   rmat = {jtot-2jvec,jvec-mvec,jvec+mvec} // Expand;
   flat = Flatten[rmat];
   prod = Sqrt[(Times @@ ((flat)!))/((jtot+1)!)] //
   		PowerExpand;
   int = Min[Select[flat,(# >= 0 && IntegerQ[#])&]];
   pos = Select[Position[rmat,int],(Length[#]==2)&];
   {p,q,r} = If[pos=={},{1,2,3},
   		RotateLeft[{1,2,3},pos[[1,2]]-pos[[1,1]]+1]];
   a = -rmat[[2,p]];
   b = -rmat[[3,q]];
   c = -rmat[[1,r]];
   d = rmat[[3,r]]-rmat[[2,p]]+1;
   e = rmat[[2,r]]-rmat[[3,q]]+1;
   sig = rmat[[3,p]]-rmat[[2,q]];
   gam = Times @@ (Gamma /@ {1-a,1-b,1-c,d,e});
   ratio = (prod/gam //. reduce);
   (-1)^sig ratio Together[
   		HypergeometricPFQ[{a,b,c},{d1,e1},1] /.
   			{d1 -> d, e1 -> e}
   		] // PowerExpand // Cancel
   ]
 
(* SixJSymbol definition *)
 
delta[l_List]:=
Sqrt[Product[(Expand[(Plus @@ l)-2l[[i]] ])!,{i,Length[l]}]/
	(((Plus @@ l)+1)!)] // PowerExpand
 
SixJSymbol[{j1_,j2_,j3_},{j4_,j5_,j6_}]:=
Block[{jtop = Rationalize[{j1,j2,j3},tol], 
   jbot = Rationalize[{j4,j5,j6},tol],
   list,alpha,beta,rmat,flat,prod,int,pos,
   a,b,c,d,e,f,g,e1,f1,g1,k,l,n,tri,p,q,r,
   sig,gam,mult, sixJ=HoldForm[
   		SixJSymbol[{j1,j2,j3},{j4,j5,j6}]]},
   list = Rationalize[
   	{{j1,j2,j3},{j4,j5,j3},{j1,j5,j6},{j4,j2,j6}},tol];
   alpha = Apply[Plus,list,1];
   beta = Rationalize[
   	{j1+j2+j4+j5,j1+j3+j4+j6,j2+j3+j5+j6},tol];
   n = PowerExpand[Times @@ (delta /@ list)];
   tri  =  And @@ Apply[Triangular, list, 1];
   If[!tri,
   		Message[ClebschGordan::tri, sixJ]; Return[0]
   ];
   rmat = Table[beta[[k]]-alpha[[l]],
   		{l,Length[alpha]},{k,Length[beta]}] // Expand;
   flat = Flatten[rmat];
   int = Min[Select[flat,(# >= 0 && IntegerQ[#])&]];
   pos = Select[Position[rmat,int],(Length[#]==2)&];
   {p,q,r} = If[pos=={},{1,2,3},
   		RotateLeft[{1,2,3},pos[[1,2]]-1]];
   a = -rmat[[1,p]];
   b = -rmat[[2,p]];
   c = -rmat[[3,p]];
   d = -rmat[[4,p]];
   e = -rmat[[1,p]]-rmat[[2,p]]-rmat[[3,q]]-rmat[[4,r]]-1;
   f = rmat[[3,q]]-rmat[[3,p]]+1;
   g = rmat[[4,r]]-rmat[[4,p]]+1;
   sig = e+1;
   gam = Times @@ (Gamma /@ {1-a,1-b,1-c,1-d,f,g});
   mult = (-1)^sig n Gamma[1-e]/gam //. reduce;
   Together[
   		HypergeometricPFQ[ {a,b,c,d}, {e1,f1,g1},1 ]  /.
   			{e1 -> e, f1 -> f, g1 -> g}] mult //
   			PowerExpand // Cancel
]
 
End[]

Protect[ ClebschGordan, ThreeJSymbol, SixJSymbol]
 
End[]
