Someone has asked on the Maxima list some time ago whether it
was possible to compute 3j coefficients (used for addition
of 2 angular momenta) using Maxima.
I am sending a Maxima program that calculates these coefficients.
It is partly inspired by a Maple program available from the Web
page of B. G. Wybourne njsym.maple.
I have tested this program against the 3j coefficients in Abramowitz
and Stegun and those from particle data group. 6j coefficients have
been tested by comparison with a table than can be found at:
http://www.strw.leidenuniv.nl/~mathar/progs/
9j coefficients are untested.
-------------- next part --------------
/* This program computes the 3j, 6j, 9j symbols of Wigner, associated
with the irreducible representations of the group SU(2).
These functions are used in quantum mechanics for addition of angular
momenta.
References are Landau and Lifshitz t. 3: Quantum Mechanics (Pergamon)
Messiah Quantum Mechanics (Dover)
Edmonds Angular momentum in quantum Mechanics (Princeton University Press)
GPL Licence 2007
*/
wigner_3j(j1,j2,m1,m2,j,m):=block([u,v,w,k,k1,k2],
/* Check that j1,j2,j are admissible */
if ((j1<0) or (j2<0) or (j<0)) then return(0),
if (not(integerp(2*j1)) or not(integerp(2*j2)) or not(integerp(2*j))) then return(0),
if (not(integerp(j1+j2+j))) then return(0),
if ((j>(j1+j2)) or (j<abs(j1-j2))) then return(0),
if (not(integerp(2*m1)) or not(integerp(2*m2)) or (not(integerp(2*m)))) then return(0),
if ((abs(m1)>j1) or (abs(m2)>j2) or (abs(m)>j)) then return(0),
if (m#(m1+m2)) then return(0),
u:sqrt((j1+j2-j)!*(j+j1-j2)!*(j+j2-j1)!*(2*j+1)/(j+j1+j2+1)!),
/* A&S Eq. (27.9.1) p. 1006 */
k1:max(j2-j-m1,j1+m2-j,0),
k2:min(j1+j2-j,j1-m1,j2+m2),
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),
w:sqrt((j1+m1)!*(j1-m1)!*(j2+m2)!*(j2-m2)!*(j+m)!*(j-m)!),
return(rootscontract(radcan(u*v*w)))
)$
racah_delta(a,b,c):=block([],
(a+b-c)!*(b+c-a)!*(c+a-b)!/(a+b+c+1)!
)$
wigner_6j(j1,j2,j3,j4,j5,j6):=block([u,t1t2,t,v],
/* Check that j1,j2,j3,j4,j5,j6 are admissible */
if ((j1<0) or (j2<0) or (j3<0)or (j4<0) or (j5<0) or (j6<0)) then return(0),
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),
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),
/* Triangle inequalities */
if (abs(j1-j2)>j3 or j3>j1+j2) then return(0),
if (abs(j5-j4)>j3 or j3>j4+j5) then return(0),
if (abs(j5-j6)>j1 or j1>j5+j6) then return(0),
if (abs(j4-j6)>j2 or j2>j4+j6) then return(0),
/* L. D. Landau E.M. Lifschitz Mecanique Quantique (Mir) p. 513 Eq. (108,10)
A. Messiah Mecanique Quantique t. 2 (Dunod) p. 915, Eq. (36) */
u:sqrt(racah_delta(j1,j2,j3)*racah_delta(j1,j5,j6)*racah_delta(j4,j2,j6)*racah_delta(j4,j5,j3)),
t1:max(j1+j2+j3,j1+j5+j6,j4+j2+j6,j4+j5+j3),
t2:min(j1+j2+j4+j5,j2+j3+j5+j6,j1+j3+j4+j6),
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),
rootscontract(radcan(u*v))
)$
/* Adapted from njsym.maple by B. G. Wybourne */
wigner_9j(a,b,c,d,e,f,h,i,j):=block([x,xmin,xmax,result],
/* Angular momentum is half-integer */
if (not(integerp(2*a)) or not(integerp(2*b)) or not(integerp(2*c)) or not(integerp(2*d))) then return(0),
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),
/* Triangle inequalities */
if (abs(a-d)>h or h>a+d) then return(0),
if (abs(i-j)>h or h>i+j) then return(0),
if (abs(b-e)>i or i>b+e) then return(0),
if (abd(d-e)>f or f>d+e) then return(0),
if (abs(c-f)>j or j>c+f) then return(0),
if (abs(c-a)>b or b>c+a) then return(0),
xmax:min(a+j,i+d,b+f),
xmin:max(abs(a-j),abs(i-d),abs(b-f)),
/* A. Messiah Mecanique Quantique t. 2 (Dunod) p. 917, Eq. (41) */
result:sum(((-1)^(2*x))*(2*x + 1)*wigner_6j(a,d,h,i,j,x)*wigner_6j(b,e,i,d,x,f)
*wigner_6j(c,f,j,x,a,b),x,xmin,xmax),
return(rootscontract(radcan(result)))
)$