shorter, potential for higher precision

np.sum has an better alogrithm but fails ...

python/damask/_rotation.py

@@ -1567,11 +1567,7 @@ class Rotation:
                          +0.000059719705868660826, -0.00001980756723965647,
                          +0.000003953714684212874, -0.00000036555001439719544])
         hmag_squared = np.sum(ho**2.,axis=-1,keepdims=True)
-        hm = hmag_squared.copy()
-        s = tfit[0] + tfit[1] * hmag_squared
-        for i in range(2,16):
-            hm *= hmag_squared
-            s  += tfit[i] * hm
+        s = sum([t*hmag_squared**i for i,t in enumerate(tfit)])                                     # np.sum fails due to higher precision
         with np.errstate(invalid='ignore'):
             ax = np.where(np.broadcast_to(np.abs(hmag_squared)<1.e-8,ho.shape[:-1]+(4,)),
                           [ 0.0, 0.0, 1.0, 0.0 ],

python/tests/test_Rotation.py

@@ -313,13 +313,13 @@ def ho2ax(ho):
     if iszero(hmag_squared):
         ax = np.array([ 0.0, 0.0, 1.0, 0.0 ])
     else:
-        hm = hmag_squared
+        hm = np.ones_like(hmag_squared)
 
         # convert the magnitude to the rotation angle
-        s = tfit[0] + tfit[1] * hmag_squared
-        for i in range(2,16):
+        s = 0.0
+        for t in tfit:
+            s += t * hm
             hm *= hmag_squared
-            s  += tfit[i] * hm
         ax = np.append(ho/np.sqrt(hmag_squared),2.0*np.arccos(np.clip(s,-1.0,1.0)))
     return ax