From 3bfa2d679cb0f95998ed9746d46fa21e8f2dc9ad Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 11 Apr 2020 11:36:37 +0200 Subject: [PATCH] simpler/correct logic for eu2om --- python/damask/_rotation.py | 24 ++++++++++++------------ src/rotations.f90 | 11 +++++------ 2 files changed, 17 insertions(+), 18 deletions(-) diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index 3336eb6c0..e58de982f 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -474,9 +474,9 @@ class Rotation: q03 = qu[0]**2+qu[3]**2 q12 = qu[1]**2+qu[2]**2 chi = np.sqrt(q03*q12) - if np.abs(chi)< 1.e-6: + if np.abs(q12) < 1.e-6: eu = np.array([np.arctan2(-P*2.0*qu[0]*qu[3],qu[0]**2-qu[3]**2), 0.0, 0.0]) - elif np.abs(q12)< 1.e-6: + elif np.abs(q03) < 1.e-6: eu = np.array([np.arctan2( 2.0*qu[1]*qu[2],qu[1]**2-qu[2]**2), np.pi, 0.0]) else: eu = np.array([np.arctan2((-P*qu[0]*qu[2]+qu[1]*qu[3])*chi, (-P*qu[0]*qu[1]-qu[2]*qu[3])*chi ), @@ -491,14 +491,14 @@ class Rotation: q12_s = qu[...,1:2]**2+qu[...,2:3]**2 chi = np.sqrt(q03_s*q12_s) - eu = np.where(np.abs(chi) < 1.0e-6, + eu = np.where(np.abs(q12_s) < 1.0e-6, np.block([np.arctan2(-P*2.0*qu[...,0:1]*qu[...,3:4],qu[...,0:1]**2-qu[...,3:4]**2), np.zeros(qu.shape[:-1]+(2,))]), np.block([np.arctan2((-P*q02+q13)*chi, (-P*q01-q23)*chi), np.arctan2( 2.0*chi, q03_s-q12_s ), np.arctan2(( P*q02+q13)*chi, (-P*q01+q23)*chi)]) ) - eu = np.where(np.logical_and(np.abs(q03_s) < 1.0e-6, np.abs(chi) > 1.0e-6), + eu = np.where(np.logical_and(np.abs(q03_s) < 1.0e-6, np.abs(q12_s) >= 1.0e-6), np.block([np.arctan2( 2.0*qu[...,1:2]*qu[...,2:3],qu[...,1:2]**2-qu[...,2:3]**2), np.ones( qu.shape[:-1]+(1,))*np.pi, np.zeros(qu.shape[:-1]+(1,))]), @@ -550,7 +550,7 @@ class Rotation: ro = np.where(np.abs(s) < 1.0e-12, [0.0,0.0,P,0.0], np.block([qu[...,1:2]/s,qu[...,2:3]/s,qu[...,3:4]/s, - np.tan(np.arccos(np.clip(qu[:,0:1],-1.0,1.0))) + np.tan(np.arccos(np.clip(qu[...,0:1],-1.0,1.0))) ]) ) ro = np.where(np.abs(qu[...,0:1]) < 1.0e-12, @@ -686,10 +686,10 @@ class Rotation: ee = 0.5*eu cPhi = np.cos(ee[1]) sPhi = np.sin(ee[1]) - qu = np.array([ cPhi*np.cos(ee[0]+ee[2]), - -P*sPhi*np.cos(ee[0]-ee[2]), - -P*sPhi*np.sin(ee[0]-ee[2]), - -P*cPhi*np.sin(ee[0]+ee[2]) ]) + qu = np.array([ cPhi*np.cos(ee[0]+ee[2]), + -P*sPhi*np.cos(ee[0]-ee[2]), + -P*sPhi*np.sin(ee[0]-ee[2]), + -P*cPhi*np.sin(ee[0]+ee[2]) ]) if qu[0] < 0.0: qu*=-1 else: ee = 0.5*eu @@ -776,9 +776,9 @@ class Rotation: ro[3] = np.tan(ro[3]*0.5) else: ax = Rotation.eu2ax(eu) - ro = np.block([ax[:,:3],np.tan(ax[:,3:4]*.5)]) - ro[ax[:,3]>=np.pi,3] = np.inf - ro[np.abs(ax[:,3])<1.e-16] = [ 0.0, 0.0, P, 0.0 ] + ro = np.block([ax[...,:3],np.tan(ax[...,3:4]*.5)]) + ro[ax[...,3]>=np.pi,3] = np.inf + ro[np.abs(ax[...,3])<1.e-16] = [ 0.0, 0.0, P, 0.0 ] return ro @staticmethod diff --git a/src/rotations.f90 b/src/rotations.f90 index 7ce366f74..1630326ad 100644 --- a/src/rotations.f90 +++ b/src/rotations.f90 @@ -432,18 +432,17 @@ pure function qu2eu(qu) result(eu) real(pReal), intent(in), dimension(4) :: qu real(pReal), dimension(3) :: eu - real(pReal) :: q12, q03, chi, chiInv + real(pReal) :: q12, q03, chi q03 = qu(1)**2+qu(4)**2 q12 = qu(2)**2+qu(3)**2 chi = sqrt(q03*q12) - degenerated: if (dEq0(chi)) then - eu = merge([atan2(-P*2.0_pReal*qu(1)*qu(4),qu(1)**2-qu(4)**2), 0.0_pReal, 0.0_pReal], & - [atan2( 2.0_pReal*qu(2)*qu(3),qu(2)**2-qu(3)**2), PI, 0.0_pReal], & - dEq0(q12)) + degenerated: if (dEq0(q12)) then + eu = [atan2(-P*2.0_pReal*qu(1)*qu(4),qu(1)**2-qu(4)**2), 0.0_pReal, 0.0_pReal] + elseif (dEq0(q03)) then + eu = [atan2( 2.0_pReal*qu(2)*qu(3),qu(2)**2-qu(3)**2), PI, 0.0_pReal] else degenerated - chiInv = 1.0_pReal/chi eu = [atan2((-P*qu(1)*qu(3)+qu(2)*qu(4))*chi, (-P*qu(1)*qu(2)-qu(3)*qu(4))*chi ), & atan2( 2.0_pReal*chi, q03-q12 ), & atan2(( P*qu(1)*qu(3)+qu(2)*qu(4))*chi, (-P*qu(1)*qu(2)+qu(3)*qu(4))*chi )]