https://github.com/JoeMcEwen/FAST-PT
Revision 137418e116f80d982427f13f4aed9273503952a8 authored by Jonathan Blazek on 17 October 2018, 13:18:51 UTC, committed by Jonathan Blazek on 17 October 2018, 13:18:51 UTC
1 parent 19472bf
Tip revision: 137418e116f80d982427f13f4aed9273503952a8 authored by Jonathan Blazek on 17 October 2018, 13:18:51 UTC
Added LICENSE
Added LICENSE
Tip revision: 137418e
Wigner_symbols.py
''' Wigner symbols.
J.E. McEwen 2016
Algorithm based on the sympy implimentation of sympy.physics.wigner,
which was based off of:
[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)
Additional routines have been added for special configurations of Wigner 3 j symbols
The code has not be written in a way to handle large input values; beware of precission issues
associated with large values aquired from factorials.
'''
import numpy as np
from scipy.misc import factorial
import sys
# store factorials in an array
def factorial_list(N):
x=np.arange(N+1)
return factorial(x)
# old way for the factorials, no longer used.
# def factorial_list(N):
#
# FL=[1]
# # N = highest factorial needed for computation
# # note_to_self : maybe it is faster to use the scipy factorial
# # function and then append to Factorial_list
#
# #if N > len(FL):
# if N > 1:
# for i in range(len(FL), int(N+1)):
# FL.append(FL[i-1]*i)
# return FL[:int(N) + 1]
# else:
# return FL
def three_j(j,m):
j_1,j_2,j_3=j
m_1,m_2,m_3=m
# symmetry conditions
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, error in three_j)")
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, error in three_m")
if m_1 + m_2 + m_3 != 0:
return 0
PF= np.int((-1) ** int(j_1 - j_2 - m_3))
M=-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
# special case identities, taken from the Mathematica website for 3-j symbols
if ( (j_1==j_2) & (j_3==0) & (m_1==-m_2) & (m_3==0) ):
return (-1)**(j_1-m_1)/np.sqrt(2*j_1+1)
#if ( (m_1==0) & (m_2==0) & (m_3==0)
max_factor=max(j_1 + j_2 + j_3 + 1, j_1 + abs(m_1), j_2 + abs(m_2),
j_3 + abs(m_3))
FL=factorial_list(max_factor)
Sqrt_Arg=(FL[int(j_1+j_2-j_3)] \
*FL[int(j_1-j_2+j_3)] \
*FL[int(-j_1+j_2+j_3)] \
*FL[int(j_1-m_1)] \
*FL[int(j_1+m_1)] \
*FL[int(j_2-m_2)] \
*FL[int(j_2+m_2)] \
*FL[int(j_3-m_3)] \
*FL[int(j_3+m_3)] )/FL[int(j_1+j_2+j_3+1)]
Sqrt_part=np.sqrt(Sqrt_Arg)
# need to fix this
#if Sqrt_part.is_complex:
# Sqrt_part=Sqrt_part.as_real_imag()[0]
i_min = max(-j_3 + j_1 + m_2, -j_3 + j_2 - m_1, 0)
i_max = min(j_2 + m_2, j_1 - m_1, j_1 + j_2 - j_3)
Sum=0
for i in range(int(i_min),int(i_max) + 1):
denom=FL[i] \
*FL[int(i+j_3-j_1-m_2)] \
*FL[int(j_2 + m_2- i)] \
*FL[int(j_1-i-m_1)] \
*FL[int(i + j_3-j_2 + m_1)] \
*FL[int(j_1 +j_2 -j_3 -i)]
Sum=Sum+np.int((-1)**i)/float(denom)
# might have to reset FL
return Sum*Sqrt_part*PF
def Delta_coef(a,b,c,prec=None):
if int(a + b- c) != (a + b - c):
raise ValueError("j values must be integer or half integer and fulfill the triangle relation")
if int(a + c - b) != (a + c - b):
raise ValueError("j values must be integer or half integer and fulfill the triangle relation")
if int(b + c - a) != (b + c - a):
raise ValueError("j values must be integer or half integer and fulfill the triangle relation")
if (a+ b - c) < 0:
return 0
if (a + c - b) < 0:
return 0
if (b + c - a) < 0:
return 0
max_factor=max(a + b -c, a+c-b, b+c -a, a+b +c + 1)
FL=factorial_list(max_factor)
Sqrt_Arg=float(FL[int(a + b -c)] \
*FL[int(a + c - b)] \
*FL[int(b + c - a)])/float(FL[int(a + b + c + 1)])
Sqrt_part=np.sqrt(Sqrt_Arg)
#if prec:
# Sqrt_part=Sqrt_part.evalf(
return Sqrt_part
def Racah(a,b,c,d,e,f, prec=None):
# the Racah symbol
PF=Delta_coef(a,b,e,prec) \
*Delta_coef(c,d,e,prec) \
*Delta_coef(a,c,f,prec) \
*Delta_coef(b,d,f, prec)
if PF==0:
return 0
i_min = max(a + b + e, c + d + e, a + c + f, b + d + f)
i_max = min(a + b+ c + d, a + d + e + f, b + c + e + f)
max_factor=max(i_max + 1, a + b + c + d, a + d + e + f,b + c + e + f)
FL=factorial_list(max_factor)
Sum=0
for i in range(int(i_min), int(i_max)+ 1):
denom=FL[int(i-a-b-e)]\
*FL[int(i-c-d-e)]\
*FL[int(i-a-c-f)]\
*FL[int(i-b-d-f)]\
*FL[int(a + b + c + d - i)]\
*FL[int(a + d + e + f - i)]\
*FL[int(b + c + e + f - i)]
Sum=Sum+((-1)**i*FL[i+1])/float(denom)
return PF*Sum*(-1)**int(a+b+c+d)
def six_j(j):
j_1,j_2,j_3,j_4,j_5,j_6=j
return (-1)**int(j_1+j_2+j_4 +j_5)*Racah(j_1, j_2, j_5, j_4, j_3, j_6,)
if __name__=="__main__":
# j=np.array([2,6,4])
# m=np.array([0,0,0])
# j_1,j_2,j_3=j
# m_1,m_2,m_3=m
# W=wigner.wigner_3j(j_1,j_2,j_3, m_1,m_2,m_3)
#
# print(W
#
# print(three_j(j,m)
# j=np.array([3,3,3,3,3,3])
# print('my six j',six_j(j)
# W=wigner.wigner_6j(3,3,3,3,3,3)
# print('6 j',W, -1/14.
#
# print(six_j(np.array([5,5,5,5,5,5])), 1/52.
print('some test cases for the 3-j symbols')
print('test 1')
print('----------------------------------------')
print('j=1,2,3 & m=0,0,0 => ', -np.sqrt(3/35.))
j=np.array([1,2,3]); m=np.array([0,0,0])
print('check' , three_j(j,m))
print('test 2')
print('----------------------------------------')
print('j=4,5,9 & m=0,0,0 => ', -21*np.sqrt(2/46189.))
j=np.array([4,5,9]); m=np.array([0,0,0])
print('check' , three_j(j,m))
print('test 3')
print('----------------------------------------')
print('j=4,5,6 & m=1,0,-1 => ', -2*np.sqrt(1/429.))
j=np.array([4,5,6]); m=np.array([1,0,-1])
print('check' , three_j(j,m))
print('some test cases for the 6-j symbols')
print('test 1')
print('----------------------------------------')
print('j=4,4,4,4,4,4 => ', -467/18018.)
j=np.array([4,4,4,4,4,4])
print('check' , six_j(j))
print('test 2')
print('----------------------------------------')
print('j=1,1,1,1,1,1 => ', 1/6.)
j=np.array([1,1,1,1,1,1])
print('check' , six_j(j))
print('test 3')
print('----------------------------------------')
print('j=1,2,3,1,3,1 => ', 1/5.)
j=np.array([1,2,3,1,2,1])
print('check' , six_j(j))
Computing file changes ...