From 653837046e0ccfc7e011f4f9a5e975e21f86d6bd Mon Sep 17 00:00:00 2001
From: Claudio	Zambaldi <c.zambaldi@mpie.de>
Date: Mon, 12 Apr 2010 11:07:25 +0000
Subject: [PATCH] new: math_QuaternionToAxisAngle

---
 code/math.f90 | 32 ++++++++++++++++++++++++++++++++
 1 file changed, 32 insertions(+)

diff --git a/code/math.f90 b/code/math.f90
index 59238fdd6..c83f430ec 100644
--- a/code/math.f90
+++ b/code/math.f90
@@ -1413,6 +1413,38 @@ pure function math_transpose3x3(A)
  ENDFUNCTION
 
 
+!********************************************************************
+! axis-angle (x, y, z, ang in deg) from quaternion (w+ix+jy+kz)
+!********************************************************************
+ PURE FUNCTION math_QuaternionToAxisAngle(Q)
+
+ use prec, only: pReal, pInt
+ implicit none
+
+ real(pReal), dimension(4), intent(in) :: Q
+ real(pReal) halfAngle, sinHalfAngle
+ real(pReal), dimension(4) :: math_QuaternionToAxisAngle  
+
+ halfAngle=acos(Q(1))
+ sinHalfAngle=sin(halfAngle)
+
+ math_QuaternionToAxisAngle(1)=Q(2)/sinHalfAngle
+ math_QuaternionToAxisAngle(2)=Q(3)/sinHalfAngle
+ math_QuaternionToAxisAngle(3)=Q(4)/sinHalfAngle
+! Remark: the above calculations gives problems 
+! for HalfAngle->0, i.e. for very small rotation angles
+! and always at inrement 0 where identical orientations 
+! are compared in the calculation of the grainrotation; 
+! the correct interpretation of these special cases 
+! is left to the postprocessing. 
+! A possible integrity check would be to check for 
+! the unit length of the resulting axis.
+
+ math_QuaternionToAxisAngle(4)=halfAngle*2.0_pReal*inDeg
+
+ ENDFUNCTION
+
+
 !****************************************************************
 ! rotation matrix from axis and angle (in radians)  
 !****************************************************************