(* :Context: Wigner` *)

(* :Title: Wigner *)

(* :Author: K.Takeda *)

(* :Package Version: 1.1 *)

(* :History:
  original version(Version1.0) : September 1996
  Version1.1 : 
*)

(* :Reference:
   J. J. Sakurai, Moderen Quantum Mechanics (Addison Wesley)
*)


BeginPackage["Wigner`"]


WignerD::usage = 
 "WignerD[{j, m, n}, {alpha, beta, gamma}] generates 
  the element of the Wigner rotation operator.
  Both Abs[m] and Abs[n] must not exceed j.
  WignerD[{j, m, n}, beta] generates the element of the reduced Wigner rotation operator."

RotateTensor::usage = 
  "RotateTensor[a_List, k_Integer, {alpha_, beta_, gamma_}] performs spherical tensor rotation according to the Wigner formula."

Begin["`Private`"]

et[{j_, m_, n_}, beta_, k_] := et[{j,m,n},beta,k] =
  Module[
         {sig, n1, n2, n3, n4,
               d1, d2, d3, d4, test,
                   num, den, cpow, spow},
         
         sig = k - n + m;
         
         n1 = (j + n)!;
         n2 = (j - n)!;
         n3 = (j + m)!;
         n4 = (j - m)!;
         
         d1 = j + n - k;
         d2 = k;
         d3 = j - k - m;
         d4 = k - n + m;
         
         test = d1<0 || d2<0 || d3<0 || d4<0;
         
         If[test, Return[0]];
         
         d1 = d1!; d2 = d2!; d3 = d3!; d4 = d4!;
         
         num = Sqrt[n1 n2 n3 n4];
         den = d1 d2 d3 d4;
         
         cpow = 2j - 2k + n - m;
         spow = 2k - n + m;
         
         (-1)^sig num/den *
                Cos[beta/2]^cpow Sin[beta/2]^spow
        ]

WignerD[{j_, m_, n_}, {alpha_, beta_, gamma_}] :=
WignerD[{j,m,n}, {alpha,beta,gamma}] =
	 WignerD[{j, m, n}, beta] *
	   Exp[-I(m alpha + n gamma)]
		
WignerD[{j_, m_, n_}, beta_] := WignerD[{j,m,n},beta] =
	    ( Plus @@
		     Table[
		           et[{j, m, n}, beta, i],
		                 {i, 0, j + n}
		          ] ) /. Cos[_]^2 + Sin[_]^2 -> 1


RotateTensor[a_List, k_Integer, {alpha_, beta_, gamma_}] :=
Table[
	Sum[a[[q+k+1]] WignerD[{k, q, p}, {alpha, beta, gamma}],
		{q, -k, k} ],
	 {p, -k, k}
	 ] /; (Length[a] == 2k+1)


End[]

EndPackage[]