<s> r"""<STR_LIT>""" <EOL> from __future__ import print_function , division <EOL> from sympy import ( Integer , pi , sqrt , sympify , Dummy , S , Sum , Ynm , <EOL> Function ) <EOL> from sympy . core . compatibility import range <EOL> _Factlist = [ <NUM_LIT:1> ] <EOL> def _calc_factlist ( nn ) : <EOL> r"""<STR_LIT>""" <EOL> if nn >= len ( _Factlist ) : <EOL> for ii in range ( len ( _Factlist ) , int ( nn + <NUM_LIT:1> ) ) : <EOL> _Factlist . append ( _Factlist [ ii - <NUM_LIT:1> ] * ii ) <EOL> return _Factlist [ : int ( nn ) + <NUM_LIT:1> ] <EOL> def wigner_3j ( j_1 , j_2 , j_3 , m_1 , m_2 , m_3 ) : <EOL> r"""<STR_LIT>""" <EOL> if int ( j_1 * <NUM_LIT:2> ) != j_1 * <NUM_LIT:2> or int ( j_2 * <NUM_LIT:2> ) != j_2 * <NUM_LIT:2> or int ( j_3 * <NUM_LIT:2> ) != j_3 * <NUM_LIT:2> : <EOL> raise ValueError ( "<STR_LIT>" ) <EOL> if int ( m_1 * <NUM_LIT:2> ) != m_1 * <NUM_LIT:2> or int ( m_2 * <NUM_LIT:2> ) != m_2 * <NUM_LIT:2> or int ( m_3 * <NUM_LIT:2> ) != m_3 * <NUM_LIT:2> : <EOL> raise ValueError ( "<STR_LIT>" ) <EOL> if m_1 + m_2 + m_3 != <NUM_LIT:0> : <EOL> return <NUM_LIT:0> <EOL> prefid = Integer ( ( - <NUM_LIT:1> ) ** int ( j_1 - j_2 - m_3 ) ) <EOL> m_3 = - m_3 <EOL> a1 = j_1 + j_2 - j_3 <EOL> if a1 < <NUM_LIT:0> : <EOL> return <NUM_LIT:0> <EOL> a2 = j_1 - j_2 + j_3 <EOL> if a2 < <NUM_LIT:0> : <EOL> return <NUM_LIT:0> <EOL> a3 = - j_1 + j_2 + j_3 <EOL> if a3 < <NUM_LIT:0> : <EOL> return <NUM_LIT:0> <EOL> if ( abs ( m_1 ) > j_1 ) or ( abs ( m_2 ) > j_2 ) or ( abs ( m_3 ) > j_3 ) : <EOL> return <NUM_LIT:0> <EOL> maxfact = max ( j_1 + j_2 + j_3 + <NUM_LIT:1> , j_1 + abs ( m_1 ) , j_2 + abs ( m_2 ) , <EOL> j_3 + abs ( m_3 ) ) <EOL> _calc_factlist ( int ( maxfact ) ) <EOL> argsqrt = Integer ( _Factlist [ int ( j_1 + j_2 - j_3 ) ] * <EOL> _Factlist [ int ( j_1 - j_2 + j_3 ) ] * <EOL> _Factlist [ int ( - j_1 + j_2 + j_3 ) ] * <EOL> _Factlist [ int ( j_1 - m_1 ) ] * <EOL> _Factlist [ int ( j_1 + m_1 ) ] * <EOL> _Factlist [ int ( j_2 - m_2 ) ] * <EOL> _Factlist [ int ( j_2 + m_2 ) ] * <EOL> _Factlist [ int ( j_3 - m_3 ) ] * <EOL> _Factlist [ int ( j_3 + m_3 ) ] ) / _Factlist [ int ( j_1 + j_2 + j_3 + <NUM_LIT:1> ) ] <EOL> ressqrt = sqrt ( argsqrt ) <EOL> if ressqrt . is_complex : <EOL> ressqrt = ressqrt . as_real_imag ( ) [ <NUM_LIT:0> ] <EOL> imin = max ( - j_3 + j_1 + m_2 , - j_3 + j_2 - m_1 , <NUM_LIT:0> ) <EOL> imax = min ( j_2 + m_2 , j_1 - m_1 , j_1 + j_2 - j_3 ) <EOL> sumres = <NUM_LIT:0> <EOL> for ii in range ( int ( imin ) , int ( imax ) + <NUM_LIT:1> ) : <EOL> 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 ) ] <EOL> sumres = sumres + Integer ( ( - <NUM_LIT:1> ) ** ii ) / den <EOL> res = ressqrt * sumres * prefid <EOL> return res <EOL> def clebsch_gordan ( j_1 , j_2 , j_3 , m_1 , m_2 , m_3 ) : <EOL> r"""<STR_LIT>""" <EOL> res = ( - <NUM_LIT:1> ) ** sympify ( j_1 - j_2 + m_3 ) * sqrt ( <NUM_LIT:2> * j_3 + <NUM_LIT:1> ) * wigner_3j ( j_1 , j_2 , j_3 , m_1 , m_2 , - m_3 ) <EOL> return res <EOL> def _big_delta_coeff ( aa , bb , cc , prec = None ) : <EOL> r"""<STR_LIT>""" <EOL> if int ( aa + bb - cc ) != ( aa + bb - cc ) : <EOL> raise ValueError ( "<STR_LIT>" ) <EOL> if int ( aa + cc - bb ) != ( aa + cc - bb ) : <EOL> raise ValueError ( "<STR_LIT>" ) <EOL> if int ( bb + cc - aa ) != ( bb + cc - aa ) : <EOL> raise ValueError ( "<STR_LIT>" ) <EOL> if ( aa + bb - cc ) < <NUM_LIT:0> : <EOL> return <NUM_LIT:0> <EOL> if ( aa + cc - bb ) < <NUM_LIT:0> : <EOL> return <NUM_LIT:0> <EOL> if ( bb + cc - aa ) < <NUM_LIT:0> : <EOL> return <NUM_LIT:0> <EOL> maxfact = max ( aa + bb - cc , aa + cc - bb , bb + cc - aa , aa + bb + cc + <NUM_LIT:1> ) <EOL> _calc_factlist ( maxfact ) <EOL> argsqrt = Integer ( _Factlist [ int ( aa + bb - cc ) ] * <EOL> _Factlist [ int ( aa + cc - bb ) ] * <EOL> _Factlist [ int ( bb + cc - aa ) ] ) / Integer ( _Factlist [ int ( aa + bb + cc + <NUM_LIT:1> ) ] ) <EOL> ressqrt = sqrt ( argsqrt ) <EOL> if prec : <EOL> ressqrt = ressqrt . evalf ( prec ) . as_real_imag ( ) [ <NUM_LIT:0> ] <EOL> return ressqrt <EOL> def racah ( aa , bb , cc , dd , ee , ff , prec = None ) : <EOL> r"""<STR_LIT>""" <EOL> 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 ) <EOL> if prefac == <NUM_LIT:0> : <EOL> return <NUM_LIT:0> <EOL> imin = max ( aa + bb + ee , cc + dd + ee , aa + cc + ff , bb + dd + ff ) <EOL> imax = min ( aa + bb + cc + dd , aa + dd + ee + ff , bb + cc + ee + ff ) <EOL> maxfact = max ( imax + <NUM_LIT:1> , aa + bb + cc + dd , aa + dd + ee + ff , <EOL> bb + cc + ee + ff ) <EOL> _calc_factlist ( maxfact ) <EOL> sumres = <NUM_LIT:0> <EOL> for kk in range ( int ( imin ) , int ( imax ) + <NUM_LIT:1> ) : <EOL> 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 ) ] <EOL> sumres = sumres + Integer ( ( - <NUM_LIT:1> ) ** kk * _Factlist [ kk + <NUM_LIT:1> ] ) / den <EOL> res = prefac * sumres * ( - <NUM_LIT:1> ) ** int ( aa + bb + cc + dd ) <EOL> return res <EOL> def wigner_6j ( j_1 , j_2 , j_3 , j_4 , j_5 , j_6 , prec = None ) : <EOL> r"""<STR_LIT>""" <EOL> res = ( - <NUM_LIT:1> ) ** int ( j_1 + j_2 + j_4 + j_5 ) * racah ( j_1 , j_2 , j_5 , j_4 , j_3 , j_6 , prec ) <EOL> return res <EOL> def wigner_9j ( j_1 , j_2 , j_3 , j_4 , j_5 , j_6 , j_7 , j_8 , j_9 , prec = None ) : <EOL> r"""<STR_LIT>""" <EOL> imin = <NUM_LIT:0> <EOL> imax = min ( j_1 + j_9 , j_2 + j_6 , j_4 + j_8 ) <EOL> sumres = <NUM_LIT:0> <EOL> for kk in range ( imin , int ( imax ) + <NUM_LIT:1> ) : <EOL> sumres = sumres + ( <NUM_LIT:2> * kk + <NUM_LIT: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 ) <EOL> return sumres <EOL> def gaunt ( l_1 , l_2 , l_3 , m_1 , m_2 , m_3 , prec = None ) : <EOL> r"""<STR_LIT>""" <EOL> if int ( l_1 ) != l_1 or int ( l_2 ) != l_2 or int ( l_3 ) != l_3 : <EOL> raise ValueError ( "<STR_LIT>" ) <EOL> if int ( m_1 ) != m_1 or int ( m_2 ) != m_2 or int ( m_3 ) != m_3 : <EOL> raise ValueError ( "<STR_LIT>" ) <EOL> sumL = l_1 + l_2 + l_3 <EOL> bigL = sumL // <NUM_LIT:2> <EOL> a1 = l_1 + l_2 - l_3 <EOL> if a1 < <NUM_LIT:0> : <EOL> return <NUM_LIT:0> <EOL> a2 = l_1 - l_2 + l_3 <EOL> if a2 < <NUM_LIT:0> : <EOL> return <NUM_LIT:0> <EOL> a3 = - l_1 + l_2 + l_3 <EOL> if a3 < <NUM_LIT:0> : <EOL> return <NUM_LIT:0> <EOL> if sumL % <NUM_LIT:2> : <EOL> return <NUM_LIT:0> <EOL> if ( m_1 + m_2 + m_3 ) != <NUM_LIT:0> : <EOL> return <NUM_LIT:0> <EOL> if ( abs ( m_1 ) > l_1 ) or ( abs ( m_2 ) > l_2 ) or ( abs ( m_3 ) > l_3 ) : <EOL> return <NUM_LIT:0> <EOL> imin = max ( - l_3 + l_1 + m_2 , - l_3 + l_2 - m_1 , <NUM_LIT:0> ) <EOL> imax = min ( l_2 + m_2 , l_1 - m_1 , l_1 + l_2 - l_3 ) <EOL> maxfact = max ( l_1 + l_2 + l_3 + <NUM_LIT:1> , imax + <NUM_LIT:1> ) <EOL> _calc_factlist ( maxfact ) <EOL> argsqrt = ( <NUM_LIT:2> * l_1 + <NUM_LIT:1> ) * ( <NUM_LIT:2> * l_2 + <NUM_LIT:1> ) * ( <NUM_LIT:2> * l_3 + <NUM_LIT: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 ] / ( <NUM_LIT:4> * pi ) <EOL> ressqrt = sqrt ( argsqrt ) <EOL> prefac = Integer ( _Factlist [ bigL ] * _Factlist [ l_2 - l_1 + l_3 ] * <EOL> _Factlist [ l_1 - l_2 + l_3 ] * _Factlist [ l_1 + l_2 - l_3 ] ) / _Factlist [ <NUM_LIT:2> * bigL + <NUM_LIT:1> ] / ( _Factlist [ bigL - l_1 ] * <EOL> _Factlist [ bigL - l_2 ] * _Factlist [ bigL - l_3 ] ) <EOL> sumres = <NUM_LIT:0> <EOL> for ii in range ( int ( imin ) , int ( imax ) + <NUM_LIT:1> ) : <EOL> 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 ] <EOL> sumres = sumres + Integer ( ( - <NUM_LIT:1> ) ** ii ) / den <EOL> res = ressqrt * prefac * sumres * ( - <NUM_LIT:1> ) ** ( bigL + l_3 + m_1 - m_2 ) <EOL> if prec is not None : <EOL> res = res . n ( prec ) <EOL> return res <EOL> class Wigner3j ( Function ) : <EOL> def doit ( self , ** hints ) : <EOL> if all ( obj . is_number for obj in self . args ) : <EOL> return wigner_3j ( * self . args ) <EOL> else : <EOL> return self <EOL> def dot_rot_grad_Ynm ( j , p , l , m , theta , phi ) : <EOL> r"""<STR_LIT>""" <EOL> j = sympify ( j ) <EOL> p = sympify ( p ) <EOL> l = sympify ( l ) <EOL> m = sympify ( m ) <EOL> theta = sympify ( theta ) <EOL> phi = sympify ( phi ) <EOL> k = Dummy ( "<STR_LIT:k>" ) <EOL> def alpha ( l , m , j , p , k ) : <EOL> return sqrt ( ( <NUM_LIT:2> * l + <NUM_LIT:1> ) * ( <NUM_LIT:2> * j + <NUM_LIT:1> ) * ( <NUM_LIT:2> * k + <NUM_LIT:1> ) / ( <NUM_LIT:4> * pi ) ) * Wigner3j ( j , l , k , S ( <NUM_LIT:0> ) , S ( <NUM_LIT:0> ) , S ( <NUM_LIT:0> ) ) * Wigner3j ( j , l , k , p , m , - m - p ) <EOL> return ( - S ( <NUM_LIT:1> ) ) ** ( m + p ) * Sum ( Ynm ( k , m + p , theta , phi ) * alpha ( l , m , j , p , k ) / <NUM_LIT:2> * ( k ** <NUM_LIT:2> - j ** <NUM_LIT:2> - l ** <NUM_LIT:2> + k - j - l ) , ( k , abs ( l - j ) , l + j ) ) </s>
