r"""
Wigner, Clebsch-Gordan, Racah, and Gaunt coefficients
Collection of functions for calculating Wigner 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] J. Rasch and A. C. H. Yu, 'Efficient Storage Scheme for
Pre-calculated Wigner 3j, 6j and Gaunt Coefficients', SIAM
J. Sci. Comput. Volume 25, Issue 4, pp. 1416-1428 (2003)
Credits and Copyright
~~~~~~~~~~~~~~~~~~~~~
This code was taken from Sage with the permission of all authors:
https://groups.google.com/forum/#!topic/sage-devel/M4NZdu-7O38
AUTHORS:
- Jens Rasch (2009-03-24): initial version for Sage
- Jens Rasch (2009-05-31): updated to sage-4.0
Copyright (C) 2008 Jens Rasch <jyr2000@gmail.com>
"""
from __future__ import print_function, division
from sympy import (Integer, pi, sqrt, sympify, Dummy, S, Sum, Ynm,
Function)
from sympy.core.compatibility import range
# This list of precomputed factorials is needed to massively
# accelerate future calculations of the various coefficients
_Factlist = [1]
def _calc_factlist(nn):
r"""
Function calculates a list of precomputed factorials in order to
massively accelerate future calculations of the various
coefficients.
INPUT:
- ``nn`` - integer, highest factorial to be computed
OUTPUT:
list of integers -- the list of precomputed factorials
EXAMPLES:
Calculate list of factorials::
sage: from sage.functions.wigner import _calc_factlist
sage: _calc_factlist(10)
[1, 1, 2, 6, 24, 120, 720, 5040, 40320, 362880, 3628800]
"""
if nn >= len(_Factlist):
for ii in range(len(_Factlist), int(nn + 1)):
_Factlist.append(_Factlist[ii - 1] * ii)
return _Factlist[:int(nn) + 1]
[docs]def wigner_3j(j_1, j_2, j_3, m_1, m_2, m_3):
r"""
Calculate the Wigner 3j symbol `\operatorname{Wigner3j}(j_1,j_2,j_3,m_1,m_2,m_3)`.
INPUT:
- ``j_1``, ``j_2``, ``j_3``, ``m_1``, ``m_2``, ``m_3`` - integer or half integer
OUTPUT:
Rational number times the square root of a rational number.
Examples
========
>>> from sympy.physics.wigner import wigner_3j
>>> wigner_3j(2, 6, 4, 0, 0, 0)
sqrt(715)/143
>>> wigner_3j(2, 6, 4, 0, 0, 1)
0
It is an error to have arguments that are not integer or half
integer values::
sage: wigner_3j(2.1, 6, 4, 0, 0, 0)
Traceback (most recent call last):
...
ValueError: j values must be integer or half integer
sage: wigner_3j(2, 6, 4, 1, 0, -1.1)
Traceback (most recent call last):
...
ValueError: m values must be integer or half integer
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`):
.. math::
\begin{aligned}
\operatorname{Wigner3j}(j_1,j_2,j_3,m_1,m_2,m_3)
&=\operatorname{Wigner3j}(j_3,j_1,j_2,m_3,m_1,m_2) \\
&=\operatorname{Wigner3j}(j_2,j_3,j_1,m_2,m_3,m_1) \\
&=(-1)^J \operatorname{Wigner3j}(j_3,j_2,j_1,m_3,m_2,m_1) \\
&=(-1)^J \operatorname{Wigner3j}(j_1,j_3,j_2,m_1,m_3,m_2) \\
&=(-1)^J \operatorname{Wigner3j}(j_2,j_1,j_3,m_2,m_1,m_3)
\end{aligned}
- invariant under space inflection, i.e.
.. math::
\operatorname{Wigner3j}(j_1,j_2,j_3,m_1,m_2,m_3)
=(-1)^J \operatorname{Wigner3j}(j_1,j_2,j_3,-m_1,-m_2,-m_3)
- 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 \neq 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]_.
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)
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")
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")
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(int(maxfact))
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[int(j_3 + m_3)]) / \
_Factlist[int(j_1 + j_2 + j_3 + 1)]
ressqrt = sqrt(argsqrt)
if ressqrt.is_complex:
ressqrt = ressqrt.as_real_imag()[0]
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(int(imin), int(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
[docs]def clebsch_gordan(j_1, j_2, j_3, m_1, m_2, m_3):
r"""
Calculates the Clebsch-Gordan coefficient
`\left\langle j_1 m_1 \; j_2 m_2 | j_3 m_3 \right\rangle`.
The reference for this function is [Edmonds74]_.
INPUT:
- ``j_1``, ``j_2``, ``j_3``, ``m_1``, ``m_2``, ``m_3`` - integer or half integer
OUTPUT:
Rational number times the square root of a rational number.
EXAMPLES::
>>> from sympy import S
>>> from sympy.physics.wigner import clebsch_gordan
>>> clebsch_gordan(S(3)/2, S(1)/2, 2, S(3)/2, S(1)/2, 2)
1
>>> clebsch_gordan(S(3)/2, S(1)/2, 1, S(3)/2, -S(1)/2, 1)
sqrt(3)/2
>>> clebsch_gordan(S(3)/2, S(1)/2, 1, -S(1)/2, S(1)/2, 0)
-sqrt(2)/2
NOTES:
The Clebsch-Gordan coefficient will be evaluated via its relation
to Wigner 3j symbols:
.. math::
\left\langle j_1 m_1 \; j_2 m_2 | j_3 m_3 \right\rangle
=(-1)^{j_1-j_2+m_3} \sqrt{2j_3+1}
\operatorname{Wigner3j}(j_1,j_2,j_3,m_1,m_2,-m_3)
See also the documentation on Wigner 3j symbols which exhibit much
higher symmetry relations than the Clebsch-Gordan coefficient.
AUTHORS:
- Jens Rasch (2009-03-24): initial version
"""
res = (-1) ** sympify(j_1 - j_2 + m_3) * sqrt(2 * j_3 + 1) * \
wigner_3j(j_1, j_2, j_3, m_1, m_2, -m_3)
return res
def _big_delta_coeff(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`` - precision of the ``sqrt()`` calculation
OUTPUT:
double - Value of the Delta coefficient
EXAMPLES::
sage: from sage.functions.wigner import _big_delta_coeff
sage: _big_delta_coeff(1,1,1)
1/2*sqrt(1/6)
"""
if int(aa + bb - cc) != (aa + bb - cc):
raise ValueError("j values must be integer or half integer and fulfill the triangle relation")
if int(aa + cc - bb) != (aa + cc - bb):
raise ValueError("j values must be integer or half integer and fulfill the triangle relation")
if int(bb + cc - aa) != (bb + cc - aa):
raise ValueError("j values must be integer or half integer and fulfill the triangle relation")
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)
if prec:
ressqrt = ressqrt.evalf(prec).as_real_imag()[0]
return ressqrt
[docs]def racah(aa, bb, cc, dd, ee, ff, prec=None):
r"""
Calculate the Racah symbol `W(a,b,c,d;e,f)`.
INPUT:
- ``a``, ..., ``f`` - integer or half integer
- ``prec`` - precision, default: ``None``. Providing a precision can
drastically speed up the calculation.
OUTPUT:
Rational number times the square root of a rational number
(if ``prec=None``), or real number if a precision is given.
Examples
========
>>> from sympy.physics.wigner import racah
>>> racah(3,3,3,3,3,3)
-1/14
NOTES:
The Racah symbol is related to the Wigner 6j symbol:
.. math::
\operatorname{Wigner6j}(j_1,j_2,j_3,j_4,j_5,j_6)
=(-1)^{j_1+j_2+j_4+j_5} W(j_1,j_2,j_5,j_4,j_3,j_6)
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]_.
AUTHORS:
- Jens Rasch (2009-03-24): initial version
"""
prefac = _big_delta_coeff(aa, bb, ee, prec) * \
_big_delta_coeff(cc, dd, ee, prec) * \
_big_delta_coeff(aa, cc, ff, prec) * \
_big_delta_coeff(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(int(imin), int(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
[docs]def wigner_6j(j_1, j_2, j_3, j_4, j_5, j_6, prec=None):
r"""
Calculate the Wigner 6j symbol `\operatorname{Wigner6j}(j_1,j_2,j_3,j_4,j_5,j_6)`.
INPUT:
- ``j_1``, ..., ``j_6`` - integer or half integer
- ``prec`` - precision, default: ``None``. Providing a precision can
drastically speed up the calculation.
OUTPUT:
Rational number times the square root of a rational number
(if ``prec=None``), or real number if a precision is given.
Examples
========
>>> from sympy.physics.wigner import wigner_6j
>>> wigner_6j(3,3,3,3,3,3)
-1/14
>>> wigner_6j(5,5,5,5,5,5)
1/52
It is an error to have arguments that are not integer or half
integer values or do not fulfill the triangle relation::
sage: wigner_6j(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 fulfill the triangle relation
sage: wigner_6j(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 fulfill the triangle relation
NOTES:
The Wigner 6j symbol is related to the Racah symbol but exhibits
more symmetries as detailed below.
.. math::
\operatorname{Wigner6j}(j_1,j_2,j_3,j_4,j_5,j_6)
=(-1)^{j_1+j_2+j_4+j_5} W(j_1,j_2,j_5,j_4,j_3,j_6)
The Wigner 6j symbol obeys the following symmetry rules:
- Wigner 6j symbols are left invariant under any permutation of
the columns:
.. math::
\begin{aligned}
\operatorname{Wigner6j}(j_1,j_2,j_3,j_4,j_5,j_6)
&=\operatorname{Wigner6j}(j_3,j_1,j_2,j_6,j_4,j_5) \\
&=\operatorname{Wigner6j}(j_2,j_3,j_1,j_5,j_6,j_4) \\
&=\operatorname{Wigner6j}(j_3,j_2,j_1,j_6,j_5,j_4) \\
&=\operatorname{Wigner6j}(j_1,j_3,j_2,j_4,j_6,j_5) \\
&=\operatorname{Wigner6j}(j_2,j_1,j_3,j_5,j_4,j_6)
\end{aligned}
- They are invariant under the exchange of the upper and lower
arguments in each of any two columns, i.e.
.. math::
\operatorname{Wigner6j}(j_1,j_2,j_3,j_4,j_5,j_6)
=\operatorname{Wigner6j}(j_1,j_5,j_6,j_4,j_2,j_3)
=\operatorname{Wigner6j}(j_4,j_2,j_6,j_1,j_5,j_3)
=\operatorname{Wigner6j}(j_4,j_5,j_3,j_1,j_2,j_6)
- additional 6 symmetries [Regge59]_ giving rise to 144 symmetries
in total
- only non-zero if any triple of `j`'s fulfill 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]_.
REFERENCES:
.. [Regge59] 'Symmetry Properties of Racah Coefficients',
T. Regge, Nuovo Cimento, Volume 11, pp. 116 (1959)
"""
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
[docs]def wigner_9j(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
`\operatorname{Wigner9j}(j_1,j_2,j_3,j_4,j_5,j_6,j_7,j_8,j_9)`.
INPUT:
- ``j_1``, ..., ``j_9`` - integer or half integer
- ``prec`` - precision, default: ``None``. Providing a precision can
drastically speed up the calculation.
OUTPUT:
Rational number times the square root of a rational number
(if ``prec=None``), or real number if a precision is given.
Examples
========
>>> from sympy.physics.wigner import wigner_9j
>>> wigner_9j(1,1,1, 1,1,1, 1,1,0 ,prec=64) # ==1/18
0.05555555...
>>> wigner_9j(1/2,1/2,0, 1/2,3/2,1, 0,1,1 ,prec=64) # ==1/6
0.1666666...
It is an error to have arguments that are not integer or half
integer values or do not fulfill the triangle relation::
sage: wigner_9j(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 fulfill the triangle relation
sage: wigner_9j(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 fulfill the triangle relation
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]_.
"""
imax = int(min(j_1 + j_9, j_2 + j_6, j_4 + j_8) * 2)
imin = imax % 2
sumres = 0
for kk in range(imin, int(imax) + 1, 2):
sumres = sumres + (kk + 1) * \
racah(j_1, j_2, j_9, j_6, j_3, kk / 2, prec) * \
racah(j_4, j_6, j_8, j_2, j_5, kk / 2, prec) * \
racah(j_1, j_4, j_9, j_8, j_7, kk / 2, prec)
return sumres
[docs]def gaunt(l_1, l_2, l_3, m_1, m_2, m_3, prec=None):
r"""
Calculate the Gaunt coefficient.
The Gaunt coefficient is defined as the integral over three
spherical harmonics:
.. math::
\begin{aligned}
\operatorname{Gaunt}(l_1,l_2,l_3,m_1,m_2,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}}
\operatorname{Wigner3j}(l_1,l_2,l_3,0,0,0)
\operatorname{Wigner3j}(l_1,l_2,l_3,m_1,m_2,m_3)
\end{aligned}
INPUT:
- ``l_1``, ``l_2``, ``l_3``, ``m_1``, ``m_2``, ``m_3`` - integer
- ``prec`` - precision, default: ``None``. Providing a precision can
drastically speed up the calculation.
OUTPUT:
Rational number times the square root of a rational number
(if ``prec=None``), or real number if a precision is given.
Examples
========
>>> from sympy.physics.wigner import gaunt
>>> gaunt(1,0,1,1,0,-1)
-1/(2*sqrt(pi))
>>> 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
NOTES:
The Gaunt coefficient obeys the following symmetry rules:
- invariant under any permutation of the columns
.. math::
\begin{aligned}
Y(l_1,l_2,l_3,m_1,m_2,m_3)
&=Y(l_3,l_1,l_2,m_3,m_1,m_2) \\
&=Y(l_2,l_3,l_1,m_2,m_3,m_1) \\
&=Y(l_3,l_2,l_1,m_3,m_2,m_1) \\
&=Y(l_1,l_3,l_2,m_1,m_3,m_2) \\
&=Y(l_2,l_1,l_3,m_2,m_1,m_3)
\end{aligned}
- invariant under space inflection, i.e.
.. math::
Y(l_1,l_2,l_3,m_1,m_2,m_3)
=Y(l_1,l_2,l_3,-m_1,-m_2,-m_3)
- 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.
`L = l_1 + l_2 + l_3 = 2n` for `n` in `\mathbb{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]_.
REFERENCES:
.. [Liberatodebrito82] 'FORTRAN program for the integral of three
spherical harmonics', A. Liberato de Brito,
Comput. Phys. Commun., Volume 25, pp. 81-85 (1982)
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")
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")
sumL = l_1 + l_2 + l_3
bigL = sumL // 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 sumL % 2:
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)
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(int(imin), int(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 * Integer((-1) ** (bigL + l_3 + m_1 - m_2))
if prec is not None:
res = res.n(prec)
return res
class Wigner3j(Function):
def doit(self, **hints):
if all(obj.is_number for obj in self.args):
return wigner_3j(*self.args)
else:
return self
[docs]def dot_rot_grad_Ynm(j, p, l, m, theta, phi):
r"""
Returns dot product of rotational gradients of spherical harmonics.
This function returns the right hand side of the following expression:
.. math ::
\vec{R}Y{_j^{p}} \cdot \vec{R}Y{_l^{m}} = (-1)^{m+p}
\sum\limits_{k=|l-j|}^{l+j}Y{_k^{m+p}} * \alpha_{l,m,j,p,k} *
\frac{1}{2} (k^2-j^2-l^2+k-j-l)
Arguments
=========
j, p, l, m .... indices in spherical harmonics (expressions or integers)
theta, phi .... angle arguments in spherical harmonics
Example
=======
>>> from sympy import symbols
>>> from sympy.physics.wigner import dot_rot_grad_Ynm
>>> theta, phi = symbols("theta phi")
>>> dot_rot_grad_Ynm(3, 2, 2, 0, theta, phi).doit()
3*sqrt(55)*Ynm(5, 2, theta, phi)/(11*sqrt(pi))
"""
j = sympify(j)
p = sympify(p)
l = sympify(l)
m = sympify(m)
theta = sympify(theta)
phi = sympify(phi)
k = Dummy("k")
def alpha(l,m,j,p,k):
return sqrt((2*l+1)*(2*j+1)*(2*k+1)/(4*pi)) * \
Wigner3j(j, l, k, S(0), S(0), S(0)) * Wigner3j(j, l, k, p, m, -m-p)
return (-S(1))**(m+p) * Sum(Ynm(k, m+p, theta, phi) * alpha(l,m,j,p,k) / 2 \
*(k**2-j**2-l**2+k-j-l), (k, abs(l-j), l+j))
