Ok, here it is:
--------------------------------------------------------------------------- --------------------
r"""
Calculate Wigner 3j, 6j, 9j, Clebsch-Gordan, Racah and Gaunt
coefficients
Collection of functions for calculating er 3j, 6j, 9j, Clebsch-Gordan,
Racah as well as Gaunt coefficients exactly, all evaluating to a
rational number times the square root of a rational number [Rasch03].
Please see the description of the individual functions for further
details and examples.
REFERENCES:
- [Rasch03] 'Efficient Storage Scheme for Precalculated Wigner 3j, 6j
and Gaunt Coefficients', J. Rasch and A. C. H. Yu, SIAM
J. Sci. Comput. Volume 25, Issue 4, pp. 1416-1428 (2003)
AUTHORS:
- Jens Rasch (2009-03-24): initial version for Sage
"""
#***********************************************************************
# Copyright (C) 2008 Jens Rasch <jyr2...@gmail.com>
#
# Distributed under the terms of the GNU General Public License (GPL)
# http://www.gnu.org/licenses/
#***********************************************************************
from sage.all import *
# from sage.calculus.calculus import sqrt,factorial
# from sage.functions.constants import pi
# This list of precomputed factorials is needed to massively
# accelerate consequetive calculations of the various coefficients
_Factlist=[1]
def _calc_factlist(nn):
if nn>=len(_Factlist):
for ii in range(len(_Factlist),nn+1):
_Factlist.append(_Factlist[ii-1]*ii)
#_Factlist.append(factorial(ii))
#return _Factlist
def test_calc_factlist(nn):
r"""
Function calculates a list of precomputed factorials in order to
massively accelerate consequetive calculations of the various
coefficients.
INPUT:
- nn Highest factorial to be computed
OUTPUT:
integer list of factorials
EXAMPLES:
Calculate list of factorials:
sage: test_calc_factlist(10)
[1, 1, 2, 6, 24, 120, 720, 5040, 40320, 362880, 3628800]
"""
_calc_factlist(nn)
return _Factlist
def Wigner3j(j_1,j_2,j_3,m_1,m_2,m_3,prec=None):
r"""
Calculate the Wigner 3j symbol
\left({j_1 \atop m_1} {j_2 \atop m_2} {j_3 \atop m_3} \right)
NOTES:
The Wigner 3j symbol obeys the following symmetry rules:
- invariant under any permutation of the columns (with the
exception of a sign change where $J:=j_1+j_2+j_3$):
\begin{eqnarray}
\left({j_1 \atop m_1} {j_2 \atop m_2} {j_3 \atop m_3}\right)
&=&
\left({j_3 \atop m_3} {j_1 \atop m_1} {j_2 \atop m_2}\right)
=\left({j_2 \atop m_2} {j_3 \atop m_3} {j_1 \atop m_1}\right)
\qquad \mbox{cyclic}
&=&
(-)^{J}\left({j_3\atop m_3} {j_2\atop m_2} {j_1\atop
m_1}\right)
=(-)^{J}\left({j_1\atop m_1} {j_3\atop m_3} {j_2\atop
m_2}\right)
=(-)^{J}\left({j_2\atop m_2} {j_1\atop m_1} {j_3\atop
m_3}\right)
\qquad\mbox{anti-cyclic}
\end{eqnarray}
- invariant under space inflection, i. e.
\begin{eqnarray}
\left({{j_1}\atop{m_1}} {{j_2}\atop{m_2}} {{j_3}\atop{m_3}}
\right)
=
(-)^{J}
\left({j_1 \atop -m_1} {j_2 \atop -m_2}{j_3 \atop -m_3}\right)
\end{eqnarray}
- symmetric with respect to the 72 additional symmetries based on
the work by [Regge58]
- zero for $j_1$, $j_2$, $j_3$ not fulfilling triangle relation
- zero for ${m_1}+{m_2}+{m_3}!= 0$
- zero for violating any one of the conditions
$j_1\ge|m_1|$, $j_2\ge|m_2|$, $j_3\ge|m_3|$
ALGORITHM:
This function uses the algorithm of [Edmonds74] to calculate the
value of the 3j symbol exactly. Note that the formula contains
alternating sums over large factorials and is therefore unsuitable
for finite precision arithmetic and only useful for a computer
algebra system [Rasch03].
INPUT:
j_1 - integer or half integer
j_2 - integer or half integer
j_3 - integer or half integer
m_1 - integer or half integer
m_2 - integer or half integer
m_3 - integer or half integer
prec - precission, default: None. Providing a precission can
drastically spead up the calculation.
OUTPUT:
The result will be a rational number times the square root of a
rational number, unless a precission is given.
EXAMPLES:
A couple of examples:
sage: Wigner3j(2,6,4,0,0,0)
sqrt(5)/sqrt(143)
sage: Wigner3j(2,6,4,0,0,1)
0
sage: Wigner3j(0.5,0.5,1,0.5,-0.5,0)
1/sqrt(6)
sage: Wigner3j(40,100,60,-10,60,-50)
95608*sqrt(21082735836735314343364163310)/(18702538494885*sqrt
(220491455010479533763))
sage: Wigner3j(2500,2500,5000,2488,2400,-4888 ,prec=64)
7.60424456883448589e-12
It is an error to have arguments that are not integer or half
integer values:
sage: Wigner3j(2.1,6,4,0,0,0)
Traceback (most recent call last):
...
ValueError: j values must be integer or half integer
sage: Wigner3j(2,6,4,1,0,-1.1)
Traceback (most recent call last):
...
ValueError: m values must be integer or half integer
REFERENCES:
- [Regge58] 'Symmetry Properties of Clebsch-Gordan Coefficients',
T. Regge, Nuovo Cimento, Volume 10, pp. 544 (1958)
- [Edmonds74] 'Angular Momentum in Quantum Mechanics',
A. R. Edmonds, Princeton University Press (1974)
- [Rasch03] 'Efficient Storage Scheme for Precalculated Wigner 3j,
6j and Gaunt Coefficients', J. Rasch and A. C. H. Yu, SIAM
J. Sci. Comput. Volume 25, Issue 4, pp. 1416-1428 (2003)
AUTHORS:
- Jens Rasch (2009-03-24): initial version
"""
if int(j_1*2)!=j_1*2 or int(j_2*2)!=j_2*2 or int(j_3*2)!=j_3*2:
raise ValueError("j values must be integer or half integer")
#return 0
if int(m_1*2)!=m_1*2 or int(m_2*2)!=m_2*2 or int(m_3*2)!=m_3*2:
raise ValueError("m values must be integer or half integer")
#return 0
if (m_1+m_2+m_3<>0):
return 0
prefid=Integer((-1)**(int(j_1-j_2-m_3)))
m_3=-m_3
a1=j_1+j_2-j_3
if (a1<0):
return 0
a2=j_1-j_2+j_3
if (a2<0):
return 0
a3=-j_1+j_2+j_3
if (a3<0):
return 0
if (abs(m_1)>j_1) or (abs(m_2)>j_2) or (abs(m_3)>j_3):
return 0
maxfact=max(j_1+j_2+j_3+1,j_1+abs(m_1),j_2+abs(m_2),j_3+abs(m_3))
_calc_factlist(maxfact)
#try:
argsqrt=Integer(_Factlist[int(j_1+j_2-j_3)]\
*_Factlist[int(j_1-j_2+j_3)]\
*_Factlist[int(-j_1+j_2+j_3)]\
*_Factlist[int(j_1-m_1)]*_Factlist[int(j_1+m_1)]\
*_Factlist[int(j_2-m_2)]*_Factlist[int(j_2+m_2)]\
*_Factlist[int(j_3-m_3)]*_Factlist[j_3+m_3])\
/_Factlist[int(j_1+j_2+j_3+1)]\
#except:
# return 0
ressqrt=sqrt(argsqrt,prec)
if type(ressqrt)==sage.rings.complex_number.ComplexNumber:
ressqrt=ressqrt.real()
imin=max(-j_3+j_1+m_2,-j_3+j_2-m_1,0)
imax=min(j_2+m_2,j_1-m_1,j_1+j_2-j_3)
sumres=0
for ii in range(imin,imax+1):
den=_Factlist[ii]*_Factlist[int(ii+j_3-j_1-m_2)]\
*_Factlist[int(j_2+m_2-ii)]*_Factlist[int(j_1-ii-m_1)]\
*_Factlist[int(ii+j_3-j_2+m_1)]\
*_Factlist[int(j_1+j_2-j_3-ii)]
sumres=sumres+Integer((-1)**ii)/den
res=ressqrt*sumres*prefid
return res
def ClebschGordan(j_1,j_2,j_3,m_1,m_2,m_3,prec=None):
r"""
Calculates the Clebsch-Gordan coefficient
$\left< j_1 m_1 \; j_2 m_2 | j_3 m_3 \right>$
NOTES:
The Clebsch-Gordan coefficient will be evaluated via its relation
to Wigner 3j symbols:
\begin{eqnarray}
\left< j_1 m_1 \; j_2 m_2 | j_3 m_3 \right>=
(-1)^(j_1-j_2+m_3) \sqrt(2j_3+1)
\left({j_1 \atop m_1} {j_2 \atop m_2} {j_3 \atop -m_3}\right)
\end{eqnarray}
See also the documentation on Wigner 3j symbols which exhibit much
higher symmetry releations that the Clebsch-Gordan coefficient.
INPUT:
j_1 - integer or half integer
j_2 - integer or half integer
j_3 - integer or half integer
m_1 - integer or half integer
m_2 - integer or half integer
m_3 - integer or half integer
prec - precission, default: None. Providing a precission can
drastically spead up the calculation.
OUTPUT:
The result will be a rational number times the square root of a
rational number, unless a precission is given.
EXAMPLES:
A couple of examples:
sage: ClebschGordan(3/2,1/2,2, 3/2,1/2,2)
1
sage: ClebschGordan(1.5,0.5,1, 1.5,-0.5,1)
sqrt(3)/2
sage: ClebschGordan(3/2,1/2,1, -1/2,1/2,0)
-sqrt(3)/sqrt(6)
REFERENCES:
- [Edmonds74] 'Angular Momentum in Quantum Mechanics',
A. R. Edmonds, Princeton University Press (1974)
- [Rasch03] 'Efficient Storage Scheme for Precalculated Wigner 3j,
6j and Gaunt Coefficients', J. Rasch and A. C. H. Yu, SIAM
J. Sci. Comput. Volume 25, Issue 4, pp. 1416-1428 (2003)
AUTHORS:
- Jens Rasch (2009-03-24): initial version
"""
res=(-1)**(int(j_1-j_2+m_3))*sqrt(2*j_3+1)\
*Wigner3j(j_1,j_2,j_3,m_1,m_2,-m_3,prec)
return res
def _bigDeltacoeff(aa,bb,cc,prec=None):
r"""
Calculates the Delta coefficient of the 3 angular momenta for
Racah symbols. Also checks that the differences are of integer
value.
INPUT:
- aa - first angular momentum, integer or half integer
- bb - second angular momentum, integer or half integer
- cc - third angular momentum,
...
read more »
