From ed50cd022b3d922e89e4d1c080c7b1d26ef2a129 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 16 Feb 2022 23:23:03 +0100 Subject: [PATCH] shorter, potential for higher precision np.sum has an better alogrithm but fails ... --- python/damask/_rotation.py | 6 +----- python/tests/test_Rotation.py | 8 ++++---- 2 files changed, 5 insertions(+), 9 deletions(-) diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index 9c6bde2dc..95c0b385a 100644 --- a/python/damask/_rotation.py +++ b/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 ], diff --git a/python/tests/test_Rotation.py b/python/tests/test_Rotation.py index fda986c2f..57f078ae5 100644 --- a/python/tests/test_Rotation.py +++ b/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