223031012
/
DAMASK_EICMD
mirror of https://github.com/achalhp/DAMASK_EICMD.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity

polishing

This commit is contained in:
Martin Diehl 2020-09-15 12:27:06 +02:00
parent 2a082b7983
commit 18c38f1284
2 changed files with 25 additions and 25 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

28
python/damask/_rotation.py
View File

@ -663,8 +663,8 @@ class Rotation:
""" """
rng = np.random.default_rng(seed) rng = np.random.default_rng(seed)
_FWHM = np.radians(FWHM) if degrees else FWHM FWHM_ = np.radians(FWHM) if degrees else FWHM
if _FWHM > np.radians(0.1): if FWHM_ > np.radians(0.1):
rotations = [] rotations = []
for i in range(N_samples): for i in range(N_samples):
while True: while True:
@ -673,11 +673,10 @@ class Rotation:
ax = np.array([a, 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.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 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] (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 if ax[3] > 0.0 else ax*-1.0,normalize=True)
R = Rotation.from_axis_angle(ax,normalize=True)
angle = R.misorientation(Rotation()).as_axis_angle()[3] # misorientation to unity 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 break
rotations.append((center @ R).as_quaternion()) rotations.append((center @ R).as_quaternion())
else: else:
@ -714,21 +713,19 @@ class Rotation:
""" """
rng = np.random.default_rng(seed) 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_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] )]) 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 = [] rotations = []
for i in range(N_samples): for i in range(N_samples):
rnd = rng.random(3) 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) 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]) 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): if FWHM_ > np.radians(0.1):
@ -737,24 +734,21 @@ class Rotation:
s = f_in_S[i_smallest] s = f_in_S[i_smallest]
a = f_in_S[i_non_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) u = np.empty(3)
while True: # rejection sampling while True: # rejection sampling
angle = (rnd[1] - 0.5)*4 *FWHM_ angle = (rnd[1] - 0.5)*4 *FWHM_
# solve cos(angle) = dot_product(fInS,u) for u. This is underdetermined, hence assume that # solve cos(angle) = dot_product(fInS,u) for u. This is underdetermined, hence assume that
# they share the smallest component. # they share the smallest component.
c = np.cos(angle) - s**2 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[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_non_smallest[0]] = np.sqrt(1.0-u[i_non_smallest[1]]**2.0-s**2.0)
u[i_smallest] = s u[i_smallest] = s
if (rnd[2] <= np.exp(-4.0*np.log(2.0)*(angle/FWHM_)**2)): if (rnd[2] <= np.exp(-4.0*np.log(2.0)*(angle/FWHM_)**2)):
ax = np.append(np.cross(u,f_in_S),angle) 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 if ax[3] > 0.0 else ax*-1.0,normalize=True) # tilt around direction of smallest component
R = R * Rotation.from_axis_angle(ax,normalize=True) # tilt around direction of smallest component
break break
else: else:
rnd = rng.random(3) rnd = rng.random(3)

22
python/tests/test_Rotation.py
View File

@ -934,17 +934,23 @@ class TestRotation:
@pytest.mark.parametrize('N_samples',[500,1000,2000]) @pytest.mark.parametrize('N_samples',[500,1000,2000])
def test_from_fiber_component(self,N_samples,FWHM): def test_from_fiber_component(self,N_samples,FWHM):
"""https://en.wikipedia.org/wiki/Full_width_at_half_maximum.""" """https://en.wikipedia.org/wiki/Full_width_at_half_maximum."""
alpha = np.array([15.0,4.6]) alpha = np.random.random(2)*np.pi
beta = np.ones(2) beta = np.zeros(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) 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=[] angles=[]
for i in range(N_samples): for i in range(N_samples):
cos = np.dot(np.dot(n.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]))
np.dot(o[i].as_matrix(),np.array([0.0, 0.0, 1.0])))
angles.append(np.arccos(np.clip(cos,-1,1))) angles.append(np.arccos(np.clip(cos,-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_samples)*2-1)
p = stats.normaltest(dist)[1]
FWHM_out = np.degrees(np.std(dist)) * (2*np.sqrt(2*np.log(2))) FWHM_out = np.degrees(np.std(dist)) * (2*np.sqrt(2*np.log(2)))
print(f'\n FWHM ratio {FWHM/FWHM_out}') print(f'\np: {p}, FWHM ratio {FWHM/FWHM_out}')
assert .85 < FWHM/FWHM_out < 1.1 assert (.85 < FWHM/FWHM_out < 1.15) and p > 0.001