From ad2badd3bebefcb73e01fbc5bc6c1315b26ac9e1 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 15 Sep 2020 22:13:28 +0200 Subject: [PATCH] [skip ci] vectorized and simplified based on Philips ideas. Test requires from_axis_angle fix --- python/damask/_rotation.py | 27 +++++++-------------------- 1 file changed, 7 insertions(+), 20 deletions(-) diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index 425e7aeef..a43c79e24 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -663,26 +663,13 @@ class Rotation: """ rng = np.random.default_rng(seed) - FWHM_ = np.radians(FWHM) if degrees else FWHM - if FWHM_ > np.radians(0.1): - rotations = [] - for i in range(N_samples): - while True: - rnd = rng.random(4) - a = rnd[0]*2.0 -1.0 # uniform - 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] - 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) - break - rotations.append((center @ R).as_quaternion()) - else: - rotations = [center.as_quaternion() for i in range(N_samples)] - - return Rotation.from_quaternion(rotations) + sigma = (np.radians(FWHM) if degrees else FWHM)/np.sqrt(2*np.log(2))*.5 + u = rng.random(N_samples) * 2.0 - 1.0 + Theta = rng.random(N_samples) * 2.0 * np.pi + omega = rng.normal(scale=sigma,size=N_samples) + ax = np.column_stack([np.sqrt(1-u**2)*np.cos(Theta),np.sqrt(1-u**2)*np.sin(Theta),u,omega]) + ax[ax[:,3]<0]*=-1 + return Rotation.from_axis_angle(ax) @ center.broadcast_to(N_samples) @staticmethod