Quaternion power to a real exponent.

Former-commit-id: e8dc72c04813e592dd1c5a0ccff2971374b3c99c

qpms/qpms_c.pyx

@@ -2,6 +2,7 @@
 # -----------------------------
 
 import numpy as np
+import cmath
 from qpms_cdefs cimport *
 cimport cython
 from cython.parallel cimport parallel, prange
@@ -617,3 +618,112 @@ cdef class trans_calculator:
         return a, b
 
     # TODO make possible to access the attributes (to show normalization etc)
+
+
+
+# Quaternions from wigner.h 
+# (mainly for testing; use moble's quaternions in python)
+
+cdef class cquat:
+    '''
+    Wrapper of the qpms_quat_t object, with the functionality
+    to evaluate Wigner D-matrix elements.
+    '''
+    cdef qpms_quat_t q
+
+    def __cinit__(self, double w, double x, double y, double z):
+        cdef qpms_quat4d_t p
+        p.c1 = w
+        p.ci = x
+        p.cj = y
+        p.ck = z
+        self.q = qpms_quat_2c_from_4d(p)
+
+    def __repr__(self): # TODO make this look like a quaternion with i,j,k
+        return repr(self.r)
+
+    def __add__(cquat self, cquat other):
+        # TODO add real numbers
+        res = cquat(0,0,0,0)
+        res.q = qpms_quat_add(self.q, other.q)
+        return res
+
+    def __mul__(cquat self, cquat other):
+        res = cquat(0,0,0,0)
+        res.q = qpms_quat_mult(self.q, other.q)
+        return res
+
+    def __neg__(cquat self):
+        res = cquat(0,0,0,0)
+        res.q = qpms_quat_rscale(-1, self.q)
+        return res
+
+    def __sub__(cquat self, cquat other):
+        res = cquat(0,0,0,0)
+        res.q = qpms_quat_add(self.q, qpms_quat_rscale(-1,other.q))
+        return res
+
+    def __abs__(self):
+        return qpms_quat_norm(self.q)
+
+    def norm(self):
+        return qpms_quat_norm(self.q)
+
+    def imnorm(self):
+        return qpms_quat_imnorm(self.q)
+
+    def exp(self):
+        res = cquat(0,0,0,0)
+        res.q = qpms_quat_exp(self.q)
+        return res
+
+    def log(self):
+        res = cquat(0,0,0,0)
+        res.q = qpms_quat_exp(self.q)
+        return res
+
+    def __pow__(cquat self, double other, _):
+        res = cquat(0,0,0,0)
+        res.q = qpms_quat_pow(self.q, other)
+        return res
+
+    def normalise(self):
+        res = cquat(0,0,0,0)
+        res.q = qpms_quat_normalise(self.q)
+        return res
+
+    property c:
+        '''
+        Quaternion representation as two complex numbers
+        '''
+        def __get__(self):
+            return (self.q.a, self.q.b)
+        def __set__(self, RaRb):
+            self.q.a = RaRb[0]
+            self.q.b = RaRb[1]
+
+    property r:
+        '''
+        Quaternion representation as four real numbers
+        '''
+        def __get__(self):
+            cdef qpms_quat4d_t p
+            p = qpms_quat_4d_from_2c(self.q)
+            return (p.c1, p.ci, p.cj, p.ck)
+        def __set__(self, wxyz):
+            cdef qpms_quat4d_t p
+            p.c1 = wxyz[0]
+            p.ci = wxyz[1]
+            p.cj = wxyz[2]
+            p.ck = wxyz[3]
+            self.q = qpms_quat_2c_from_4d(p)
+
+    def wignerDelem(self, qpms_l_t l, qpms_m_t mp, qpms_m_t m):
+        '''
+        Returns an element of a bosonic Wigner matrix.
+        '''
+        # don't crash on bad l, m here
+        if (abs(m) > l or abs(mp) > l):
+            return 0
+        return qpms_wignerD_elem(self.q, l, mp, m)
+

qpms/qpms_cdefs.pxd

@@ -34,6 +34,9 @@ cdef extern from "qpms_types.h":
         QPMS_NORMALISATION_SPHARM
         QPMS_NORMALISATION_SPHARM_CS
         QPMS_NORMALISATION_UNDEF
+    ctypedef int qpms_lm_t
+    ctypedef int qpms_l_t
+    ctypedef int qpms_m_t
     # maybe more if needed
 
 # Point generators from lattices.h
@@ -69,9 +72,32 @@ cdef extern from "lattices.h":
     PGen PGen_1D_new_minMaxR(double period, double offset, double minR, bint inc_minR,
             double maxR, bint inc_maxR, PGen_1D_incrementDirection incdir)
 
-
 ctypedef double complex cdouble
 
+cdef extern from "wigner.h":
+    struct qpms_quat_t:
+        cdouble a
+        cdouble b
+    struct qpms_quat4d_t:
+        double c1
+        double ci
+        double cj
+        double ck
+    qpms_quat_t qpms_quat_2c_from_4d(qpms_quat4d_t q)
+    qpms_quat4d_t qpms_quat_4d_from_2c(qpms_quat_t q)
+    qpms_quat_t qpms_quat_mult(qpms_quat_t p, qpms_quat_t q)
+    qpms_quat_t qpms_quat_add(qpms_quat_t p, qpms_quat_t q)
+    qpms_quat_t qpms_quat_rscale(double s, qpms_quat_t q) 
+    qpms_quat_t qpms_quat_conj(qpms_quat_t q) 
+    double qpms_quat_norm(qpms_quat_t q) 
+    double qpms_quat_imnorm(qpms_quat_t q)
+    qpms_quat_t qpms_quat_normalise(qpms_quat_t q) 
+    qpms_quat_t qpms_quat_log(qpms_quat_t q)
+    qpms_quat_t qpms_quat_exp(qpms_quat_t q)
+    qpms_quat_t qpms_quat_pow(qpms_quat_t q, double exponent)
+    cdouble qpms_wignerD_elem(qpms_quat_t q, qpms_l_t l,
+                           qpms_m_t mp, qpms_m_t m)
+
 
 #cdef extern from "numpy/arrayobject.h":
 #    cdef enum NPY_TYPES:

qpms/wigner.h

@@ -28,6 +28,14 @@ static inline qpms_quat_t qpms_quat_2c_from_4d (qpms_quat4d_t q) {
 	return q2c;
 }
 
+/// Conversion from the 2*complex to the 4*double quaternion.
+// TODO is this really correct? 
+// I.e. do the axis from moble's text match this convention?
+static inline qpms_quat4d_t qpms_quat_4d_from_2c (qpms_quat_t q) {
+	qpms_quat4d_t q4d = {creal(q.a), cimag(q.b), creal(q.b), cimag(q.a)};
+	return q4d;
+}
+
 /// Quaternion multiplication.
 /**
  * \f[ (P Q)_a = P_a Q_a - \bar P_b Q_b, \f]
@@ -48,6 +56,18 @@ static inline qpms_quat_t qpms_quat_add(qpms_quat_t p, qpms_quat_t q) {
 	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);
+	const double vn = sqrt(q4.ci*q4.ci + q4.cj*q4.cj + q4.ck *q4.ck);
+	const double ea = exp(q4.c1);
+	const double cv = ea*sin(vn)/vn; // "vector" part common prefactor
+	const qpms_quat4d_t r4 = {ea * cos(vn), cv*q4.ci, cv*q4.cj, cv*q4.ck};
+	return qpms_quat_2c_from_4d(r4);
+}
+
+	
+
 /// Quaternion scaling with a real number.
 static inline qpms_quat_t qpms_quat_rscale(double s, qpms_quat_t q) {
 	qpms_quat_t r = {s * q.a, s * q.b};
@@ -61,22 +81,44 @@ static const qpms_quat_t qpms_quat_j = {0, 1};
 static const qpms_quat_t qpms_quat_k = {0, I};
 
 /// Quaternion conjugation.
-static inline qpms_quat_t qpms_quat_conj(qpms_quat_t q) {
+static inline qpms_quat_t qpms_quat_conj(const qpms_quat_t q) {
 	qpms_quat_t r = {conj(q.a), -q.b};
 	return r;
 }
 
 /// Quaternion norm.
-static inline double qpms_quat_norm(qpms_quat_t q) {
+static inline double qpms_quat_norm(const qpms_quat_t q) {
 	return sqrt(creal(q.a * conj(q.a) + q.b * conj(q.b)));
 }
 
+/// 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);
+	return sqrt(z*z + x*x + y*y);
+}
+
 /// Quaternion normalisation to unit norm.
 static inline qpms_quat_t qpms_quat_normalise(qpms_quat_t q) {
 	double n = qpms_quat_norm(q);
 	return qpms_quat_rscale(1/n, q);
 }
 
+/// Logarithm of a quaternion.
+static inline qpms_quat_t qpms_quat_log(const qpms_quat_t q) {
+	const double n = qpms_quat_norm(q);
+	const double vc = acos(creal(q.a)/n) / qpms_quat_imnorm(q);
+	const qpms_quat_t r = {log(n) + cimag(q.a)*vc*I,
+		q.b*vc};
+	return r;
+}
+
+/// Quaternion power to a real exponent.
+static inline qpms_quat_t qpms_quat_pow(const qpms_quat_t q, const double exponent) {
+	const qpms_quat_t qe = qpms_quat_rscale(exponent, 
+			qpms_quat_log(q));
+	return qpms_quat_exp(qe);
+}
+
 /// Wigner D matrix element from a rotator quaternion for integer l.
 /**
  * The D matrix are calculated using formulae (3), (4), (6), (7) from

setup.py

@@ -36,7 +36,7 @@ qpms_c = Extension('qpms_c',
         )
 
 setup(name='qpms',
-        version = "0.2.993",
+        version = "0.2.994",
         packages=['qpms'],
         setup_requires=['cython>0.28'],
         install_requires=['cython>=0.21','quaternion','spherical_functions',