From 6b5b0fae22d9ed2c5ad489d0dd651dc65c923129 Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Sun, 22 Sep 2019 12:10:39 -0700
Subject: [PATCH] mixed up real and aimag part in quaternion

+ some tests
---
 src/quaternions.f90 | 111 ++++++++++++++++++++++++++++++--------------
 src/rotations.f90   |  13 ++----
 2 files changed, 82 insertions(+), 42 deletions(-)

diff --git a/src/quaternions.f90 b/src/quaternions.f90
index 37363ad12..314c88466 100644
--- a/src/quaternions.f90
+++ b/src/quaternions.f90
@@ -88,38 +88,52 @@ module quaternions
 
   end type
 
-interface assignment (=)
-  module procedure assign_quat__
-  module procedure assign_vec__
-end interface assignment (=)
-
-interface quaternion
-  module procedure init__
-end interface quaternion
-
-interface abs
-  procedure abs__
-end interface abs
-
-interface dot_product
-  procedure dot_product__
-end interface dot_product
-
-interface conjg
-  module procedure conjg__
-end interface conjg
-
-interface exp
-  module procedure exp__
-end interface exp
-
-interface log
-  module procedure log__
-end interface log
+  interface assignment (=)
+    module procedure assign_quat__
+    module procedure assign_vec__
+  end interface assignment (=)
+  
+  interface quaternion
+    module procedure init__
+  end interface quaternion
+  
+  interface abs
+    procedure abs__
+  end interface abs
+  
+  interface dot_product
+    procedure dot_product__
+  end interface dot_product
+  
+  interface conjg
+    module procedure conjg__
+  end interface conjg
+  
+  interface exp
+    module procedure exp__
+  end interface exp
+  
+  interface log
+    module procedure log__
+  end interface log
+  
+  private :: &
+    unitTest
 
 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
 !---------------------------------------------------------------------------------------------------
@@ -434,29 +448,58 @@ end function asArray
 
 
 !---------------------------------------------------------------------------------------------------
-!> quaternion as plain array
+!> real part of a quaternion
 !---------------------------------------------------------------------------------------------------
 pure function real__(self)
 
-  real(pReal), dimension(3)     :: real__
+  real(pReal)                   :: real__
   class(quaternion), intent(in) :: self
 
-  real__ = [self%x,self%y,self%z]
+  real__ = self%w
 
 end function real__
 
 
 !---------------------------------------------------------------------------------------------------
-!> quaternion as plain array
+!> imagninary part of a quaternion
 !---------------------------------------------------------------------------------------------------
 pure function aimag__(self)
 
-  real(pReal)                   :: aimag__
+  real(pReal), dimension(3)     :: aimag__
   class(quaternion), intent(in) :: self
 
-  aimag__ = self%w
+  aimag__ = [self%x,self%y,self%z]
 
 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
diff --git a/src/rotations.f90 b/src/rotations.f90
index c3d9a233e..bbbd51ee5 100644
--- a/src/rotations.f90
+++ b/src/rotations.f90
@@ -88,10 +88,10 @@ contains
 !--------------------------------------------------------------------------------------------------
 !> @brief doing self test
 !--------------------------------------------------------------------------------------------------
-subroutine rotations_init()
+subroutine rotations_init
 
   write(6,'(/,a)')   ' <<<+-  rotations init  -+>>>'
-  call unitTest()
+  call unitTest
 
 end subroutine rotations_init
 
@@ -265,7 +265,7 @@ pure function rotVector(self,v,active) result(vRot)
     else
       q = conjg(self%q) * (quaternion([0.0_pReal, v_normed(1), v_normed(2), v_normed(3)]) * self%q )
     endif
-    vRot = q%real()*norm2(v)
+    vRot = q%aimag()*norm2(v)
   endif
 
 end function rotVector
@@ -1233,8 +1233,6 @@ subroutine unitTest
       if(qu(1)<0.0_pReal) qu = qu * (-1.0_pReal)
     endif
     
- 
-    
     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(ax2qu(qu2ax(qu))-qu),1.0e-12_pReal)) msg = trim(msg)//'ax2qu/qu2ax,'
@@ -1267,11 +1265,10 @@ subroutine unitTest
     ho = qu2ho(qu)
     if(dNeq0(norm2(ho2qu(cu2ho(ho2cu(ho)))-qu),1.0e-7_pReal))  msg = trim(msg)//'cu2ho/ho2cu,'
 
-    if(len_trim(msg) /= 0) &
-      call IO_error(401,ext_msg=msg)
+    if(len_trim(msg) /= 0) call IO_error(401,ext_msg=msg)
 
   enddo
-      
+
 end subroutine unitTest