Merge pull request #148 from petermao/bugfix/Wigner6jIntegerSum

Bugfix/wigner6j integer sum
nikolasibalic · Sep 30, 2023 · d1341a9 · d1341a9
2 parents 63f37a8 + a3682e5
commit d1341a9
Showing 1 changed file with 44 additions and 18 deletions.
diff --git 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