share/contrib/clebsch_gordan.mac

   1 /* This program computes the 3j, 6j, 9j symbols of Wigner, associated
   2 with the irreducible representations of the group SU(2).
   3 These functions are used in quantum mechanics for addition of angular
   4 momenta.
   5 References are Landau and Lifshitz t. 3: Quantum Mechanics (Pergamon)
   6                Messiah Quantum Mechanics (Dover)
   7                Edmonds Angular momentum in quantum Mechanics (Princeton University Press)
   8 Page and equation numbers refer to the french edition of Landau and Lifshitz "Mecanique Quantique" (Mir) and Messiah "Mecanique Quantique" (Dunod).
   9 GPL (v2)  Licence 2007,2011
  10 */
  11
  12
  13 /* computes the Clebsch-Gordan coefficient <j1,j2,m1,m2|j,m> */
  14 clebsch_gordan(j1,j2,m1,m2,j,m):=block([u,v,w,k,k1,k2],
  15 /* Check that j1,j2,j are admissible */
  16 if ((j1<0) or (j2<0) or (j<0)) then return(0),
  17 if (not(integerp(2*j1)) or not(integerp(2*j2)) or not(integerp(2*j))) then return(0),
  18 if (not(integerp(j1+j2+j))) then return(0),
  19 if ((j>(j1+j2)) or (j<abs(j1-j2))) then return(0),
  20 if (not(integerp(2*m1)) or not(integerp(2*m2)) or (not(integerp(2*m)))) then return(0),
  21 if ((abs(m1)>j1) or (abs(m2)>j2) or (abs(m)>j)) then return(0),
  22 if (m#(m1+m2)) then return(0),
  23 u:sqrt((j1+j2-j)!*(j+j1-j2)!*(j+j2-j1)!*(2*j+1)/(j+j1+j2+1)!),
  24 /* A&S Eq. (27.9.1) p. 1006. Note that A&S misleadingly uses the terminology Wigner coefficients */
  25 k1:max(j2-j-m1,j1+m2-j,0),
  26 k2:min(j1+j2-j,j1-m1,j2+m2),
  27 v:sum((-1)^k/(k!*(j1+j2-j-k)!*(j1-m1-k)!*(j2+m2-k)!*(j-j2+m1+k)!*(j-j1-m2+k)!),k,k1,k2),
  28 w:sqrt((j1+m1)!*(j1-m1)!*(j2+m2)!*(j2-m2)!*(j+m)!*(j-m)!),
  29 return(rootscontract(radcan(u*v*w)))
  30 )$
  31
  32 /* Computes Racah's triangle function */
  33 racah_delta(a,b,c):=block([],
  34 (a+b-c)!*(b+c-a)!*(c+a-b)!/(a+b+c+1)!
  35 )$
  36
  37 /* Computes Wigner's 3j symbol. Messiah (op.cit. p.908, Eq.12) */
  38
  39 wigner_3j(j1,j2,j3,m1,m2,m3):=block([],
  40 (-1)^(j1-j2-m3)/sqrt(2*j3+1)*clebsch_gordan(j1,j2,m1,m2,j3,-m3)
  41 )$
  42
  43   /* Computes Racah's V coefficient. Messiah (op. cit. p. 908, footnote 4)*/
  44 racah_v(a,b,c,a1,b1,c1):=block([],
  45 (-1)^(c-c1)/sqrt(2*c+1)*clebsch_gordan(a,b,a1,b1,c,-c1)
  46 )$
  47
  48
  49 /* Computes Wigner's 6j symbol */
  50 wigner_6j(j1,j2,j3,j4,j5,j6):=block([u,t1t2,t,v],
  51 /* Check that j1,j2,j3,j4,j5,j6 are admissible */
  52 if ((j1<0) or (j2<0) or (j3<0)or (j4<0) or (j5<0) or (j6<0)) then return(0),
  53 if (not(integerp(2*j1)) or not(integerp(2*j2)) or not(integerp(2*j3)) or not(integerp(2*j4)) or not(integerp(2*j5)) or not(integerp(2*j6))) then return (0),
  54 if (not(integerp(j1+j2+j3)) or not(integerp(j4+j5+j3)) or not(integerp(j4+j2+j6)) or (not(integerp(j1+j5+j6)))) then return(0),
  55 /* Triangle inequalities */
  56 if (abs(j1-j2)>j3 or j3>j1+j2) then return(0),
  57 if (abs(j5-j4)>j3 or j3>j4+j5) then return(0),
  58 if (abs(j5-j6)>j1 or j1>j5+j6) then return(0),
  59 if (abs(j4-j6)>j2 or j2>j4+j6) then return(0),
  60 /* L. D. Landau E.M. Lifschitz Mecanique Quantique (Mir) p. 513 Eq. (108,10)
  61 with J1=j4,J2,j5,J3=j6
  62    A. Messiah Mecanique Quantique t. 2 (Dunod) p. 915, Eq. (36) */
  63 u:sqrt(racah_delta(j1,j2,j3)*racah_delta(j1,j5,j6)*racah_delta(j4,j2,j6)*racah_delta(j4,j5,j3)),
  64 t1:max(j1+j2+j3,j1+j5+j6,j4+j2+j6,j4+j5+j3),
  65 t2:min(j1+j2+j4+j5,j2+j3+j5+j6,j1+j3+j4+j6),
  66 v:sum((-1)^t*(t+1)!/((t-j1-j2-j3)!*(t-j1-j5-j6)!*(t-j4-j2-j6)!*(t-j4-j5-j3)!*(j1+j2+j4+j5-t)!*(j2+j3+j5+j6-t)!*(j1+j3+j4+j6-t)!),t,t1,t2),
  67 rootscontract(radcan(u*v))
  68 )$
  69
  70 /* Racah's W coefficient, Messiah (op. cit.) p. 913, Eq. (30) */
  71 racah_w(j1,j2,j5,j4,j3,j6):=block([],
  72 (-1)^(j1+j2+j4+j5)*wigner_6j(j1,j2,j3,j4,j5,j6)
  73 )$
  74
  75 /* Computes Wigner's 9j symbol. Adapted from njsym.maple by B. G. Wybourne */
  76 wigner_9j(a,b,c,d,e,f,h,i,j):=block([x,xmin,xmax,result],
  77 /* Angular momentum is half-integer */
  78  if (not(integerp(2*a)) or not(integerp(2*b)) or not(integerp(2*c)) or not(integerp(2*d))) then return(0),
  79 if (not(integerp(2*e)) or not(integerp(2*f)) or not(integerp(2*h)) or not(integerp(2*i)) or not(integerp(2*j))) then return(0),
  80 /* Triangle inequalities */
  81  if (abs(a-d)>h or h>a+d) then return(0),
  82  if (abs(i-j)>h or h>i+j) then return(0),
  83  if (abs(b-e)>i or i>b+e) then return(0),
  84  if (abd(d-e)>f or f>d+e) then return(0),
  85  if (abs(c-f)>j or j>c+f) then return(0),
  86  if (abs(c-a)>b or b>c+a) then return(0),
  87  xmax:min(a+j,i+d,b+f),
  88  xmin:max(abs(a-j),abs(i-d),abs(b-f)),
  89 /* A. Messiah  Mecanique Quantique t. 2 (Dunod) p. 917, Eq. (41) */
  90  result:sum(((-1)^(2*x))*(2*x + 1)*wigner_6j(a,d,h,i,j,x)*wigner_6j(b,e,i,d,x,f)
  91 *wigner_6j(c,f,j,x,a,b),x,xmin,xmax),
  92 return(rootscontract(radcan(result)))
  93 )$