diff --git a/arc/wigner.py b/arc/wigner.py index 5fac5a8..0d9491f 100644 --- a/arc/wigner.py +++ b/arc/wigner.py @@ -187,7 +187,7 @@ def Wigner3j(j1, j2, j3, m1, m2, m3): ) -def Wigner6j(j1, j2, j3, J1, J2, J3): +def Wigner6j(j1, j2, j3, J1, J2, J3, verbose=False): r""" Evaluates Wigner 6-j symbol @@ -228,26 +228,52 @@ def Wigner6j(j1, j2, j3, J1, J2, J3): # Check if the 4 triads ( (j1 j2 j3), (j1 J2 J3), (J1 j2 J3), (J1 J2 j3) ) # satisfy the triangular inequalities - if ( - (abs(j1 - j2) > j3) - | (j1 + j2 < j3) - | (abs(j1 - J2) > J3) - | (j1 + J2 < J3) - | (abs(J1 - j2) > J3) - | (J1 + j2 < J3) - | (abs(J1 - J2) > j3) - | (J1 + J2 < j3) - ): + IsTriangle = True + msg = '' + if ((abs(j1 - j2) > j3) | (j1 + j2 < j3)): + IsTriangle = False + msg += '(%.1f, %.1f, %.1f) is not triangular\n'%(j1,j2,j3) + if ((abs(j1 - J2) > J3) | (j1 + J2 < J3)): + IsTriangle = False + msg += '(%.1f, %.1f, %.1f) is not triangular\n'%(j1,J2,J3) + if ((abs(J1 - j2) > J3) | (J1 + j2 < J3)): + IsTriangle = False + msg += '(%.1f, %.1f, %.1f) is not triangular\n'%(J1,j2,J3) + if ((abs(J1 - J2) > j3) | (J1 + J2 < j3)): + IsTriangle = False + msg += '(%.1f, %.1f, %.1f) is not triangular\n'%(J1,J2,j3) + if not IsTriangle: + msg = "WARNING!!\n" + msg + msg += "For the 6j-Symbol:\n ⎰%3.1f %3.1f %3.1f⎱\n ⎱%3.1f %3.1f %3.1f⎰"%(j1,j2,j3,J1,J2,J3) + if verbose: + print(msg) + return 0 raise ValueError("6j-Symbol is not triangular!") # Check if the sum of the elements of each traid is an integer - if ( - (2 * (j1 + j2 + j3) != roundPy2(2 * (j1 + j2 + j3))) - | (2 * (j1 + J2 + J3) != roundPy2(2 * (j1 + J2 + J3))) - | (2 * (J1 + j2 + J3) != roundPy2(2 * (J1 + j2 + J3))) - | (2 * (J1 + J2 + j3) != roundPy2(2 * (J1 + J2 + j3))) - ): - raise ValueError("6j-Symbol is not triangular!") + SumIsInteger = True + msg = '' + if (2 * roundPy2(j1 + j2 + j3) != roundPy2(2 * (j1 + j2 + j3))): + SumIsInteger = False + msg += '%.1f + %.1f + %.1f is not an integer\n'%(j1,j2,j3) + if (2 * roundPy2(j1 + J2 + J3) != roundPy2(2 * (j1 + J2 + J3))): + SumIsInteger = False + msg += '%.1f + %.1f + %.1f is not an integer\n'%(j1,J2,J3) + if (2 * roundPy2(J1 + j2 + J3) != roundPy2(2 * (J1 + j2 + J3))): + SumIsInteger = False + msg += '%.1f + %.1f + %.1f is not an integer\n'%(J1,j2,J3) + if (2 * roundPy2(J1 + J2 + j3) != roundPy2(2 * (J1 + J2 + j3))): + SumIsInteger = False + msg += '%.1f + %.1f + %.1f is not an integer\n'%(J1,J2,j3) + if not SumIsInteger: + msg = "WARNING!!\n" + msg + msg += "For the 6j-Symbol:\n ⎰%3.1f %3.1f %3.1f⎱\n ⎱%3.1f %3.1f %3.1f⎰"%(j1,j2,j3,J1,J2,J3) + msg += "\n6j-Symbol is undefined when any triad has a non-integer sum" + if verbose: + print(msg) + return np.nan + raise ValueError("6j-Symbol is undefined when any triad has a non-integer sum") + raise ValueError("6j-Symbol is defined only when all triads have integer sums") # if possible, use precalculated values global wignerPrecal