From 1f8d146b13d1a84fb2317478f36393f54b2af223 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= <marek@necada.org>
Date: Tue, 23 Jul 2019 20:38:35 +0300
Subject: [PATCH] Quaternion fix constants, approximate equality tests etc.

Former-commit-id: 91528cd0141572528a149bbbed074e17a1c19f58
---
 qpms/wigner.h | 65 +++++++++++++++++++++++++++++++++++++++++++++++----
 1 file changed, 60 insertions(+), 5 deletions(-)

diff --git a/qpms/wigner.h b/qpms/wigner.h
index 5328b56..5449c0a 100644
--- a/qpms/wigner.h
+++ b/qpms/wigner.h
@@ -45,6 +45,14 @@ static inline qpms_quat_t qpms_quat_add(qpms_quat_t p, qpms_quat_t q) {
 	return r;
 }
 
+/// Quaternion substraction.
+static inline qpms_quat_t qpms_quat_sub(qpms_quat_t p, qpms_quat_t q) {
+	qpms_quat_t r;
+	r.a = p.a-q.a;
+	r.b = p.b-q.b;
+	return r;
+}
+
 /// Exponential function of a quaternion \f$e^Q$\f.
 static inline qpms_quat_t qpms_quat_exp(const qpms_quat_t q) {
 	const qpms_quat4d_t q4 = qpms_quat_4d_from_2c(q);
@@ -62,10 +70,14 @@ static inline qpms_quat_t qpms_quat_rscale(double s, qpms_quat_t q) {
 }
 
 // quaternion "basis"
-static const qpms_quat_t qpms_quat_1 = {1, 0};
-static const qpms_quat_t qpms_quat_i = {I, 0};
-static const qpms_quat_t qpms_quat_j = {0, 1};
-static const qpms_quat_t qpms_quat_k = {0, I};
+/// Quaternion real unit.
+static const qpms_quat_t QPMS_QUAT_1 = {1, 0};
+/// Quaternion imaginary unit i.
+static const qpms_quat_t QPMS_QUAT_I = {0, I};
+/// Quaternion imaginury unik j.
+static const qpms_quat_t QPMS_QUAT_J = {0, 1};
+/// Quaternion imaginary unit k.
+static const qpms_quat_t QPMS_QUAT_K = {I, 0};
 
 /// Quaternion conjugation.
 static inline qpms_quat_t qpms_quat_conj(const qpms_quat_t q) {
@@ -78,6 +90,11 @@ static inline double qpms_quat_norm(const qpms_quat_t q) {
 	return sqrt(creal(q.a * conj(q.a) + q.b * conj(q.b)));
 }
 
+/// Test approximate equality of quaternions.
+static inline bool qpms_quat_isclose(const qpms_quat_t p, const qpms_quat_t q, double atol) {
+	return qpms_quat_norm(qpms_quat_sub(p,q)) <= atol;
+}
+
 /// Norm of the quaternion imaginary (vector) part.
 static inline double qpms_quat_imnorm(const qpms_quat_t q) {
 	const double z = cimag(q.a), x = cimag(q.b), y = creal(q.b);
@@ -113,7 +130,6 @@ static inline qpms_quat_t qpms_quat_pow(const qpms_quat_t q, const double expone
 	return qpms_quat_exp(qe);
 }
 
-
 /// Quaternion inversion.
 /** \f[ q^{-1} = \frac{q*}{|q|}. \f] */
 static inline qpms_quat_t qpms_quat_inv(const qpms_quat_t q) {
@@ -195,6 +211,11 @@ static inline qpms_irot3_t qpms_irot3_pow(const qpms_irot3_t p, int n) {
 	return r;
 }
 
+/// Test approximate equality of irot3.
+static inline bool qpms_irot3_isclose(const qpms_irot3_t p, const qpms_irot3_t q, double atol) {
+	return qpms_quat_isclose(p.rot, q.rot, atol) && p.det == q.det;
+}
+
 /// Apply an improper rotation onto a 3d cartesian vector.
 static inline cart3_t qpms_irot3_apply_cart3(const qpms_irot3_t p, const cart3_t v) {
 #ifndef NDEBUG
@@ -203,4 +224,38 @@ static inline cart3_t qpms_irot3_apply_cart3(const qpms_irot3_t p, const cart3_t
 	return cart3_scale(p.det, qpms_quat_rot_cart3(p.rot, v));
 }
 
+// Some basic transformations with irot3 type
+/// Identity
+static const qpms_irot3_t QPMS_IROT3_IDENTITY = {QPMS_QUAT_1, 1};
+/// Spatial inversion.
+static const qpms_irot3_t QPMS_IROT3_INVERSION = {QPMS_QUAT_1, -1};
+/// yz-plane mirror symmetry
+static const qpms_irot3_t QPMS_IROT3_XFLIP = {QPMS_QUAT_I, -1};
+/// xz-plane mirror symmetry
+static const qpms_irot3_t QPMS_IROT3_YFLIP = {QPMS_QUAT_J, -1};
+/// xy-plane mirror symmetry
+static const qpms_irot3_t QPMS_IROT3_ZFLIP = {QPMS_QUAT_K, -1};
+
+/// versor representing rotation around z-axis.
+static inline qpms_quat_t qpms_quat_zrot_angle(double angle) {
+	qpms_quat_t q =  {cexp(I*(angle/2)), 0};
+	return q;
+}
+
+/// versor representing rotation \f$ C_N^k \f$, i.e. of angle \f$ 2\pi k / N\f$ around z axis.
+static inline qpms_quat_t qpms_quat_zrot_Nk(double N, double k) {
+	return qpms_quat_zrot_angle(M_PI * k / N);
+}
+
+/// Rotation around z-axis.
+static inline qpms_irot3_t qpms_irot3_zrot_angle(double angle) {
+	qpms_irot3_t q = {qpms_quat_zrot_angle(angle), 1};
+	return q;
+}
+
+/// Rotation \f$ C_N^k \f$, i.e. of angle \f$ 2\pi k / N\f$ around z axis.
+static inline qpms_irot3_t qpms_irot3_zrot_Nk(double N, double k) {
+	return qpms_irot3_zrot_angle(M_PI * k / N);
+}
+
 #endif //QPMS_WIGNER_H