fixed fiber component sampling and testing

python/damask/_rotation.py

@@ -687,55 +687,49 @@ class Rotation:
 
 
     @staticmethod
-    def from_fiber_component(alpha,beta,FWHM=0.0,N=500,degrees=True,seed=None):
+    def from_fiber_component(alpha,beta,sigma=0.0,N=500,degrees=True,seed=None):
         """
         Calculate set of rotations with Gaussian distribution around direction.
 
-        References
-        ----------
-        K. Helming, Texturapproximation durch Modellkomponenten
-        Cuvillier Verlag, 1996
-
         Parameters
         ----------
         alpha : numpy.ndarray of size 2
-            tbd.
+            Polar coordinates (phi from x,theta from z) of fiber direction in crystal frame.
         beta : numpy.ndarray of size 2
-            tbd.
+            Polar coordinates (phi from x,theta from z) of fiber direction in sample frame.
-        FWHM : float, optional
+        sigma : float, optional
-            Full width at half maximum of the Gaussian distribution.
+            Standard deviation of (Gaussian) misorientation distribution.
             Defaults to 0.
-        N_samples : int, optional
+        N : int, optional
             Number of samples, defaults to 500.
         degrees : boolean, optional
-            FWHM, alpha, and beta are given in degrees.
+            sigma, alpha, and beta are 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)
-        FWHM_,alpha_,beta_ = map(np.radians,(FWHM,alpha,beta)) if degrees else (FWHM,alpha,beta)
+        sigma_,alpha_,beta_ = map(np.radians,(sigma,alpha,beta)) if degrees else (sigma,alpha,beta)
 
         d_cr  = np.array([np.sin(alpha_[0])*np.cos(alpha_[1]), np.sin(alpha_[0])*np.sin(alpha_[1]), np.cos(alpha_[0])])
-        d_lab = np.array([np.sin(beta_[0] )*np.cos(beta_[1] ), np.sin(beta_[0] )*np.sin(beta_[1] ), np.cos(beta_[0] )])
+        d_lab = np.array([np.sin( beta_[0])*np.cos( beta_[1]), np.sin( beta_[0])*np.sin( beta_[1]), np.cos( beta_[0])])
         ax_align = np.append(np.cross(d_lab,d_cr), np.arccos(np.dot(d_lab,d_cr)))
         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
-        a = abs(np.random.normal(scale=FWHM_,size=N))
+        a = abs(np.random.normal(scale=sigma_,size=N))
         p = np.vstack((np.sqrt(1-u**2)*np.cos(v),
                        np.sqrt(1-u**2)*np.sin(v),
                        u,
                        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_cr,(N,3)),b.reshape(N,1)))
+        f = np.hstack((np.broadcast_to(d_lab,(N,3)),b.reshape(N,1)))
         f[f[:,3]<0] *= -1.
-        return Rotation.from_axis_angle(p) \
+        return R_align.broadcast_to(N) \
-             * Rotation.from_axis_angle(f) \
+             @ Rotation.from_axis_angle(p,normalize=True) \
-             * R_align
+             @ Rotation.from_axis_angle(f)
-
 
 
 ####################################################################################################

python/tests/test_Rotation.py

@@ -925,14 +925,14 @@ class TestRotation:
         dist = angles * (np.random.randint(0,2,N_samples)*2-1)
 
         p = stats.normaltest(dist)[1]
-        FWHM_out = np.std(dist) * (2*np.sqrt(2*np.log(2)))
+        FWHM_out = np.std(dist) # * (2*np.sqrt(2*np.log(2)))
         print(f'\np: {p}, FWHM ratio {FWHM/FWHM_out}')
         assert (.9 < FWHM/FWHM_out < 1.1) and p > 0.001
 
 
-    @pytest.mark.parametrize('FWHM',[10,15,20])
+    @pytest.mark.parametrize('sigma',[5,10,15,20])
-    @pytest.mark.parametrize('N_samples',[500,1000,2000])
+    @pytest.mark.parametrize('N',[1000,10000,100000])
-    def test_from_fiber_component(self,N_samples,FWHM):
+    def test_from_fiber_component(self,N,sigma):
         """https://en.wikipedia.org/wiki/Full_width_at_half_maximum."""
         alpha = np.random.random(2)*np.pi
         beta  = np.random.random(2)*np.pi
@@ -942,12 +942,11 @@ class TestRotation:
         ax = np.append(np.cross(f_in_C,f_in_S), - np.arccos(np.dot(f_in_C,f_in_S)))
         n = Rotation.from_axis_angle(ax if ax[3] > 0.0 else ax*-1.0 ,normalize=True)                # rotation to align fiber axis in crystal and sample system
 
-        o = Rotation.from_fiber_component(alpha,beta,np.radians(FWHM),N_samples,False)
+        o = Rotation.from_fiber_component(alpha,beta,np.radians(sigma),N,False)
-        angles=[np.arccos(np.clip(np.dot(n@f_in_S,o[i]@f_in_S),-1,1)) for i in range(N_samples)]
+        angles = np.arccos(np.clip(np.dot(o@np.broadcast_to(f_in_S,(N,3)),n@f_in_S),-1,1))
-        dist = np.array(angles) * (np.random.randint(0,2,N_samples)*2-1)
+        dist   = np.array(angles) * (np.random.randint(0,2,N)*2-1)
 
         p = stats.normaltest(dist)[1]
-        FWHM_out = np.degrees(np.std(dist)) * (2*np.sqrt(2*np.log(2)))
+        sigma_out = np.degrees(np.std(dist))
-        print(f'\np: {p}, FWHM ratio {FWHM/FWHM_out}')
+        print(f'\np: {p}, sigma ratio {sigma/sigma_out}')
-        assert (.85 < FWHM/FWHM_out < 1.15) and p > 0.001
+        assert (.85 < sigma/sigma_out < 1.15) and p > 0.001
-