Sophegrid expansion

VERSION

@@ -1 +1 @@
-v1.1.3
+v1.1.4

deerlab/utils.py

@@ -785,7 +785,7 @@ def read_pickle(filename):
 
 
 # --------------------------------------------------------------------------------------
-def sophegrid(octants,maxphi,size):
+def sophegrid(octants,maxphi,size,closed_phi=False):
     """
     Construct spherical grid over spherical angles based on input parameters. 
     The grid implemented in this function is often called the SOPHE grid [1]_. 
@@ -799,6 +799,9 @@ def sophegrid(octants,maxphi,size):
         Largest value of angle phi (radians).
     size : integer  
         Number of orientations between theta=0 and theta=pi/2. 
+    closed_phi : bool
+        Set to true if grid point at maxPhi should be included, false otherwise. Default is false.
+
 
     Returns
     -------
@@ -834,7 +837,10 @@ def sophegrid(octants,maxphi,size):
         weights = np.zeros(nOrientations)
 
         sindth2 = np.sin(dtheta/2)
-        w1 = 1.0
+        if closed_phi:
+            w1=0.5
+        else:
+            w1 = 1.0
 
         # North pole (z orientation)
         phi[0] = 0
@@ -861,10 +867,11 @@ def sophegrid(octants,maxphi,size):
         weights[idx] = sindth2*dPhi*np.concatenate([[w1], np.ones(nPhi-2), [0.5]])
 
         # Border removal
-        rmv = np.cumsum(nOct*np.arange(1,size)+1)
-        phi = np.delete(phi,rmv)
-        theta = np.delete(theta,rmv)
-        weights = np.delete(weights,rmv)
+        if not closed_phi:
+            rmv = np.cumsum(nOct*np.arange(1,size)+1)
+            phi = np.delete(phi,rmv)
+            theta = np.delete(theta,rmv)
+            weights = np.delete(weights,rmv)
 
         # For C1, add lower hemisphere
         if octants==8:

docsrc/source/changelog.rst

@@ -24,6 +24,10 @@ Release Notes
 - |fix| : Something which was not working as expected or leading to errors has been fixed.
 - |api| : This will require changes in your scripts or code.
 
+Release ``v1.1.4`` - tba
+------------------------------------------
+- |enhancement| : Expanded sophgrid to allow for closed phi integral. (:pr:`482`)
+
 Release ``v1.1.3`` - July 2024
 ------------------------------------------
 - |fix| : Removes unnecessary files from the docs

test/test_utils.py

@@ -13,4 +13,11 @@ def test_sophegrid():
     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]))
+    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]))
+
+    phi, theta, weights = sophegrid(1,np.pi/2,5,closed_phi=True)
+
+    assert np.allclose(weights.sum(),1)
+    assert np.allclose(phi, np.array([0, 0, 1.5708, 0, 0.7854, 1.5708, 0, 0.5236, 1.0472, 1.5708, 0, 0.3927, 0.7854, 1.1781, 1.5708]),rtol=1e-4)
+    assert np.allclose(theta, np.array([0, 0.3927, 0.3927, 0.7854, 0.7854, 0.7854, 1.1781, 1.1781, 1.1781, 1.1781, 1.5708, 1.5708, 1.5708, 1.5708, 1.5708]),rtol=1e-4)
+    assert np.allclose(weights*4*np.pi, np.array([0.24146, 0.93818, 0.93818, 0.86676, 1.7335, 0.86676, 0.75499, 1.51, 1.51, 0.75499, 0.30645, 0.61289, 0.61289, 0.61289, 0.30645]),rtol=1e-4)