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