mixed up real and aimag part in quaternion

+ some tests
This commit is contained in:
Martin Diehl 2019-09-22 12:10:39 -07:00
parent 25c9bb1cd7
commit 6b5b0fae22
2 changed files with 82 additions and 42 deletions

View File

@ -88,38 +88,52 @@ module quaternions
end type end type
interface assignment (=) interface assignment (=)
module procedure assign_quat__ module procedure assign_quat__
module procedure assign_vec__ module procedure assign_vec__
end interface assignment (=) end interface assignment (=)
interface quaternion interface quaternion
module procedure init__ module procedure init__
end interface quaternion end interface quaternion
interface abs interface abs
procedure abs__ procedure abs__
end interface abs end interface abs
interface dot_product interface dot_product
procedure dot_product__ procedure dot_product__
end interface dot_product end interface dot_product
interface conjg interface conjg
module procedure conjg__ module procedure conjg__
end interface conjg end interface conjg
interface exp interface exp
module procedure exp__ module procedure exp__
end interface exp end interface exp
interface log interface log
module procedure log__ module procedure log__
end interface log end interface log
private :: &
unitTest
contains contains
!--------------------------------------------------------------------------------------------------
!> @brief doing self test
!--------------------------------------------------------------------------------------------------
subroutine quaternions_init()
write(6,'(/,a)') ' <<<+- quaternions init -+>>>'
call unitTest
end subroutine quaternions_init
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
!> constructor for a quaternion from a 4-vector !> constructor for a quaternion from a 4-vector
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
@ -434,29 +448,58 @@ end function asArray
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
!> quaternion as plain array !> real part of a quaternion
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
pure function real__(self) pure function real__(self)
real(pReal), dimension(3) :: real__ real(pReal) :: real__
class(quaternion), intent(in) :: self class(quaternion), intent(in) :: self
real__ = [self%x,self%y,self%z] real__ = self%w
end function real__ end function real__
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
!> quaternion as plain array !> imagninary part of a quaternion
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
pure function aimag__(self) pure function aimag__(self)
real(pReal) :: aimag__ real(pReal), dimension(3) :: aimag__
class(quaternion), intent(in) :: self class(quaternion), intent(in) :: self
aimag__ = self%w aimag__ = [self%x,self%y,self%z]
end function aimag__ end function aimag__
!--------------------------------------------------------------------------------------------------
!> @brief check correctness of (some) quaternions functions
!--------------------------------------------------------------------------------------------------
subroutine unitTest
real(pReal), dimension(4) :: qu
type(quaternion) :: q, q_2
call random_number(qu)
q = qu
write(6,*) q%asArray() == qu
write(6,*) q%real() == qu(1)
write(6,*) q%aimag() == qu(2:4)
q_2 = q%homomorphed()
write(6,*) q == q_2*(-1.0_pReal)
write(6,*) q_2%real() == qu(1)*(-1.0_pReal)
write(6,*) q_2%aimag() == qu(2:4)*(-1.0_pReal)
q_2 = conjg(q)
write(6,*) q_2%real() == q%real()
write(6,*) q_2%aimag() == q%aimag()*(-1.0_pReal)
end subroutine unitTest
end module quaternions end module quaternions

View File

@ -88,10 +88,10 @@ contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief doing self test !> @brief doing self test
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine rotations_init() subroutine rotations_init
write(6,'(/,a)') ' <<<+- rotations init -+>>>' write(6,'(/,a)') ' <<<+- rotations init -+>>>'
call unitTest() call unitTest
end subroutine rotations_init end subroutine rotations_init
@ -265,7 +265,7 @@ pure function rotVector(self,v,active) result(vRot)
else else
q = conjg(self%q) * (quaternion([0.0_pReal, v_normed(1), v_normed(2), v_normed(3)]) * self%q ) q = conjg(self%q) * (quaternion([0.0_pReal, v_normed(1), v_normed(2), v_normed(3)]) * self%q )
endif endif
vRot = q%real()*norm2(v) vRot = q%aimag()*norm2(v)
endif endif
end function rotVector end function rotVector
@ -1233,8 +1233,6 @@ subroutine unitTest
if(qu(1)<0.0_pReal) qu = qu * (-1.0_pReal) if(qu(1)<0.0_pReal) qu = qu * (-1.0_pReal)
endif endif
if(dNeq0(norm2(om2qu(qu2om(qu))-qu),1.0e-12_pReal)) msg = trim(msg)//'om2qu/qu2om,' if(dNeq0(norm2(om2qu(qu2om(qu))-qu),1.0e-12_pReal)) msg = trim(msg)//'om2qu/qu2om,'
if(dNeq0(norm2(eu2qu(qu2eu(qu))-qu),1.0e-12_pReal)) msg = trim(msg)//'eu2qu/qu2eu,' if(dNeq0(norm2(eu2qu(qu2eu(qu))-qu),1.0e-12_pReal)) msg = trim(msg)//'eu2qu/qu2eu,'
if(dNeq0(norm2(ax2qu(qu2ax(qu))-qu),1.0e-12_pReal)) msg = trim(msg)//'ax2qu/qu2ax,' if(dNeq0(norm2(ax2qu(qu2ax(qu))-qu),1.0e-12_pReal)) msg = trim(msg)//'ax2qu/qu2ax,'
@ -1267,11 +1265,10 @@ subroutine unitTest
ho = qu2ho(qu) ho = qu2ho(qu)
if(dNeq0(norm2(ho2qu(cu2ho(ho2cu(ho)))-qu),1.0e-7_pReal)) msg = trim(msg)//'cu2ho/ho2cu,' if(dNeq0(norm2(ho2qu(cu2ho(ho2cu(ho)))-qu),1.0e-7_pReal)) msg = trim(msg)//'cu2ho/ho2cu,'
if(len_trim(msg) /= 0) & if(len_trim(msg) /= 0) call IO_error(401,ext_msg=msg)
call IO_error(401,ext_msg=msg)
enddo enddo
end subroutine unitTest end subroutine unitTest