From 23d2337fb29bc0df81c5dced473135c36e780cec Mon Sep 17 00:00:00 2001
From: Philip Eisenlohr <eisenlohr@egr.msu.edu>
Date: Tue, 15 Nov 2022 16:11:29 -0500
Subject: [PATCH] add option to normalize quaternions

---
 python/damask/_rotation.py    | 13 +++++++++----
 python/tests/test_Rotation.py |  7 ++++---
 2 files changed, 13 insertions(+), 7 deletions(-)

diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py
index b72cec402..e87bcf936 100644
--- a/python/damask/_rotation.py
+++ b/python/damask/_rotation.py
@@ -751,6 +751,7 @@ class Rotation:
     @staticmethod
     def from_quaternion(q: Union[Sequence[FloatSequence], np.ndarray],
                         accept_homomorph: bool = False,
+                        normalize: bool = False,
                         P: Literal[1, -1] = -1) -> 'Rotation':
         """
         Initialize from quaternion.
@@ -762,6 +763,8 @@ class Rotation:
         accept_homomorph : bool, optional
             Allow homomorphic variants, i.e. q_0 < 0 (negative real hemisphere).
             Defaults to False.
+        normalize: bool, optional
+            Allow ǀqǀ ≠ 1. Defaults to False.
         P : int ∈ {-1,1}, optional
             Sign convention. Defaults to -1.
 
@@ -773,12 +776,14 @@ class Rotation:
             raise ValueError('P ∉ {-1,1}')
 
         qu[...,1:4] *= -P
+
         if accept_homomorph:
             qu[qu[...,0] < 0.0] *= -1
-        else:
-            if np.any(qu[...,0] < 0.0):
-                raise ValueError('quaternion with negative first (real) component')
-        if not np.all(np.isclose(np.linalg.norm(qu,axis=-1), 1.0,rtol=0.0)):
+        elif np.any(qu[...,0] < 0.0):
+            raise ValueError('quaternion with negative first (real) component')
+        if normalize:
+            qu /= np.linalg.norm(qu,axis=-1,keepdims=True)
+        elif not np.all(np.isclose(np.linalg.norm(qu,axis=-1), 1.0,rtol=1e-8)):
             raise ValueError('quaternion is not of unit length')
 
         return Rotation(qu)
diff --git a/python/tests/test_Rotation.py b/python/tests/test_Rotation.py
index 8b1bd0de8..d39e20f93 100644
--- a/python/tests/test_Rotation.py
+++ b/python/tests/test_Rotation.py
@@ -760,11 +760,12 @@ class TestRotation:
 
     @pytest.mark.parametrize('P',[1,-1])
     @pytest.mark.parametrize('accept_homomorph',[True,False])
-    def test_quaternion(self,set_of_rotations,P,accept_homomorph):
-        c = np.array([1,P*-1,P*-1,P*-1]) * (-1 if accept_homomorph else 1)
+    @pytest.mark.parametrize('normalize',[True,False])
+    def test_quaternion(self,set_of_rotations,P,accept_homomorph,normalize):
+        c = np.array([1,P*-1,P*-1,P*-1]) * (-1 if accept_homomorph else 1) * (0.9 if normalize else 1.0)
         for rot in set_of_rotations:
             m = rot.as_cubochoric()
-            o = Rotation.from_quaternion(rot.as_quaternion()*c,accept_homomorph,P).as_cubochoric()
+            o = Rotation.from_quaternion(rot.as_quaternion()*c,accept_homomorph,normalize,P).as_cubochoric()
             ok = np.allclose(m,o,atol=atol)
             if np.count_nonzero(np.isclose(np.abs(o),np.pi**(2./3.)*.5)):
                 ok |= np.allclose(m*-1.,o,atol=atol)