Fixing Sophgrid bug (#464)

* Fixing Sophgrid bug * Add unit test for sophgrid
JeschkeLab · Oct 5, 2023 · 3c077ff · 3c077ff
1 parent b9b76e4
commit 3c077ff
Show file tree

Hide file tree

Showing 2 changed files with 20 additions and 4 deletions.
diff --git a/deerlab/utils.py b/deerlab/utils.py
@@ -833,15 +833,15 @@ def sophegrid(octants,maxphi,size):
         weights = np.zeros(nOrientations)
 
         sindth2 = np.sin(dtheta/2)
-        w1 = 0.5
+        w1 = 1.0
 
         # North pole (z orientation)
         phi[0] = 0
         theta[0] = 0
         weights[0] = maxphi*(1 - np.cos(dtheta/2))
 
         # All but equatorial slice
-        Start = 2
+        Start = 1
         for iSlice in np.arange(2,size):
             nPhi = nOct*(iSlice-1) + 1
             dPhi = maxphi/(nPhi - 1)
@@ -854,13 +854,13 @@ def sophegrid(octants,maxphi,size):
         # Equatorial slice
         nPhi = nOct*(size - 1) + 1
         dPhi = maxphi/(nPhi - 1)
-        idx = Start + (np.arange(0,nPhi) - 1)
+        idx = Start + np.arange(0,nPhi)
         phi[idx] = np.linspace(0,maxphi,nPhi)
         theta[idx] = np.pi/2
         weights[idx] = sindth2*dPhi*np.concatenate([[w1], np.ones(nPhi-2), [0.5]])
 
         # Border removal
-        rmv = np.cumsum(nOct*np.arange(1,size-1)+1)+1 
+        rmv = np.cumsum(nOct*np.arange(1,size)+1)
         phi = np.delete(phi,rmv)
         theta = np.delete(theta,rmv)
         weights = np.delete(weights,rmv)

diff --git a/test/test_utils.py b/test/test_utils.py
@@ -0,0 +1,16 @@
+import numpy as np
+from deerlab import sophegrid
+
+
+# ======================================================================
+# Test sophegrid function
+
+def test_sophegrid():
+    # Test that the function returns the values and weights summing to 1
+    # Comparison taken from EasySpin 6.0
+    phi, theta, weights = sophegrid(4,np.pi*2,3)
+
+    assert np.allclose(weights.sum(),1)
+    assert np.allclose(phi, np.array([0,0,1.57079632679490,3.14159265358979,4.71238898038469,0,0.785398163397448,1.57079632679490,2.35619449019235,3.14159265358979,3.92699081698724,4.71238898038469,5.49778714378214]))
+    assert np.allclose(theta, np.array([0,0.785398163397448,0.785398163397448,0.785398163397448,0.785398163397448,1.57079632679490,1.57079632679490,1.57079632679490,1.57079632679490,1.57079632679490,1.57079632679490,1.57079632679490,1.57079632679490]))
+    assert np.allclose(weights*4*np.pi, np.array([0.956558005801449,1.70021769237074,1.70021769237074,1.70021769237074,1.70021769237074,0.601117729884346,0.601117729884346,0.601117729884346,0.601117729884346,0.601117729884346,0.601117729884346,0.601117729884346,0.601117729884346]))