home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 9 / 09.iso / e / e032 / 3.ddi / FILES / PRELOAD.PAK / CLEBSCHG.M < prev    next >
Encoding:
Text File  |  1992-07-29  |  5.7 KB  |  196 lines
  1. (* :Title: Clebsch-Gordan Coefficients *)
  2.  
  3. (* :Author: Paul C. Abbott  *)
  4.  
  5.  
  6. (* :Summary:
  7.   This package computes the Clebsch Gordan, 3-j coefficients (Wigner),
  8.   and 6-j coefficients (Wigner) for both numeric and symbolic arguments.
  9. *)
  10.  
  11. (* :Context: ClebschGordan` *)
  12.  
  13.  
  14. (* :Package Version: 1.1 *)
  15.  
  16. (* :Copyright: Copyright 1990 Wolfram Research, Inc. *)
  17.  
  18. (* :History:
  19.   Version 1.0 by Stephen Wolfram (Wolfram Research), May 1987.
  20.   Extended to Symbolic arguments by Paul C. Abbott, February 1990.
  21.   Modified for new behaviour of Hypergeometric functions
  22.       by Paul C. Abbott, November 1991.
  23.   Modified to rationalize real input, February 1992.
  24. *)
  25.  
  26. (* :Keywords:
  27.   Clebsch-Gordan coefficients, Wigner coefficients, 3-j symbols,
  28.   6-j symbols, Racah coefficients, spin statistics, vector addition
  29.   coefficients, spherical harmonics, angular momentum
  30. *)
  31.  
  32. (* :Source:
  33.   Edmonds, A. R.: "Angular momenta in quantum mechanics",
  34.   Princeton University Press, 1974
  35. *)
  36.  
  37. (* :Source:
  38.   Rao, K. R. and Chiu, C. B.: Journal of Physics, A,
  39.   Vol. 22 pp3779-88, 1989
  40. *)
  41.  
  42. (* :Mathematica Version: 2.0 *)
  43.  
  44. Begin["System`"]
  45.  
  46. Unprotect[ ClebschGordan, ThreeJSymbol, SixJSymbol]
  47.  
  48. Map[ Clear, {ClebschGordan, ThreeJSymbol, SixJSymbol}]
  49.  
  50. ClebschGordan::usage=
  51. "ClebschGordan[{j1,m1},{j2,m2},{j3,m3}] evaluates the Clebsch-Gordan
  52. coefficient through its relation to the 3-j symbol.  The ji must
  53. satisfy triangularity conditions and m1+m2 = m3.";
  54.  
  55. ThreeJSymbol::usage=
  56. "ThreeJSymbol[{j1,m1},{j2,m2},{j3,m3}] evaluates the 3-j symbol.  The
  57. ji may be symbolic and must satisfy triangularity conditions. Also note
  58. that m1+m2+m3 = 0.";
  59.  
  60. SixJSymbol::usage=
  61. "The SixJSymbol[{a,b,e},{d,e,f}] returns the value of the 6-j symbol.";
  62.  
  63. Begin["`Private`"];
  64.  
  65. (*  Triangularity conditions *)
  66.  
  67. Triangular[s1_,s2_,s3_]:=
  68. Abs[s1-s2] <= s3 <= s1 + s2 &&
  69.    Abs[s1-s3] <= s2 <= s1 + s3 &&
  70.    Abs[s2-s3] <= s1 <= s2 + s3
  71.  
  72. (*  Tolerance for rationalization *)
  73.  
  74. tol = 10^-5;
  75.  
  76. (*  Physical conditions *)
  77.  
  78. Physical[ji_,mi_]:= 
  79. Which[
  80.     NumberQ[mi] && NumberQ[ji], 
  81.         ji >= Abs[mi] && IntegerQ[ji-mi], 
  82.     NumberQ[ji-mi], 
  83.         IntegerQ[ji-mi] && ji-mi >= 0,
  84.     NumberQ[ji+mi], 
  85.         IntegerQ[ji+mi] && ji+mi >= 0, 
  86.     True, 
  87.         True
  88. ]
  89.  
  90. (* Factorial and Gamma reductions *)
  91.  
  92. reduce =
  93. {Factorial[a_] :> Gamma[a+1],
  94.    Gamma[k_Integer+a_] :> Gamma[a] Pochhammer[a,k],
  95.    Gamma[k_Rational+a_] :> Gamma[a+1/2] Pochhammer[a+1/2,k-1/2]}
  96.  
  97. (* ClebschGordan definition *)
  98.  
  99. ClebschGordan::tri= "`1` is not triangular."
  100. ClebschGordan::phy= "`1` is not physical."
  101.  
  102. ClebschGordan[{j1_,m1_},{j2_,m2_},{j3_,m3_}] :=
  103.     (-1)^Rationalize[j1-j2+m3,tol]*
  104.         Sqrt[Rationalize[2j3+1,tol]] *
  105.             ThreeJSymbol[{j1,m1},{j2,m2},{j3,-m3}]
  106.  
  107. (* ThreeJSymbol definition *)
  108.    
  109. ThreeJSymbol[{j1_,m1_},{j2_,m2_},{j3_,m3_}]:=
  110. Block[{jvec = Rationalize[{j1,j2,j3},tol], 
  111.         mvec = Rationalize[{m1,m2,m3},tol],
  112.    jtot,rmat,flat,prod,int,pos,
  113.    p,q,r,a,b,c,d,e,d1,e1,sig,gam,ratio,
  114.    tri, phy, threeJ},
  115.    tri = Triangular @ (Sequence @@ jvec);
  116.    phy = (Expand[Plus @@ mvec] == 0 &&
  117.       And @@ 
  118.         Apply[Physical, Transpose[{jvec,mvec}], 1]);
  119.    threeJ=HoldForm[ThreeJSymbol[{j1,m1},{j2,m2},{j3,m3}]];
  120.    jtot = Plus @@ jvec;
  121.    If[!tri, Message[ClebschGordan::tri, threeJ]; Return[0] ];
  122.    If[!phy, Message[ClebschGordan::phy, threeJ]; Return[0] ];
  123.    rmat = {jtot-2jvec,jvec-mvec,jvec+mvec} // Expand;
  124.    flat = Flatten[rmat];
  125.    prod = Sqrt[(Times @@ ((flat)!))/((jtot+1)!)] //
  126.            PowerExpand;
  127.    int = Min[Select[flat,(# >= 0 && IntegerQ[#])&]];
  128.    pos = Select[Position[rmat,int],(Length[#]==2)&];
  129.    {p,q,r} = If[pos=={},{1,2,3},
  130.            RotateLeft[{1,2,3},pos[[1,2]]-pos[[1,1]]+1]];
  131.    a = -rmat[[2,p]];
  132.    b = -rmat[[3,q]];
  133.    c = -rmat[[1,r]];
  134.    d = rmat[[3,r]]-rmat[[2,p]]+1;
  135.    e = rmat[[2,r]]-rmat[[3,q]]+1;
  136.    sig = rmat[[3,p]]-rmat[[2,q]];
  137.    gam = Times @@ (Gamma /@ {1-a,1-b,1-c,d,e});
  138.    ratio = (prod/gam //. reduce);
  139.    (-1)^sig ratio Together[
  140.            HypergeometricPFQ[{a,b,c},{d1,e1},1] /.
  141.                {d1 -> d, e1 -> e}
  142.            ] // PowerExpand // Cancel
  143.    ]
  144.  
  145. (* SixJSymbol definition *)
  146.  
  147. delta[l_List]:=
  148. Sqrt[Product[(Expand[(Plus @@ l)-2l[[i]] ])!,{i,Length[l]}]/
  149.     (((Plus @@ l)+1)!)] // PowerExpand
  150.  
  151. SixJSymbol[{j1_,j2_,j3_},{j4_,j5_,j6_}]:=
  152. Block[{jtop = Rationalize[{j1,j2,j3},tol], 
  153.    jbot = Rationalize[{j4,j5,j6},tol],
  154.    list,alpha,beta,rmat,flat,prod,int,pos,
  155.    a,b,c,d,e,f,g,e1,f1,g1,k,l,n,tri,p,q,r,
  156.    sig,gam,mult, sixJ=HoldForm[
  157.            SixJSymbol[{j1,j2,j3},{j4,j5,j6}]]},
  158.    list = Rationalize[
  159.        {{j1,j2,j3},{j4,j5,j3},{j1,j5,j6},{j4,j2,j6}},tol];
  160.    alpha = Apply[Plus,list,1];
  161.    beta = Rationalize[
  162.        {j1+j2+j4+j5,j1+j3+j4+j6,j2+j3+j5+j6},tol];
  163.    n = PowerExpand[Times @@ (delta /@ list)];
  164.    tri  =  And @@ Apply[Triangular, list, 1];
  165.    If[!tri,
  166.            Message[ClebschGordan::tri, sixJ]; Return[0]
  167.    ];
  168.    rmat = Table[beta[[k]]-alpha[[l]],
  169.            {l,Length[alpha]},{k,Length[beta]}] // Expand;
  170.    flat = Flatten[rmat];
  171.    int = Min[Select[flat,(# >= 0 && IntegerQ[#])&]];
  172.    pos = Select[Position[rmat,int],(Length[#]==2)&];
  173.    {p,q,r} = If[pos=={},{1,2,3},
  174.            RotateLeft[{1,2,3},pos[[1,2]]-1]];
  175.    a = -rmat[[1,p]];
  176.    b = -rmat[[2,p]];
  177.    c = -rmat[[3,p]];
  178.    d = -rmat[[4,p]];
  179.    e = -rmat[[1,p]]-rmat[[2,p]]-rmat[[3,q]]-rmat[[4,r]]-1;
  180.    f = rmat[[3,q]]-rmat[[3,p]]+1;
  181.    g = rmat[[4,r]]-rmat[[4,p]]+1;
  182.    sig = e+1;
  183.    gam = Times @@ (Gamma /@ {1-a,1-b,1-c,1-d,f,g});
  184.    mult = (-1)^sig n Gamma[1-e]/gam //. reduce;
  185.    Together[
  186.            HypergeometricPFQ[ {a,b,c,d}, {e1,f1,g1},1 ]  /.
  187.                {e1 -> e, f1 -> f, g1 -> g}] mult //
  188.                PowerExpand // Cancel
  189. ]
  190.  
  191. End[]
  192.  
  193. Protect[ ClebschGordan, ThreeJSymbol, SixJSymbol]
  194.  
  195. End[]
  196.