diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index 6a94405ee..425e7aeef 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -663,8 +663,8 @@ class Rotation: """ rng = np.random.default_rng(seed) - _FWHM = np.radians(FWHM) if degrees else FWHM - if _FWHM > np.radians(0.1): + FWHM_ = np.radians(FWHM) if degrees else FWHM + if FWHM_ > np.radians(0.1): rotations = [] for i in range(N_samples): while True: @@ -673,11 +673,10 @@ class Rotation: ax = np.array([a, np.sqrt(1.0-a**2.0)*np.cos(rnd[1]*2.0*np.pi), # random axis np.sqrt(1.0-a**2.0)*np.sin(rnd[1]*2.0*np.pi), # random axis - (rnd[2]-0.5) * 4 *_FWHM]) # rotation by [0, +-2 FWHM] - if ax[3] < 0.0: ax*=-1.0 - R = Rotation.from_axis_angle(ax,normalize=True) + (rnd[2]-0.5) * 4 *FWHM_]) # rotation by [0, +-2 FWHM] + R = Rotation.from_axis_angle(ax if ax[3] > 0.0 else ax*-1.0,normalize=True) angle = R.misorientation(Rotation()).as_axis_angle()[3] # misorientation to unity - if rnd[3] <= np.exp(-4.0*np.log(2.0)*(angle/_FWHM)**2): # rejection sampling (Gaussian) + if rnd[3] <= np.exp(-4.0*np.log(2.0)*(angle/FWHM_)**2): # rejection sampling (Gaussian) break rotations.append((center @ R).as_quaternion()) else: @@ -714,21 +713,19 @@ class Rotation: """ rng = np.random.default_rng(seed) - FWHM_,alpha_,beta_ = map(np.radians,(FWHM,alpha,beta)) if degrees else (FWHM,alpha,beta) + FWHM_,alpha_,beta_ = np.radians((FWHM,alpha,beta)) if degrees else (FWHM,alpha,beta) f_in_C = np.array([np.sin(alpha_[0])*np.cos(alpha_[1]), np.sin(alpha_[0])*np.sin(alpha_[1]), np.cos(alpha_[0])]) f_in_S = np.array([np.sin(beta_[0] )*np.cos(beta_[1] ), np.sin(beta_[0] )*np.sin(beta_[1] ), np.cos(beta_[0] )]) + ax = np.append(np.cross(f_in_C,f_in_S), - np.arccos(np.dot(f_in_C,f_in_S))) + R_align = 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 rotations = [] for i in range(N_samples): rnd = rng.random(3) - ax = np.append(np.cross(f_in_C,f_in_S), - np.arccos(np.dot(f_in_C,f_in_S))) - if ax[3] < 0.0: ax *= -1.0 - R = Rotation.from_axis_angle(ax,normalize=True) # rotation to align fiber axis in crystal and sample system - ax = np.append(f_in_S,rnd[0]*np.pi*2.0) if ax[3] > np.pi: ax = np.append(ax[:3]*-1,np.pi*2-ax[3]) - R = R @ Rotation.from_axis_angle(ax) # rotation (0..360deg) perpendicular to fiber axis + R = R_align @ Rotation.from_axis_angle(ax) # rotation (0..360deg) perpendicular to fiber axis if FWHM_ > np.radians(0.1): @@ -737,24 +734,21 @@ class Rotation: s = f_in_S[i_smallest] a = f_in_S[i_non_smallest] - x = sum([x**2 for x in a]) + x = np.sum(a**2) u = np.empty(3) while True: # rejection sampling angle = (rnd[1] - 0.5)*4 *FWHM_ # solve cos(angle) = dot_product(fInS,u) for u. This is underdetermined, hence assume that # they share the smallest component. - c = np.cos(angle) - s**2 - u[i_non_smallest[1]] = -(2.0*c*a[1] + np.sqrt(4*((c*a[1])**2.0-x*(c**2.0-a[0]**2*(1.0-s**2)))))/(2*x) u[i_non_smallest[0]] = np.sqrt(1.0-u[i_non_smallest[1]]**2.0-s**2.0) u[i_smallest] = s if (rnd[2] <= np.exp(-4.0*np.log(2.0)*(angle/FWHM_)**2)): ax = np.append(np.cross(u,f_in_S),angle) - if ax[3]<0.0: ax *= -1.0 - R = R * Rotation.from_axis_angle(ax,normalize=True) # tilt around direction of smallest component + R = R @ Rotation.from_axis_angle(ax if ax[3] > 0.0 else ax*-1.0,normalize=True) # tilt around direction of smallest component break else: rnd = rng.random(3) diff --git a/python/tests/test_Rotation.py b/python/tests/test_Rotation.py index 409e569c0..b00f661ae 100644 --- a/python/tests/test_Rotation.py +++ b/python/tests/test_Rotation.py @@ -934,17 +934,23 @@ class TestRotation: @pytest.mark.parametrize('N_samples',[500,1000,2000]) def test_from_fiber_component(self,N_samples,FWHM): """https://en.wikipedia.org/wiki/Full_width_at_half_maximum.""" - alpha = np.array([15.0,4.6]) - beta = np.ones(2) - n = Rotation.from_quaternion([0.9914448613738086,-0.01046806021254377,0.13010575149028156,7.146345858741878e-08]) - o = Rotation.from_fiber_component(alpha,beta,FWHM,N_samples,True) + alpha = np.random.random(2)*np.pi + beta = np.zeros(2) + + f_in_C = np.array([np.sin(alpha[0])*np.cos(alpha[1]), np.sin(alpha[0])*np.sin(alpha[1]), np.cos(alpha[0])]) + f_in_S = np.array([np.sin(beta[0] )*np.cos(beta[1] ), np.sin(beta[0] )*np.sin(beta[1] ), np.cos(beta[0] )]) + 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) angles=[] for i in range(N_samples): - cos = np.dot(np.dot(n.as_matrix(),np.array([0.0,0.0,1.0])), - np.dot(o[i].as_matrix(),np.array([0.0, 0.0, 1.0]))) + cos = np.dot(n@np.array([0.0,0.0,1.0]),o[i]@np.array([0.0, 0.0, 1.0])) angles.append(np.arccos(np.clip(cos,-1,1))) dist = np.array(angles) * (np.random.randint(0,2,N_samples)*2-1) + p = stats.normaltest(dist)[1] FWHM_out = np.degrees(np.std(dist)) * (2*np.sqrt(2*np.log(2))) - print(f'\n FWHM ratio {FWHM/FWHM_out}') - assert .85 < FWHM/FWHM_out < 1.1 + print(f'\np: {p}, FWHM ratio {FWHM/FWHM_out}') + assert (.85 < FWHM/FWHM_out < 1.15) and p > 0.001 +