https://github.com/JoeMcEwen/FAST-PT
Revision c56761670c0e4c154b7f2794c79c06c840ba5b0a authored by Jonathan Blazek on 01 March 2019, 22:19:49 UTC, committed by Jonathan Blazek on 01 March 2019, 22:19:49 UTC
1 parent 137418e
Tip revision: c56761670c0e4c154b7f2794c79c06c840ba5b0a authored by Jonathan Blazek on 01 March 2019, 22:19:49 UTC
merging changes from develop branch to prepare for merge with dev. develop will be deprecated.
merging changes from develop branch to prepare for merge with dev. develop will be deprecated.
Tip revision: c567616
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 ...