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

bugfix: valid range for unit quaternion range is [-1,+1]

This commit is contained in:
Martin Diehl 2019-02-04 00:06:38 +01:00
parent 9a4e9e62b6
commit 5320803842
1 changed files with 6 additions and 6 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

12
src/rotations.f90
View File

@ -353,8 +353,6 @@ end function eu2qu
!> @brief orientation matrix to Euler angles !> @brief orientation matrix to Euler angles
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
pure function om2eu(om) result(eu) pure function om2eu(om) result(eu)
use prec, only: &
dEq
use math, only: & use math, only: &
PI PI
@ -363,7 +361,7 @@ pure function om2eu(om) result(eu)
real(pReal), dimension(3) :: eu real(pReal), dimension(3) :: eu
real(pReal) :: zeta real(pReal) :: zeta
if (dEq(abs(om(3,3)),1.0_pReal,1.0e-15_pReal)) then if (abs(om(3,3))>1.0_pReal) then
eu = [ atan2( om(1,2),om(1,1)), 0.5*PI*(1-om(3,3)),0.0_pReal ] eu = [ atan2( om(1,2),om(1,1)), 0.5*PI*(1-om(3,3)),0.0_pReal ]
else else
zeta = 1.0_pReal/sqrt(1.0_pReal-om(3,3)**2.0_pReal) zeta = 1.0_pReal/sqrt(1.0_pReal-om(3,3)**2.0_pReal)
@ -774,7 +772,7 @@ pure function qu2ax(qu) result(ax)
real(pReal) :: omega, s real(pReal) :: omega, s
omega = 2.0 * acos(math_clip(qu%w,0.0_pReal,1.0_pReal)) omega = 2.0 * acos(math_clip(qu%w,-1.0_pReal,1.0_pReal))
! if the angle equals zero, then we return the rotation axis as [001] ! if the angle equals zero, then we return the rotation axis as [001]
if (dEq0(omega)) then if (dEq0(omega)) then
ax = [ 0.0, 0.0, 1.0, 0.0 ] ax = [ 0.0, 0.0, 1.0, 0.0 ]
@ -798,6 +796,8 @@ pure function qu2ro(qu) result(ro)
IEEE_positive_inf IEEE_positive_inf
use prec, only: & use prec, only: &
dEq0 dEq0
use math, only: &
math_clip
type(quaternion), intent(in) :: qu type(quaternion), intent(in) :: qu
real(pReal), dimension(4) :: ro real(pReal), dimension(4) :: ro
@ -810,7 +810,7 @@ pure function qu2ro(qu) result(ro)
else else
s = norm2([qu%x,qu%y,qu%z]) s = norm2([qu%x,qu%y,qu%z])
ro = merge ( [ 0.0_pReal, 0.0_pReal, P, 0.0_pReal], & ro = merge ( [ 0.0_pReal, 0.0_pReal, P, 0.0_pReal], &
[ qu%x/s, qu%y/s, qu%z/s, tan(acos(qu%w))], & [ qu%x/s, qu%y/s, qu%z/s, tan(acos(math_clip(qu%w,-1.0_pReal,1.0_pReal))], &
s < thr) !ToDo: not save (PGI compiler) s < thr) !ToDo: not save (PGI compiler)
end if end if
@ -833,7 +833,7 @@ pure function qu2ho(qu) result(ho)
real(pReal) :: omega, f real(pReal) :: omega, f
omega = 2.0 * acos(math_clip(qu%w,0.0_pReal,1.0_pReal)) omega = 2.0 * acos(math_clip(qu%w,-1.0_pReal,1.0_pReal))
if (dEq0(omega)) then if (dEq0(omega)) then
ho = [ 0.0, 0.0, 0.0 ] ho = [ 0.0, 0.0, 0.0 ]