fixed spherical component sampling and testing

python/damask/_rotation.py

@@ -639,7 +639,7 @@ class Rotation:
 
 
     @staticmethod
-    def from_spherical_component(center,FWHM,N_samples=500,degrees=True,seed=None):
+    def from_spherical_component(center,sigma,N=500,degrees=True,seed=None):
         """
         Calculate set of rotations with Gaussian distribution around center.
 
@@ -652,25 +652,26 @@ class Rotation:
         ----------
         center : Rotation
             Central Rotation.
-        FWHM : float
+        sigma : float
-            Full width at half maximum of the Gaussian distribution.
+            Standard deviation of (Gaussian) misorientation distribution.
-        N_samples : int, optional
+        N : int, optional
             Number of samples, defaults to 500.
         degrees : boolean, optional
-            FWHM is given in degrees.
+            sigma is given in degrees.
         seed : {None, int, array_like[ints], SeedSequence, BitGenerator, Generator}, optional
-            A seed to initialize the BitGenerator. Defaults to None.
+            A seed to initialize the BitGenerator. Defaults to None, i.e. unpredictable entropy
-            If None, then fresh, unpredictable entropy will be pulled from the OS.
+            will be pulled from the OS.
 
         """
         rng = np.random.default_rng(seed)
-        sigma = (np.radians(FWHM) if degrees else FWHM)/np.sqrt(2*np.log(2))*.5
+        sigma = np.radians(sigma) if degrees else sigma
-        u = rng.random(N_samples) * 2.0 - 1.0
+        u,Theta = (rng.random((N,2)) * 2.0 * np.array([1,np.pi]) - np.array([1.0, 0])).T
-        Theta = rng.random(N_samples) * 2.0 * np.pi
+        omega = rng.normal(scale=sigma,size=N)
-        omega = rng.normal(scale=sigma,size=N_samples)
+        p = np.column_stack([np.sqrt(1-u**2)*np.cos(Theta),
-        ax = np.column_stack([np.sqrt(1-u**2)*np.cos(Theta),np.sqrt(1-u**2)*np.sin(Theta),u,omega])
+                             np.sqrt(1-u**2)*np.sin(Theta),
-        ax[ax[:,3]<0]*=-1
+                             u, omega])
-        return  Rotation.from_axis_angle(ax) @ center.broadcast_to(N_samples)
+        p[p[:,3]<0] *= -1
+        return  Rotation.from_axis_angle(p) @ center
 
 
     @staticmethod
@@ -705,14 +706,13 @@ class Rotation:
         if np.isclose(ax_align[3],0.0): ax_align[:3] = np.array([1,0,0])
         R_align  = Rotation.from_axis_angle(ax_align if ax_align[3] > 0.0 else -ax_align,normalize=True)          # rotation to align fiber axis in crystal and sample system
 
-        u,v,b = (np.random.random((N,3)) * 2 * np.array([1,np.pi,np.pi]) - np.array([1,0,np.pi])).T
+        u,Theta,b = (np.random.random((N,3)) * 2 * np.array([1,np.pi,np.pi]) - np.array([1,0,np.pi])).T
-        a = abs(np.random.normal(scale=sigma_,size=N))
+        omega = abs(np.random.normal(scale=sigma_,size=N))
-        p = np.vstack((np.sqrt(1-u**2)*np.cos(v),
+        p = np.column_stack([np.sqrt(1-u**2)*np.cos(Theta),
-                       np.sqrt(1-u**2)*np.sin(v),
+                             np.sqrt(1-u**2)*np.sin(Theta),
-                       u,
+                             u, omega])
-                       a)).T
         p[:,:3] = np.einsum('ij,...j->...i',np.eye(3)-np.outer(d_lab,d_lab),p[:,:3])
-        f = np.hstack((np.broadcast_to(d_lab,(N,3)),b.reshape(N,1)))
+        f = np.column_stack((np.broadcast_to(d_lab,(N,3)),b))
         f[f[:,3]<0] *= -1.
         return R_align.broadcast_to(N) \
              @ Rotation.from_axis_angle(p,normalize=True) \

python/tests/test_Rotation.py

@@ -914,20 +914,18 @@ class TestRotation:
         assert np.isclose(avg_angle,10+(angle-10)/2.)
 
 
-    @pytest.mark.parametrize('FWHM',[5,10,15])
+    @pytest.mark.parametrize('sigma',[5,10,15,20])
-    @pytest.mark.parametrize('N_samples',[600,1200,2400])
+    @pytest.mark.parametrize('N',[1000,10000,100000])
-    def test_spherical_component(self,N_samples,FWHM):
+    def test_spherical_component(self,N,sigma):
-        """https://en.wikipedia.org/wiki/Full_width_at_half_maximum."""
         c = Rotation.from_random()
-        o = Rotation.from_spherical_component(c,FWHM,N_samples)
+        o = Rotation.from_spherical_component(c,sigma,N)
-        m = c.broadcast_to(N_samples).misorientation(o)
+        _, angles = c.misorientation(o).as_axis_angle(pair=True,degrees=True)
-        _, angles = m.as_axis_angle(pair=True,degrees=True)
+        angles[::2] *= -1                                                                           # flip angle for every second to symmetrize distribution
-        dist = angles * (np.random.randint(0,2,N_samples)*2-1)
 
-        p = stats.normaltest(dist)[1]
+        p = stats.normaltest(angles)[1]
-        FWHM_out = np.std(dist) # * (2*np.sqrt(2*np.log(2)))
+        sigma_out = np.std(angles)
-        print(f'\np: {p}, FWHM ratio {FWHM/FWHM_out}')
+        print(f'\np: {p}, sigma ratio {sigma/sigma_out}')
-        assert (.9 < FWHM/FWHM_out < 1.1) and p > 0.001
+        assert (.9 < sigma/sigma_out < 1.1) and p > 0.001
 
 
     @pytest.mark.parametrize('sigma',[5,10,15,20])
@@ -949,4 +947,4 @@ class TestRotation:
         p = stats.normaltest(dist)[1]
         sigma_out = np.degrees(np.std(dist))
         print(f'\np: {p}, sigma ratio {sigma/sigma_out}')
-        assert (.85 < sigma/sigma_out < 1.15) and p > 0.001
+        assert (.9 < sigma/sigma_out < 1.1) and p > 0.001