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, integer or half integer - prec - precission of the sqrt() calculation OUTPUT: double - Value of the Delta coefficient """ if(int(aa+bb-cc)!=(aa+bb-cc)): raise ValueError("j values must be integer or half integer and fullfil the triangle relation") #return 0 if(int(aa+cc-bb)!=(aa+cc-bb)): raise ValueError("j values must be integer or half integer and fullfil the triangle relation") #return 0 if(int(bb+cc-aa)!=(bb+cc-aa)): raise ValueError("j values must be integer or half integer and fullfil the triangle relation") #return 0 if(aa+bb-cc)<0: return 0 if(aa+cc-bb)<0: return 0 if(bb+cc-aa)<0: return 0 maxfact=max(aa+bb-cc,aa+cc-bb,bb+cc-aa,aa+bb+cc+1) _calc_factlist(maxfact) argsqrt=Integer(_Factlist[int(aa+bb-cc)]*_Factlist[int(aa+cc-bb)]\ *_Factlist[int(bb+cc-aa)])\ /Integer(_Factlist[int(aa+bb+cc+1)]) ressqrt=sqrt(argsqrt,prec) if type(ressqrt)==sage.rings.complex_number.ComplexNumber: res=ressqrt.real() else: res=ressqrt return res def test_bigDeltacoeff(aa,bb,cc,prec=None): r""" Test for 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, integer or half integer - prec - precission of the sqrt() calculation OUTPUT: double - Value of the Delta coefficient EXAMPLES: Simple examples: sage: test_bigDeltacoeff(1,1,1) 1/(2*sqrt(6)) """ return _bigDeltacoeff(aa,bb,cc,prec) def Racah(aa,bb,cc,dd,ee,ff,prec=None): r""" Calculate the Racah symbol W(a,b,c,d;e,f) NOTES: The Racah symbol is related to the Wigner 6j symbol: \begin{eqnarray} \left({j_1 \atop j_4} {j_2 \atop j_5} {j_3 \atop j_6} \right) = (-1)^{j_1+j_2+j_4+j_5} W(j_1,j_2,j_5,j_4,j_3,j_6) \end{eqnarray} Please see the 6j symbol for its much richer symmetries and for additional properties. ALGORITHM: This function uses the algorithm of [Edmonds74] to calculate the value of the 6j 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: a - integer or half integer b - integer or half integer c - integer or half integer d - integer or half integer e - integer or half integer f - 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 and testcases: sage: Racah(3,3,3,3,3,3) -1/14 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 """ prefac=_bigDeltacoeff(aa,bb,ee,prec)*_bigDeltacoeff(cc,dd,ee,prec) \ *_bigDeltacoeff(aa,cc,ff,prec)*_bigDeltacoeff (bb,dd,ff,prec) if prefac==0: return 0 imin=max(aa+bb+ee,cc+dd+ee,aa+cc+ff,bb+dd+ff) imax=min(aa+bb+cc+dd,aa+dd+ee+ff,bb+cc+ee+ff) maxfact=max(imax+1,aa+bb+cc+dd,aa+dd+ee+ff,bb+cc+ee+ff) _calc_factlist(maxfact) sumres=0 for kk in range(imin,imax+1): den=_Factlist[int(kk-aa-bb-ee)]*_Factlist[int(kk-cc-dd-ee)]\ *_Factlist[int(kk-aa-cc-ff)]\ *_Factlist[int(kk-bb-dd-ff)]*_Factlist[int(aa+bb+cc+dd-kk)] \ *_Factlist[int(aa+dd+ee+ff-kk)]\ *_Factlist[int(bb+cc+ee+ff-kk)] sumres=sumres+Integer((-1)**kk*_Factlist[kk+1])/den res=prefac*sumres*(-1)**(int(aa+bb+cc+dd)) return res def Wigner6j(j_1,j_2,j_3,j_4,j_5,j_6,prec=None): r""" Calculate the Wigner 6j symbol \left({j_1 \atop j_4} {j_2 \atop j_5} {j_3 \atop j_6} \right) NOTES: The Wigner 6j symbol is related to the Racah symbol but exhibits more symmetries as detailed below. \begin{eqnarray} \left({j_1 \atop j_4} {j_2 \atop j_5} {j_3 \atop j_6} \right) = (-1)^{j_1+j_2+j_4+j_5} W(j_1,j_2,j_5,j_4,j_3,j_6) \end{eqnarray} The Wigner 6j symbol obeys the following symmetry rules: - Wigner $6j$ symbols are left invariant under any permutation of the columns: \begin{eqnarray} \left({j_1 \atop j_4} {j_2 \atop j_5} {j_3 \atop j_6} \right) &=& \left({j_3 \atop j_6} {j_1 \atop j_4} {j_2 \atop j_5} \right) =\left({j_2 \atop j_5} {j_3 \atop j_6} {j_2 \atop j_4} \right) \qquad \mbox{cyclic} \\ &=& \left({j_3 \atop j_6} {j_2 \atop j_5} {j_1 \atop j_4} \right) =\left({j_1 \atop j_4} {j_3 \atop j_6} {j_2 \atop j_5} \right) =\left({j_2 \atop j_5} {j_1 \atop j_4} {j_3 \atop j_6} \right) \qquad \mbox{anti-cyclic} \end{eqnarray} - They are invariant under the exchange of the upper and lower arguments in each of any two columns, i. e. \begin{eqnarray} \left({j_1 \atop j_4} {j_2 \atop j_5} {j_3 \atop j_6} \right) = \left({j_1 \atop j_4} {j_5 \atop j_2} {j_6 \atop j_3} \right) = \left({j_4 \atop j_1} {j_2 \atop j_5} {j_6 \atop j_3} \right) = \left({j_4 \atop j_1} {j_5 \atop j_2} {j_3 \atop j_6} \right) \end{eqnarray} - additional 6 symmetries [Regge59] giving rise to 144 symmetries in total - only non-zero if any triple of $j$'s fulfil a triangle relation ALGORITHM: This function uses the algorithm of [Edmonds74] to calculate the value of the 6j 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 j_4 - integer or half integer j_5 - integer or half integer j_6 - 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 and testcases: sage: Wigner6j(3,3,3,3,3,3) -1/14 sage: Wigner6j(5,5,5,5,5,5) 1/52 sage: Wigner6j(6,6,6,6,6,6) 309/10868 sage: Wigner6j(8,8,8,8,8,8) -12219/965770 sage: Wigner6j(30,30,30,30,30,30) 36082186869033479581/87954851694828981714124 sage: Wigner6j(0.5,0.5,1,0.5,0.5,1) 1/6 sage: Wigner6j(200,200,200,200,200,200, prec=1000)*1.0 0.000155903212413242 It is an error to have arguments that are not integer or half integer values or do not fulfill the triangle relation: sage: Wigner6j(2.5,2.5,2.5,2.5,2.5,2.5) Traceback (most recent call last): ... ValueError: j values must be integer or half integer and fullfil the triangle relation sage: Wigner6j(0.5,0.5,1.1,0.5,0.5,1.1) Traceback (most recent call last): ... ValueError: j values must be integer or half integer and fullfil the triangle relation REFERENCES: - [Regge59] 'Symmetry Properties of Racah Coefficients', T. Regge, Nuovo Cimento, Volume 11, pp. 116 (1959) - [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) """ res=(-1)**(int(j_1+j_2+j_4+j_5))*Racah (j_1,j_2,j_5,j_4,j_3,j_6,prec) return res def Wigner9j(j_1,j_2,j_3,j_4,j_5,j_6,j_7,j_8,j_9,prec=None): r""" Calculate the Wigner 9j symbol \left\{ \begin{matrix} {j_1} & {j_2} & {j_3} \\ {j_4} & {j_5} & {j_6} \\ {j_7} & {j_8} & {j_9} \end{matrix} \right\} 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 j_4 - integer or half integer j_5 - integer or half integer j_6 - integer or half integer j_7 - integer or half integer j_8 - integer or half integer j_9 - 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 and testcases, note that for speed reasons a precission is given: sage: Wigner9j(1,1,1, 1,1,1, 1,1,0 ,prec=64) # ==1/18 0.0555555555555555555 sage: Wigner9j(1,1,1, 1,1,1, 1,1,1) 0 sage: Wigner9j(1,1,1, 1,1,1, 1,1,2 ,prec=64) # ==1/18 0.0555555555555555556 sage: Wigner9j(1,2,1, 2,2,2, 1,2,1 ,prec=64) # ==-1/150 -0.00666666666666666667 sage: Wigner9j(3,3,2, 2,2,2, 3,3,2 ,prec=64) # ==157/14700 0.0106802721088435374 sage: Wigner9j(3,3,2, 3,3,2, 3,3,2 ,prec=64) # ==3221*sqrt(70)/ (246960*sqrt(105)) - 365/(3528*sqrt(70)*sqrt(105)) 0.00944247746651111739 sage: Wigner9j(3,3,1, 3.5,3.5,2, 3.5,3.5,1 ,prec=64) # ==3221*sqrt(70)/(246960*sqrt(105)) - 365/(3528*sqrt(70)*sqrt(105)) 0.0110216678544351364 sage: Wigner9j(100,80,50, 50,100,70, 60,50,100 ,prec=1000)*1.0 1.05597798065761e-7 sage: Wigner9j(30,30,10, 30.5,30.5,20, 30.5,30.5,10 ,prec=1000) *1.0 # ==(80944680186359968990/95103769817469)*sqrt (1/682288158959699477295) 0.0000325841699408828 sage: Wigner9j(64,62.5,114.5, 61.5,61,112.5, 113.5,110.5,60, prec=1000)*1.0 -3.41407910055520e-39 sage: Wigner9j(15,15,15, 15,3,15, 15,18,10, prec=1000)*1.0 -0.0000778324615309539 sage: Wigner9j(1.5,1,1.5, 1,1,1, 1.5,1,1.5) 0 It is an error to have arguments that are not integer or half integer values or do not fulfill the triangle relation: sage: Wigner9j(0.5,0.5,0.5, 0.5,0.5,0.5, 0.5,0.5,0.5,prec=64) Traceback (most recent call last): ... ValueError: j values must be integer or half integer and fullfil the triangle relation sage: Wigner9j(1,1,1, 0.5,1,1.5, 0.5,1,2.5,prec=64) Traceback (most recent call last): ... ValueError: j values must be integer or half integer and fullfil the triangle relation 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) """ imin=0 imax=min(j_1+j_9,j_2+j_6,j_4+j_8) sumres=0 for kk in range(imin,imax+1): sumres=sumres+(2*kk+1)*Racah(j_1,j_2,j_9,j_6,j_3,kk,prec)\ *Racah(j_4,j_6,j_8,j_2,j_5,kk,prec)\ *Racah(j_1,j_4,j_9,j_8,j_7,kk,prec) return sumres def Gaunt(l_1,l_2,l_3,m_1,m_2,m_3,prec=None): r""" Calculate the Gaunt coefficient which is defined as the integral over three spherical harmonics: \begin{eqnarray} Y{j_1 \atop m_1} {j_2 \atop m_2} {j_3 \atop m_3} &=& \int Y_{l_1,m_1}(\Omega) Y_{l_2,m_2}(\Omega) Y_{l_3,m_3}(\Omega)\; d\Omega \\ &=& \sqrt{\frac{(2l_1+1)(2l_2+1)(2l_3+1)}{4\pi}} \left({l_1 \atop 0 } {l_2 \atop 0 } {l_3 \atop 0} \right) \left({l_1 \atop m_1} {l_2 \atop m_2} {l_3 \atop m_3} \right) \end{eqnarray} NOTES: The Gaunt coefficient obeys the following symmetry rules: - invariant under any permutation of the columns \begin{eqnarray} Y{j_1 \atop m_1} {j_2 \atop m_2} {j_3 \atop m_3} &=& Y{j_3 \atop m_3} {j_1 \atop m_1} {j_2 \atop m_2} =Y{j_2 \atop m_2} {j_3 \atop m_3} {j_1 \atop m_1} \qquad \mbox{cyclic} \\ &=& Y{j_3 \atop m_3} {j_2 \atop m_2} {j_1 \atop m_1} =Y{j_1 \atop m_1} {j_3 \atop m_3} {j_2 \atop m_2} =Y{j_2 \atop m_2} {j_1 \atop m_1} {j_3 \atop m_3} \qquad\mbox{anti-cyclic} \end{eqnarray} - invariant under space inflection, i.e. \begin{eqnarray} Y{j_1 \atop m_1} {j_2 \atop m_2} {j_3 \atop m_3} = Y{j_1 \atop -m_1} {j_2 \atop -m_2} {j_3 \atop -m_3} \end{eqnarray} - symmetric with respect to the 72 Regge symmetries as inherited for the $3j$ symbols [Regge58] - zero for $l_1$, $l_2$, $l_3$ not fulfilling triangle relation - zero for violating any one of the conditions: $l_1\ge|m_1|$, $l_2\ge|m_2|$, $l_3\ge|m_3| - non-zero only for an even sum of the $l_i$, i. e. $J=l_1+l_2+l_3=2n$ for $n$ in $\mathbf{N}$ ALGORITHM: This function uses the algorithm of [Liberatodebrito82] to calculate the value of the Gaunt coefficient 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 j_2 - integer j_3 - integer m_1 - integer m_2 - integer m_3 - 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: sage: Gaunt(1,0,1,1,0,-1) -1/(2*sqrt(pi)) sage: Gaunt(1,0,1,1,0,0) 0 sage: Gaunt(29,29,34,10,-5,-5) 1821867940156*sqrt(22134)/(215552371055153321*sqrt(pi)) sage: Gaunt(20,20,40,1,-1,0) 28384503878959800/(74029560764440771*sqrt(pi)) sage: Gaunt(12,15,5,2,3,-5) 91*sqrt(36890)/(124062*sqrt(pi)) sage: Gaunt(10,10,12,9,3,-12) -98*sqrt(6279)/(62031*sqrt(pi)) sage: Gaunt(1000,1000,1200,9,3,-12).n(64) 0.00689500421922113448 It is an error to use non-integer values for $l$ and $m$: sage: Gaunt(1.2,0,1.2,0,0,0) Traceback (most recent call last): ... ValueError: l values must be integer sage: Gaunt(1,0,1,1.1,0,-1.1) Traceback (most recent call last): ... ValueError: m values must be integer REFERENCES: - [Regge58] 'Symmetry Properties of Clebsch-Gordan Coefficients', T. Regge, Nuovo Cimento, Volume 10, pp. 544 (1958) - [Liberatodebrito82] 'FORTRAN program for the integral of three spherical harmonics', A. Liberato de Brito, Comput. Phys. Commun., Volume 25, pp. 81-85 (1982) - [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 """ if int(l_1)!=l_1 or int(l_2)!=l_2 or int(l_3)!=l_3: raise ValueError("l values must be integer") #return 0 if int(m_1)!=m_1 or int(m_2)!=m_2 or int(m_3)!=m_3: raise ValueError("m values must be integer") bigL=(l_1+l_2+l_3)/2 a1=l_1+l_2-l_3 if (a1<0): return 0 a2=l_1-l_2+l_3 if (a2<0): return 0 a3=-l_1+l_2+l_3 if (a3<0): return 0 if Mod(2*bigL,2)<>0: # if int(int(2*bigL)/2)<>int(bigL): return 0 if (m_1+m_2+m_3)<>0: return 0 if (abs(m_1)>l_1) or (abs(m_2)>l_2) or (abs(m_3)>l_3): return 0 imin=max(-l_3+l_1+m_2,-l_3+l_2-m_1,0) imax=min(l_2+m_2,l_1-m_1,l_1+l_2-l_3) maxfact=max(l_1+l_2+l_3+1,imax+1) _calc_factlist(maxfact) argsqrt=(2*l_1+1)*(2*l_2+1)*(2*l_3+1)*_Factlist[l_1-m_1]\ *_Factlist[l_1+m_1]*_Factlist[l_2-m_2]*_Factlist [l_2+m_2]\ *_Factlist[l_3-m_3]*_Factlist[l_3+m_3]/(4*pi) ressqrt=sqrt(argsqrt,prec) if type(ressqrt)==sage.rings.complex_number.ComplexNumber: ressqrt=ressqrt.real() prefac=Integer(_Factlist[bigL]*_Factlist[l_2-l_1+l_3]\ *_Factlist[l_1-l_2+l_3]\ *_Factlist[l_1+l_2-l_3])/_Factlist[2*bigL+1]\ /(_Factlist[bigL-l_1]*_Factlist[bigL-l_2]*_Factlist[bigL-l_3]) sumres=0 for ii in range(imin,imax+1): den=_Factlist[ii]*_Factlist[ii+l_3-l_1-m_2]\ *_Factlist[l_2+m_2-ii]*_Factlist[l_1-ii-m_1]\ *_Factlist[ii+l_3-l_2+m_1]*_Factlist[l_1+l_2-l_3-ii] sumres=sumres+Integer((-1)**ii)/den res=ressqrt*prefac*sumres*(-1)**(bigL+l_3+m_1-m_2) return res --~--~---------~--~----~------------~-------~--~----~ To post to this group, send email to sage-devel@googlegroups.com To unsubscribe from this group, send email to sage-devel-unsubscr...@googlegroups.com For more options, visit this group at http://groups.google.com/group/sage-devel URLs: http://www.sagemath.org -~----------~----~----~----~------~----~------~--~---
