From 3a08a8bbe27ec3b6b74bbda5516ec97747f5b6ed Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Fri, 10 Jan 2020 12:07:30 -0500 Subject: [PATCH 01/12] always using intrinsic init when assigning quaternions as output variables --- src/quaternions.f90 | 68 +++++++++++++++++---------------------------- 1 file changed, 26 insertions(+), 42 deletions(-) diff --git a/src/quaternions.f90 b/src/quaternions.f90 index 8efb985ed..e55d3804e 100644 --- a/src/quaternions.f90 +++ b/src/quaternions.f90 @@ -142,32 +142,29 @@ type(quaternion) pure function init__(array) real(pReal), intent(in), dimension(4) :: array - init__%w=array(1) - init__%x=array(2) - init__%y=array(3) - init__%z=array(4) + init__%w = array(1) + init__%x = array(2) + init__%y = array(3) + init__%z = array(4) end function init__ !--------------------------------------------------------------------------------------------------- -!> assing a quaternion +!> assigning a quaternion !--------------------------------------------------------------------------------------------------- elemental pure subroutine assign_quat__(self,other) type(quaternion), intent(out) :: self type(quaternion), intent(in) :: other - self%w = other%w - self%x = other%x - self%y = other%y - self%z = other%z - + self = [other%w,other%x,other%y,other%z] + end subroutine assign_quat__ !--------------------------------------------------------------------------------------------------- -!> assing a 4-vector +!> assigning a 4-vector !--------------------------------------------------------------------------------------------------- pure subroutine assign_vec__(self,other) @@ -189,11 +186,8 @@ type(quaternion) elemental pure function add__(self,other) class(quaternion), intent(in) :: self,other - add__%w = self%w + other%w - add__%x = self%x + other%x - add__%y = self%y + other%y - add__%z = self%z + other%z - + add__ = [self%w + other%w,self%x + other%x,self%y + other%y,self%z + other%z] + end function add__ @@ -204,11 +198,8 @@ type(quaternion) elemental pure function pos__(self) class(quaternion), intent(in) :: self - pos__%w = self%w - pos__%x = self%x - pos__%y = self%y - pos__%z = self%z - + pos__ = [self%w,self%x,self%y,self%z] + end function pos__ @@ -219,26 +210,20 @@ type(quaternion) elemental pure function sub__(self,other) class(quaternion), intent(in) :: self,other - sub__%w = self%w - other%w - sub__%x = self%x - other%x - sub__%y = self%y - other%y - sub__%z = self%z - other%z - + sub__ = [self%w - other%w,self%x - other%x,self%y - other%y,self%z - other%z] + end function sub__ !--------------------------------------------------------------------------------------------------- -!> unary positive operator +!> unary negative operator !--------------------------------------------------------------------------------------------------- type(quaternion) elemental pure function neg__(self) class(quaternion), intent(in) :: self - neg__%w = -self%w - neg__%x = -self%x - neg__%y = -self%y - neg__%z = -self%z - + neg__ = [-self%w,-self%x,-self%y,-self%z] + end function neg__ @@ -258,18 +243,15 @@ end function mul_quat__ !--------------------------------------------------------------------------------------------------- -!> multiplication of quaternions with scalar +!> multiplication of quaternion with scalar !--------------------------------------------------------------------------------------------------- type(quaternion) elemental pure function mul_scal__(self,scal) class(quaternion), intent(in) :: self real(pReal), intent(in) :: scal - mul_scal__%w = self%w*scal - mul_scal__%x = self%x*scal - mul_scal__%y = self%y*scal - mul_scal__%z = self%z*scal - + mul_scal__ = [self%w,self%x,self%y,self%z]*scal + end function mul_scal__ @@ -418,7 +400,7 @@ type(quaternion) elemental pure function conjg__(a) class(quaternion), intent(in) :: a - conjg__ = quaternion([a%w, -a%x, -a%y, -a%z]) + conjg__ = [a%w, -a%x, -a%y, -a%z] end function conjg__ @@ -430,7 +412,7 @@ type(quaternion) elemental pure function quat_homomorphed(self) class(quaternion), intent(in) :: self - quat_homomorphed = quaternion(-[self%w,self%x,self%y,self%z]) + quat_homomorphed = -[self%w,self%x,self%y,self%z] end function quat_homomorphed @@ -444,6 +426,7 @@ pure function asArray(self) class(quaternion), intent(in) :: self asArray = [self%w,self%x,self%y,self%z] + if (self%w < 0) asArray = -asArray end function asArray @@ -484,6 +467,7 @@ subroutine unitTest type(quaternion) :: q, q_2 call random_number(qu) + if (qu(1) < 0.0_pReal) qu = -qu q = qu q_2 = q + q @@ -492,10 +476,10 @@ subroutine unitTest q_2 = q - q if(any(dNeq0(q_2%asArray()))) call IO_error(401,ext_msg='sub__') - q_2 = q * 5.0_preal + q_2 = q * 5.0_pReal if(any(dNeq(q_2%asArray(),5.0_pReal*qu))) call IO_error(401,ext_msg='mul__') - q_2 = q / 0.5_preal + q_2 = q / 0.5_pReal if(any(dNeq(q_2%asArray(),2.0_pReal*qu))) call IO_error(401,ext_msg='div__') q_2 = q From 79cafebffe5ae92c029ab8af987c8401ec3ec8f1 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 11 Jan 2020 03:08:39 +0100 Subject: [PATCH 02/12] following https://www.python.org/dev/peps/pep-0257/ --- src/quaternions.f90 | 106 ++++++++++++++++++++++++-------------------- 1 file changed, 57 insertions(+), 49 deletions(-) diff --git a/src/quaternions.f90 b/src/quaternions.f90 index e55d3804e..3df230e31 100644 --- a/src/quaternions.f90 +++ b/src/quaternions.f90 @@ -31,7 +31,8 @@ !> @author Marc De Graef, Carnegie Mellon University !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH !> @brief general quaternion math, not limited to unit quaternions -!> @details w is the real part, (x, y, z) are the imaginary parts. +!> @details w is the real part, (x, y, z) are the imaginary parts. +!> @details https://users.aalto.fi/~ssarkka/pub/quat.pdf !--------------------------------------------------------------------------------------------------- module quaternions use prec @@ -117,6 +118,14 @@ module quaternions interface log module procedure log__ end interface log + + interface real + module procedure real__ + end interface real + + interface aimag + module procedure aimag__ + end interface aimag private :: & unitTest @@ -125,18 +134,18 @@ contains !-------------------------------------------------------------------------------------------------- -!> @brief doing self test +!> @brief do self test !-------------------------------------------------------------------------------------------------- subroutine quaternions_init - write(6,'(/,a)') ' <<<+- quaternions init -+>>>' + write(6,'(/,a)') ' <<<+- quaternions init -+>>>'; flush(6) call unitTest end subroutine quaternions_init !--------------------------------------------------------------------------------------------------- -!> constructor for a quaternion from a 4-vector +!> construct a quaternion from a 4-vector !--------------------------------------------------------------------------------------------------- type(quaternion) pure function init__(array) @@ -151,7 +160,7 @@ end function init__ !--------------------------------------------------------------------------------------------------- -!> assigning a quaternion +!> assign a quaternion !--------------------------------------------------------------------------------------------------- elemental pure subroutine assign_quat__(self,other) @@ -164,7 +173,7 @@ end subroutine assign_quat__ !--------------------------------------------------------------------------------------------------- -!> assigning a 4-vector +!> assign a 4-vector !--------------------------------------------------------------------------------------------------- pure subroutine assign_vec__(self,other) @@ -180,7 +189,7 @@ end subroutine assign_vec__ !--------------------------------------------------------------------------------------------------- -!> addition of two quaternions +!> add a quaternion !--------------------------------------------------------------------------------------------------- type(quaternion) elemental pure function add__(self,other) @@ -192,7 +201,7 @@ end function add__ !--------------------------------------------------------------------------------------------------- -!> unary positive operator +!> return (unary positive operator) !--------------------------------------------------------------------------------------------------- type(quaternion) elemental pure function pos__(self) @@ -204,7 +213,7 @@ end function pos__ !--------------------------------------------------------------------------------------------------- -!> subtraction of two quaternions +!> subtract a quaternion !--------------------------------------------------------------------------------------------------- type(quaternion) elemental pure function sub__(self,other) @@ -216,7 +225,7 @@ end function sub__ !--------------------------------------------------------------------------------------------------- -!> unary negative operator +!> negate (unary negative operator) !--------------------------------------------------------------------------------------------------- type(quaternion) elemental pure function neg__(self) @@ -228,7 +237,7 @@ end function neg__ !--------------------------------------------------------------------------------------------------- -!> multiplication of two quaternions +!> multiply with a quaternion !--------------------------------------------------------------------------------------------------- type(quaternion) elemental pure function mul_quat__(self,other) @@ -243,7 +252,7 @@ end function mul_quat__ !--------------------------------------------------------------------------------------------------- -!> multiplication of quaternion with scalar +!> multiply with a scalar !--------------------------------------------------------------------------------------------------- type(quaternion) elemental pure function mul_scal__(self,scal) @@ -256,7 +265,7 @@ end function mul_scal__ !--------------------------------------------------------------------------------------------------- -!> division of two quaternions +!> divide by a quaternion !--------------------------------------------------------------------------------------------------- type(quaternion) elemental pure function div_quat__(self,other) @@ -268,7 +277,7 @@ end function div_quat__ !--------------------------------------------------------------------------------------------------- -!> divisiont of quaternions by scalar +!> divide by a scalar !--------------------------------------------------------------------------------------------------- type(quaternion) elemental pure function div_scal__(self,scal) @@ -281,7 +290,7 @@ end function div_scal__ !--------------------------------------------------------------------------------------------------- -!> equality of two quaternions +!> test equality !--------------------------------------------------------------------------------------------------- logical elemental pure function eq__(self,other) @@ -294,7 +303,7 @@ end function eq__ !--------------------------------------------------------------------------------------------------- -!> inequality of two quaternions +!> test inequality !--------------------------------------------------------------------------------------------------- logical elemental pure function neq__(self,other) @@ -306,20 +315,7 @@ end function neq__ !--------------------------------------------------------------------------------------------------- -!> quaternion to the power of a scalar -!--------------------------------------------------------------------------------------------------- -type(quaternion) elemental pure function pow_scal__(self,expon) - - class(quaternion), intent(in) :: self - real(pReal), intent(in) :: expon - - pow_scal__ = exp(log(self)*expon) - -end function pow_scal__ - - -!--------------------------------------------------------------------------------------------------- -!> quaternion to the power of a quaternion +!> raise to the power of a quaternion !--------------------------------------------------------------------------------------------------- type(quaternion) elemental pure function pow_quat__(self,expon) @@ -332,7 +328,20 @@ end function pow_quat__ !--------------------------------------------------------------------------------------------------- -!> exponential of a quaternion +!> raise to the power of a scalar +!--------------------------------------------------------------------------------------------------- +type(quaternion) elemental pure function pow_scal__(self,expon) + + class(quaternion), intent(in) :: self + real(pReal), intent(in) :: expon + + pow_scal__ = exp(log(self)*expon) + +end function pow_scal__ + + +!--------------------------------------------------------------------------------------------------- +!> take exponential !> ToDo: Lacks any check for invalid operations !--------------------------------------------------------------------------------------------------- type(quaternion) elemental pure function exp__(self) @@ -340,7 +349,7 @@ type(quaternion) elemental pure function exp__(self) class(quaternion), intent(in) :: self real(pReal) :: absImag - absImag = norm2([self%x, self%y, self%z]) + absImag = norm2(aimag(self)) exp__ = exp(self%w) * [ cos(absImag), & self%x/absImag * sin(absImag), & @@ -351,7 +360,7 @@ end function exp__ !--------------------------------------------------------------------------------------------------- -!> logarithm of a quaternion +!> take logarithm !> ToDo: Lacks any check for invalid operations !--------------------------------------------------------------------------------------------------- type(quaternion) elemental pure function log__(self) @@ -359,7 +368,7 @@ type(quaternion) elemental pure function log__(self) class(quaternion), intent(in) :: self real(pReal) :: absImag - absImag = norm2([self%x, self%y, self%z]) + absImag = norm2(aimag(self)) log__ = [log(abs(self)), & self%x/absImag * acos(self%w/abs(self)), & @@ -370,7 +379,7 @@ end function log__ !--------------------------------------------------------------------------------------------------- -!> norm of a quaternion +!> return norm !--------------------------------------------------------------------------------------------------- real(pReal) elemental pure function abs__(a) @@ -382,7 +391,7 @@ end function abs__ !--------------------------------------------------------------------------------------------------- -!> dot product of two quaternions +!> calculate dot product !--------------------------------------------------------------------------------------------------- real(pReal) elemental pure function dot_product__(a,b) @@ -394,7 +403,7 @@ end function dot_product__ !--------------------------------------------------------------------------------------------------- -!> conjugate complex of a quaternion +!> take conjugate complex !--------------------------------------------------------------------------------------------------- type(quaternion) elemental pure function conjg__(a) @@ -406,7 +415,7 @@ end function conjg__ !--------------------------------------------------------------------------------------------------- -!> homomorphed quaternion of a quaternion +!> homomorph !--------------------------------------------------------------------------------------------------- type(quaternion) elemental pure function quat_homomorphed(self) @@ -418,7 +427,7 @@ end function quat_homomorphed !--------------------------------------------------------------------------------------------------- -!> quaternion as plain array +!> return as plain array !--------------------------------------------------------------------------------------------------- pure function asArray(self) @@ -432,7 +441,7 @@ end function asArray !--------------------------------------------------------------------------------------------------- -!> real part of a quaternion +!> real part (scalar) !--------------------------------------------------------------------------------------------------- pure function real__(self) @@ -445,7 +454,7 @@ end function real__ !--------------------------------------------------------------------------------------------------- -!> imaginary part of a quaternion +!> imaginary part (3-vector) !--------------------------------------------------------------------------------------------------- pure function aimag__(self) @@ -463,37 +472,36 @@ end function aimag__ subroutine unitTest real(pReal), dimension(4) :: qu - type(quaternion) :: q, q_2 call random_number(qu) if (qu(1) < 0.0_pReal) qu = -qu q = qu - + q_2 = q + q if(any(dNeq(q_2%asArray(),2.0_pReal*qu))) call IO_error(401,ext_msg='add__') - + q_2 = q - q if(any(dNeq0(q_2%asArray()))) call IO_error(401,ext_msg='sub__') - + q_2 = q * 5.0_pReal if(any(dNeq(q_2%asArray(),5.0_pReal*qu))) call IO_error(401,ext_msg='mul__') - + q_2 = q / 0.5_pReal if(any(dNeq(q_2%asArray(),2.0_pReal*qu))) call IO_error(401,ext_msg='div__') - + q_2 = q if(q_2 /= q) call IO_error(401,ext_msg='eq__') if(any(dNeq(q%asArray(),qu))) call IO_error(401,ext_msg='eq__') if(dNeq(q%real(), qu(1))) call IO_error(401,ext_msg='real()') if(any(dNeq(q%aimag(), qu(2:4)))) call IO_error(401,ext_msg='aimag()') - + q_2 = q%homomorphed() if(q /= q_2* (-1.0_pReal)) call IO_error(401,ext_msg='homomorphed') if(dNeq(q_2%real(), qu(1)* (-1.0_pReal))) call IO_error(401,ext_msg='homomorphed/real') if(any(dNeq(q_2%aimag(),qu(2:4)*(-1.0_pReal)))) call IO_error(401,ext_msg='homomorphed/aimag') - + q_2 = conjg(q) if(dNeq(q_2%real(), q%real())) call IO_error(401,ext_msg='conjg/real') if(any(dNeq(q_2%aimag(),q%aimag()*(-1.0_pReal)))) call IO_error(401,ext_msg='conjg/aimag') From aefd401e8cdda0ad4164f09728595f4c2467d7b1 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 11 Jan 2020 03:11:45 +0100 Subject: [PATCH 03/12] this is a quaternion class it is meant to represent any quaternion, not only unit quaternions/rotations that follow a specific convention. Need to check in rotations.f90 where the homomorph should happen --- src/quaternions.f90 | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/quaternions.f90 b/src/quaternions.f90 index 3df230e31..9a14fa211 100644 --- a/src/quaternions.f90 +++ b/src/quaternions.f90 @@ -435,7 +435,6 @@ pure function asArray(self) class(quaternion), intent(in) :: self asArray = [self%w,self%x,self%y,self%z] - if (self%w < 0) asArray = -asArray end function asArray @@ -475,7 +474,7 @@ subroutine unitTest type(quaternion) :: q, q_2 call random_number(qu) - if (qu(1) < 0.0_pReal) qu = -qu + qu = (qu-0.5_pReal) * 2.0_pReal q = qu q_2 = q + q From c7180c329571a952b37bca68953910c82588ef79 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 11 Jan 2020 03:50:17 +0100 Subject: [PATCH 04/12] some more tests for quaternion operations --- src/quaternions.f90 | 17 +++++++++++++++-- 1 file changed, 15 insertions(+), 2 deletions(-) diff --git a/src/quaternions.f90 b/src/quaternions.f90 index 9a14fa211..252135eec 100644 --- a/src/quaternions.f90 +++ b/src/quaternions.f90 @@ -475,7 +475,10 @@ subroutine unitTest call random_number(qu) qu = (qu-0.5_pReal) * 2.0_pReal - q = qu + q = quaternion(qu) + + q_2= qu + if(any(dNeq(q%asArray(),q_2%asArray()))) call IO_error(401,ext_msg='assign_vec__') q_2 = q + q if(any(dNeq(q_2%asArray(),2.0_pReal*qu))) call IO_error(401,ext_msg='add__') @@ -489,8 +492,13 @@ subroutine unitTest q_2 = q / 0.5_pReal if(any(dNeq(q_2%asArray(),2.0_pReal*qu))) call IO_error(401,ext_msg='div__') + q_2 = q * 0.3_pReal + if(dNeq0(abs(q)) .and. q_2 == q) call IO_error(401,ext_msg='eq__') + q_2 = q - if(q_2 /= q) call IO_error(401,ext_msg='eq__') + if(q_2 /= q) call IO_error(401,ext_msg='neq__') + + if(dNeq(abs(q),norm2(qu))) call IO_error(401,ext_msg='abs__') if(any(dNeq(q%asArray(),qu))) call IO_error(401,ext_msg='eq__') if(dNeq(q%real(), qu(1))) call IO_error(401,ext_msg='real()') @@ -504,6 +512,11 @@ subroutine unitTest q_2 = conjg(q) if(dNeq(q_2%real(), q%real())) call IO_error(401,ext_msg='conjg/real') if(any(dNeq(q_2%aimag(),q%aimag()*(-1.0_pReal)))) call IO_error(401,ext_msg='conjg/aimag') + + if (norm2(aimag(q)) * abs(real(q)) > 0.0_pReal) then + if (dNeq0(abs(q-exp(log(q))),1.0e-12_pReal)) call IO_error(401,ext_msg='exp/log') + if (dNeq0(abs(q-log(exp(q))),1.0e-12_pReal)) call IO_error(401,ext_msg='log/exp') + endif end subroutine unitTest From 115716b8c21df19868bc5f52c53c0ba2ce38724f Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 11 Jan 2020 03:58:12 +0100 Subject: [PATCH 05/12] polishing/use existing functions --- src/quaternions.f90 | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/src/quaternions.f90 b/src/quaternions.f90 index 252135eec..a49f53007 100644 --- a/src/quaternions.f90 +++ b/src/quaternions.f90 @@ -195,7 +195,8 @@ type(quaternion) elemental pure function add__(self,other) class(quaternion), intent(in) :: self,other - add__ = [self%w + other%w,self%x + other%x,self%y + other%y,self%z + other%z] + add__ = [ self%w, self%x, self%y ,self%z] & + + [other%w, other%x, other%y,other%z] end function add__ @@ -207,7 +208,7 @@ type(quaternion) elemental pure function pos__(self) class(quaternion), intent(in) :: self - pos__ = [self%w,self%x,self%y,self%z] + pos__ = self * (+1.0_pReal) end function pos__ @@ -219,7 +220,8 @@ type(quaternion) elemental pure function sub__(self,other) class(quaternion), intent(in) :: self,other - sub__ = [self%w - other%w,self%x - other%x,self%y - other%y,self%z - other%z] + sub__ = [ self%w, self%x, self%y ,self%z] & + - [other%w, other%x, other%y,other%z] end function sub__ @@ -231,7 +233,7 @@ type(quaternion) elemental pure function neg__(self) class(quaternion), intent(in) :: self - neg__ = [-self%w,-self%x,-self%y,-self%z] + neg__ = self * (-1.0_pReal) end function neg__ @@ -333,7 +335,7 @@ end function pow_quat__ type(quaternion) elemental pure function pow_scal__(self,expon) class(quaternion), intent(in) :: self - real(pReal), intent(in) :: expon + real(pReal), intent(in) :: expon pow_scal__ = exp(log(self)*expon) @@ -421,7 +423,7 @@ type(quaternion) elemental pure function quat_homomorphed(self) class(quaternion), intent(in) :: self - quat_homomorphed = -[self%w,self%x,self%y,self%z] + quat_homomorphed = - self end function quat_homomorphed From de95ca590608c150f081e8fa4d9f247feb188118 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 11 Jan 2020 04:15:51 +0100 Subject: [PATCH 06/12] inverse of a quaternion --- src/quaternions.f90 | 38 ++++++++++++++++++++++++++++---------- 1 file changed, 28 insertions(+), 10 deletions(-) diff --git a/src/quaternions.f90 b/src/quaternions.f90 index a49f53007..0ce7765e6 100644 --- a/src/quaternions.f90 +++ b/src/quaternions.f90 @@ -82,12 +82,13 @@ module quaternions procedure, public :: conjg__ procedure, public :: exp__ procedure, public :: log__ - - procedure, public :: homomorphed => quat_homomorphed - procedure, public :: asArray procedure, public :: real => real__ procedure, public :: aimag => aimag__ + procedure, public :: homomorphed + procedure, public :: asArray + procedure, public :: inverse + end type interface assignment (=) @@ -137,7 +138,6 @@ contains !> @brief do self test !-------------------------------------------------------------------------------------------------- subroutine quaternions_init - write(6,'(/,a)') ' <<<+- quaternions init -+>>>'; flush(6) call unitTest @@ -419,13 +419,13 @@ end function conjg__ !--------------------------------------------------------------------------------------------------- !> homomorph !--------------------------------------------------------------------------------------------------- -type(quaternion) elemental pure function quat_homomorphed(self) +type(quaternion) elemental pure function homomorphed(self) class(quaternion), intent(in) :: self - quat_homomorphed = - self + homomorphed = - self -end function quat_homomorphed +end function homomorphed !--------------------------------------------------------------------------------------------------- @@ -467,6 +467,18 @@ pure function aimag__(self) end function aimag__ +!--------------------------------------------------------------------------------------------------- +!> inverse +!--------------------------------------------------------------------------------------------------- +type(quaternion) elemental pure function inverse(self) + + class(quaternion), intent(in) :: self + + inverse = conjg(self)/abs(self)**2.0_pReal + +end function inverse + + !-------------------------------------------------------------------------------------------------- !> @brief check correctness of (some) quaternions functions !-------------------------------------------------------------------------------------------------- @@ -478,7 +490,7 @@ subroutine unitTest call random_number(qu) qu = (qu-0.5_pReal) * 2.0_pReal q = quaternion(qu) - + q_2= qu if(any(dNeq(q%asArray(),q_2%asArray()))) call IO_error(401,ext_msg='assign_vec__') @@ -515,9 +527,15 @@ subroutine unitTest if(dNeq(q_2%real(), q%real())) call IO_error(401,ext_msg='conjg/real') if(any(dNeq(q_2%aimag(),q%aimag()*(-1.0_pReal)))) call IO_error(401,ext_msg='conjg/aimag') + if(abs(q) > 0.0_pReal) then + q_2 = q * q%inverse() + if( dNeq(real(q_2), 1.0_pReal,1.0e-15_pReal)) call IO_error(401,ext_msg='inverse/real') + if(any(dNeq0(aimag(q_2), 1.0e-15_pReal))) call IO_error(401,ext_msg='inverse/aimag') + endif + if (norm2(aimag(q)) * abs(real(q)) > 0.0_pReal) then - if (dNeq0(abs(q-exp(log(q))),1.0e-12_pReal)) call IO_error(401,ext_msg='exp/log') - if (dNeq0(abs(q-log(exp(q))),1.0e-12_pReal)) call IO_error(401,ext_msg='log/exp') + if (dNeq0(abs(q-exp(log(q))),1.0e-13_pReal)) call IO_error(401,ext_msg='exp/log') + if (dNeq0(abs(q-log(exp(q))),1.0e-13_pReal)) call IO_error(401,ext_msg='log/exp') endif end subroutine unitTest From f02851959765844654bd43265cfd34b1c5cc6a6a Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 11 Jan 2020 04:44:30 +0100 Subject: [PATCH 07/12] some facts from wikipedia as tests --- src/quaternions.f90 | 50 ++++++++++++++++++++++++++------------------- 1 file changed, 29 insertions(+), 21 deletions(-) diff --git a/src/quaternions.f90 b/src/quaternions.f90 index 0ce7765e6..d474f3994 100644 --- a/src/quaternions.f90 +++ b/src/quaternions.f90 @@ -79,9 +79,9 @@ module quaternions procedure, public :: abs__ procedure, public :: dot_product__ - procedure, public :: conjg__ procedure, public :: exp__ procedure, public :: log__ + procedure, public :: conjg => conjg__ procedure, public :: real => real__ procedure, public :: aimag => aimag__ @@ -138,6 +138,7 @@ contains !> @brief do self test !-------------------------------------------------------------------------------------------------- subroutine quaternions_init + write(6,'(/,a)') ' <<<+- quaternions init -+>>>'; flush(6) call unitTest @@ -259,7 +260,7 @@ end function mul_quat__ type(quaternion) elemental pure function mul_scal__(self,scal) class(quaternion), intent(in) :: self - real(pReal), intent(in) :: scal + real(pReal), intent(in) :: scal mul_scal__ = [self%w,self%x,self%y,self%z]*scal @@ -284,7 +285,7 @@ end function div_quat__ type(quaternion) elemental pure function div_scal__(self,scal) class(quaternion), intent(in) :: self - real(pReal), intent(in) :: scal + real(pReal), intent(in) :: scal div_scal__ = [self%w,self%x,self%y,self%z]/scal @@ -492,50 +493,57 @@ subroutine unitTest q = quaternion(qu) q_2= qu - if(any(dNeq(q%asArray(),q_2%asArray()))) call IO_error(401,ext_msg='assign_vec__') + if(any(dNeq(q%asArray(),q_2%asArray()))) call IO_error(401,ext_msg='assign_vec__') q_2 = q + q - if(any(dNeq(q_2%asArray(),2.0_pReal*qu))) call IO_error(401,ext_msg='add__') + if(any(dNeq(q_2%asArray(),2.0_pReal*qu))) call IO_error(401,ext_msg='add__') q_2 = q - q - if(any(dNeq0(q_2%asArray()))) call IO_error(401,ext_msg='sub__') + if(any(dNeq0(q_2%asArray()))) call IO_error(401,ext_msg='sub__') q_2 = q * 5.0_pReal - if(any(dNeq(q_2%asArray(),5.0_pReal*qu))) call IO_error(401,ext_msg='mul__') + if(any(dNeq(q_2%asArray(),5.0_pReal*qu))) call IO_error(401,ext_msg='mul__') q_2 = q / 0.5_pReal - if(any(dNeq(q_2%asArray(),2.0_pReal*qu))) call IO_error(401,ext_msg='div__') + if(any(dNeq(q_2%asArray(),2.0_pReal*qu))) call IO_error(401,ext_msg='div__') q_2 = q * 0.3_pReal - if(dNeq0(abs(q)) .and. q_2 == q) call IO_error(401,ext_msg='eq__') + if(dNeq0(abs(q)) .and. q_2 == q) call IO_error(401,ext_msg='eq__') q_2 = q - if(q_2 /= q) call IO_error(401,ext_msg='neq__') + if(q_2 /= q) call IO_error(401,ext_msg='neq__') - if(dNeq(abs(q),norm2(qu))) call IO_error(401,ext_msg='abs__') + if(dNeq(abs(q),norm2(qu))) call IO_error(401,ext_msg='abs__') + if(dNeq(abs(q)**2.0_pReal, real(q*q%conjg()))) call IO_error(401,ext_msg='abs__/*conjg') - if(any(dNeq(q%asArray(),qu))) call IO_error(401,ext_msg='eq__') - if(dNeq(q%real(), qu(1))) call IO_error(401,ext_msg='real()') - if(any(dNeq(q%aimag(), qu(2:4)))) call IO_error(401,ext_msg='aimag()') + if(any(dNeq(q%asArray(),qu))) call IO_error(401,ext_msg='eq__') + if(dNeq(q%real(), qu(1))) call IO_error(401,ext_msg='real()') + if(any(dNeq(q%aimag(), qu(2:4)))) call IO_error(401,ext_msg='aimag()') q_2 = q%homomorphed() - if(q /= q_2* (-1.0_pReal)) call IO_error(401,ext_msg='homomorphed') - if(dNeq(q_2%real(), qu(1)* (-1.0_pReal))) call IO_error(401,ext_msg='homomorphed/real') - if(any(dNeq(q_2%aimag(),qu(2:4)*(-1.0_pReal)))) call IO_error(401,ext_msg='homomorphed/aimag') + if(q /= q_2* (-1.0_pReal)) call IO_error(401,ext_msg='homomorphed') + if(dNeq(q_2%real(), qu(1)* (-1.0_pReal))) call IO_error(401,ext_msg='homomorphed/real') + if(any(dNeq(q_2%aimag(),qu(2:4)*(-1.0_pReal)))) call IO_error(401,ext_msg='homomorphed/aimag') q_2 = conjg(q) - if(dNeq(q_2%real(), q%real())) call IO_error(401,ext_msg='conjg/real') - if(any(dNeq(q_2%aimag(),q%aimag()*(-1.0_pReal)))) call IO_error(401,ext_msg='conjg/aimag') + if(dNeq(abs(q),abs(q_2))) call IO_error(401,ext_msg='conjg/abs') + if(q /= conjg(q_2)) call IO_error(401,ext_msg='conjg/involution') + if(dNeq(q_2%real(), q%real())) call IO_error(401,ext_msg='conjg/real') + if(any(dNeq(q_2%aimag(),q%aimag()*(-1.0_pReal)))) call IO_error(401,ext_msg='conjg/aimag') if(abs(q) > 0.0_pReal) then q_2 = q * q%inverse() if( dNeq(real(q_2), 1.0_pReal,1.0e-15_pReal)) call IO_error(401,ext_msg='inverse/real') if(any(dNeq0(aimag(q_2), 1.0e-15_pReal))) call IO_error(401,ext_msg='inverse/aimag') + + q_2 = q/abs(q) + q_2 = conjg(q_2) - inverse(q_2) + if(any(dNeq0(q_2%asArray(),1.0e-15_pReal))) call IO_error(401,ext_msg='inverse/conjg') endif if (norm2(aimag(q)) * abs(real(q)) > 0.0_pReal) then - if (dNeq0(abs(q-exp(log(q))),1.0e-13_pReal)) call IO_error(401,ext_msg='exp/log') - if (dNeq0(abs(q-log(exp(q))),1.0e-13_pReal)) call IO_error(401,ext_msg='log/exp') + if (dNeq0(abs(q-exp(log(q))),1.0e-13_pReal)) call IO_error(401,ext_msg='exp/log') + if (dNeq0(abs(q-log(exp(q))),1.0e-13_pReal)) call IO_error(401,ext_msg='log/exp') endif end subroutine unitTest From 3a6819f548897d20a04a1df5498b764ee79d5413 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 11 Jan 2020 05:14:17 +0100 Subject: [PATCH 08/12] check for invalid operations --- src/quaternions.f90 | 24 +++++++++++++----------- 1 file changed, 13 insertions(+), 11 deletions(-) diff --git a/src/quaternions.f90 b/src/quaternions.f90 index d474f3994..bdd44b93b 100644 --- a/src/quaternions.f90 +++ b/src/quaternions.f90 @@ -345,7 +345,6 @@ end function pow_scal__ !--------------------------------------------------------------------------------------------------- !> take exponential -!> ToDo: Lacks any check for invalid operations !--------------------------------------------------------------------------------------------------- type(quaternion) elemental pure function exp__(self) @@ -354,17 +353,18 @@ type(quaternion) elemental pure function exp__(self) absImag = norm2(aimag(self)) - exp__ = exp(self%w) * [ cos(absImag), & - self%x/absImag * sin(absImag), & - self%y/absImag * sin(absImag), & - self%z/absImag * sin(absImag)] + exp__ = merge(exp(self%w) * [ cos(absImag), & + self%x/absImag * sin(absImag), & + self%y/absImag * sin(absImag), & + self%z/absImag * sin(absImag)], & + IEEE_value(1.0_pReal,IEEE_SIGNALING_NAN), & + dNeq0(absImag)) end function exp__ !--------------------------------------------------------------------------------------------------- !> take logarithm -!> ToDo: Lacks any check for invalid operations !--------------------------------------------------------------------------------------------------- type(quaternion) elemental pure function log__(self) @@ -373,10 +373,12 @@ type(quaternion) elemental pure function log__(self) absImag = norm2(aimag(self)) - log__ = [log(abs(self)), & - self%x/absImag * acos(self%w/abs(self)), & - self%y/absImag * acos(self%w/abs(self)), & - self%z/absImag * acos(self%w/abs(self))] + log__ = merge([log(abs(self)), & + self%x/absImag * acos(self%w/abs(self)), & + self%y/absImag * acos(self%w/abs(self)), & + self%z/absImag * acos(self%w/abs(self))], & + IEEE_value(1.0_pReal,IEEE_SIGNALING_NAN), & + dNeq0(absImag)) end function log__ @@ -541,7 +543,7 @@ subroutine unitTest if(any(dNeq0(q_2%asArray(),1.0e-15_pReal))) call IO_error(401,ext_msg='inverse/conjg') endif - if (norm2(aimag(q)) * abs(real(q)) > 0.0_pReal) then + if (norm2(aimag(q)) > 0.0_pReal) then if (dNeq0(abs(q-exp(log(q))),1.0e-13_pReal)) call IO_error(401,ext_msg='exp/log') if (dNeq0(abs(q-log(exp(q))),1.0e-13_pReal)) call IO_error(401,ext_msg='log/exp') endif From 842666cc20a4aafd1608f22e70ab1b152c1f9396 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 11 Jan 2020 11:23:53 +0100 Subject: [PATCH 09/12] no overlap with Marc's code --- src/quaternions.f90 | 33 ++------------------------------- 1 file changed, 2 insertions(+), 31 deletions(-) diff --git a/src/quaternions.f90 b/src/quaternions.f90 index bdd44b93b..886c7d071 100644 --- a/src/quaternions.f90 +++ b/src/quaternions.f90 @@ -1,38 +1,9 @@ -! ################################################################### -! Copyright (c) 2013-2015, Marc De Graef/Carnegie Mellon University -! Modified 2017-2019, Martin Diehl/Max-Planck-Institut für Eisenforschung GmbH -! All rights reserved. -! -! Redistribution and use in source and binary forms, with or without modification, are -! permitted provided that the following conditions are met: -! -! - Redistributions of source code must retain the above copyright notice, this list -! of conditions and the following disclaimer. -! - Redistributions in binary form must reproduce the above copyright notice, this -! list of conditions and the following disclaimer in the documentation and/or -! other materials provided with the distribution. -! - Neither the names of Marc De Graef, Carnegie Mellon University nor the names -! of its contributors may be used to endorse or promote products derived from -! this software without specific prior written permission. -! -! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" -! AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE -! IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE -! ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE -! LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL -! DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR -! SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER -! CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, -! OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE -! USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. -! ################################################################### - !--------------------------------------------------------------------------------------------------- -!> @author Marc De Graef, Carnegie Mellon University !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH +!> @author Philip Eisenlohr, Michigan State University !> @brief general quaternion math, not limited to unit quaternions !> @details w is the real part, (x, y, z) are the imaginary parts. -!> @details https://users.aalto.fi/~ssarkka/pub/quat.pdf +!> @details https://en.wikipedia.org/wiki/Quaternion !--------------------------------------------------------------------------------------------------- module quaternions use prec From e762cb4dfd5352ac18b2d7f53c926dacac1ad779 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 11 Jan 2020 12:36:35 +0100 Subject: [PATCH 10/12] issue with gfortran < 9 the false branch of merge seems to be evaluated which results in a signaling NaN --- src/quaternions.f90 | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/quaternions.f90 b/src/quaternions.f90 index 886c7d071..a67f345a0 100644 --- a/src/quaternions.f90 +++ b/src/quaternions.f90 @@ -513,12 +513,14 @@ subroutine unitTest q_2 = conjg(q_2) - inverse(q_2) if(any(dNeq0(q_2%asArray(),1.0e-15_pReal))) call IO_error(401,ext_msg='inverse/conjg') endif - + +#if !(defined(__GFORTRAN__) && __GNUC__ < 9) if (norm2(aimag(q)) > 0.0_pReal) then if (dNeq0(abs(q-exp(log(q))),1.0e-13_pReal)) call IO_error(401,ext_msg='exp/log') if (dNeq0(abs(q-log(exp(q))),1.0e-13_pReal)) call IO_error(401,ext_msg='log/exp') endif - +#endif + end subroutine unitTest From ac112d2d36f783b942186c5d6fa003bb62f194d3 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 11 Jan 2020 13:55:56 +0100 Subject: [PATCH 11/12] tolerance needed for optimized code --- src/quaternions.f90 | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/quaternions.f90 b/src/quaternions.f90 index a67f345a0..0ca404880 100644 --- a/src/quaternions.f90 +++ b/src/quaternions.f90 @@ -487,7 +487,8 @@ subroutine unitTest if(q_2 /= q) call IO_error(401,ext_msg='neq__') if(dNeq(abs(q),norm2(qu))) call IO_error(401,ext_msg='abs__') - if(dNeq(abs(q)**2.0_pReal, real(q*q%conjg()))) call IO_error(401,ext_msg='abs__/*conjg') + if(dNeq(abs(q)**2.0_pReal, real(q*q%conjg()),1.0e-14_pReal)) & + call IO_error(401,ext_msg='abs__/*conjg') if(any(dNeq(q%asArray(),qu))) call IO_error(401,ext_msg='eq__') if(dNeq(q%real(), qu(1))) call IO_error(401,ext_msg='real()') From 300f1b7015cf644e4c175aea43483d651493bff8 Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Sat, 11 Jan 2020 11:36:22 -0500 Subject: [PATCH 12/12] added options to return "natural" versions of asQ, asRodrig, and asAxisAngle --- python/damask/orientation.py | 40 +++++++++++++++++++++++------------- 1 file changed, 26 insertions(+), 14 deletions(-) diff --git a/python/damask/orientation.py b/python/damask/orientation.py index a63444155..a86ba9331 100644 --- a/python/damask/orientation.py +++ b/python/damask/orientation.py @@ -170,9 +170,18 @@ class Rotation: ################################################################################################ # convert to different orientation representations (numpy arrays) - def asQuaternion(self): - """Unit quaternion: (q, p_1, p_2, p_3).""" - return self.quaternion.asArray() + def asQuaternion(self, + quaternion = False): + """ + Unit quaternion [q, p_1, p_2, p_3] unless quaternion == True: damask.quaternion object. + + Parameters + ---------- + quaternion : bool, optional + return quaternion as DAMASK object. + + """ + return self.quaternion if quaternion else self.quaternion.asArray() def asEulers(self, degrees = False): @@ -190,33 +199,36 @@ class Rotation: return eu def asAxisAngle(self, - degrees = False): + degrees = False, + pair = False): """ - Axis angle pair: ([n_1, n_2, n_3], ω). + Axis angle representation [n_1, n_2, n_3, ω] unless pair == True: ([n_1, n_2, n_3], ω). Parameters ---------- degrees : bool, optional return rotation angle in degrees. + pair : bool, optional + return tuple of axis and angle. """ ax = qu2ax(self.quaternion.asArray()) if degrees: ax[3] = np.degrees(ax[3]) - return ax + return (ax[:3],np.degrees(ax[3])) if pair else ax def asMatrix(self): """Rotation matrix.""" return qu2om(self.quaternion.asArray()) def asRodrigues(self, - vector=False): + vector = False): """ - Rodrigues-Frank vector: ([n_1, n_2, n_3], tan(ω/2)). + Rodrigues-Frank vector representation [n_1, n_2, n_3, tan(ω/2)] unless vector == True: [n_1, n_2, n_3] * tan(ω/2). Parameters ---------- vector : bool, optional - return as array of length 3, i.e. scale the unit vector giving the rotation axis. + return as actual Rodrigues--Frank vector, i.e. rotation axis scaled by tan(ω/2). """ ro = qu2ro(self.quaternion.asArray()) @@ -252,8 +264,8 @@ class Rotation: acceptHomomorph = False, P = -1): - qu = quaternion if isinstance(quaternion, np.ndarray) and quaternion.dtype == np.dtype(float) \ - else np.array(quaternion,dtype=float) + qu = quaternion if isinstance(quaternion, np.ndarray) and quaternion.dtype == np.dtype(float) \ + else np.array(quaternion,dtype=float) if P > 0: qu[1:4] *= -1 # convert from P=1 to P=-1 if qu[0] < 0.0: if acceptHomomorph: @@ -1193,9 +1205,9 @@ class Orientation: ref = orientations[0] for o in orientations: closest.append(o.equivalentOrientations( - ref.disorientation(o, - SST = False, # select (o[ther]'s) sym orientation - symmetries = True)[2]).rotation) # with lowest misorientation + ref.disorientation(o, + SST = False, # select (o[ther]'s) sym orientation + symmetries = True)[2]).rotation) # with lowest misorientation return Orientation(Rotation.fromAverage(closest,weights),ref.lattice)