quaternions.h

Go to the documentation of this file.

 
#ifndef QPMS_WIGNER_H
#define QPMS_WIGNER_H
 
#include "qpms_types.h"
#include "vectors.h"
#include "tiny_inlines.h"
 
#define QPMS_QUAT_ATOL (1e-10)
 
// TODO is this really correct? 
// I.e. do the axis from moble's text match this convention?
static inline qpms_quat_t qpms_quat_2c_from_4d (qpms_quat4d_t q) {
    qpms_quat_t q2c = {q.c1 + I * q.ck, q.cj + I * q.ci};
    return q2c;
}
 
// 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;
}
 
 
static inline qpms_quat_t qpms_quat_mult(qpms_quat_t p, qpms_quat_t q) {
    qpms_quat_t r;
    r.a = p.a * q.a - conj(p.b) * q.b;
    r.b = p.b * q.a + conj(p.a) * q.b;
    return r;
}
 
static inline qpms_quat_t qpms_quat_add(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;
}
 
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;
}
 
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 = vn ? (ea*sin(vn)/vn) : ea; // "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);
}
 
static inline qpms_quat_t qpms_quat_rscale(double s, qpms_quat_t q) {
    qpms_quat_t r = {s * q.a, s * q.b};
    return r;
}
 
// quaternion "basis"
static const qpms_quat_t QPMS_QUAT_1 = {1, 0};
static const qpms_quat_t QPMS_QUAT_I = {0, I};
static const qpms_quat_t QPMS_QUAT_J = {0, 1};
static const qpms_quat_t QPMS_QUAT_K = {I, 0};
 
static inline qpms_quat_t qpms_quat_conj(const qpms_quat_t q) {
    qpms_quat_t r = {conj(q.a), -q.b};
    return r;
}
 
static inline double qpms_quat_norm(const qpms_quat_t q) {
    return sqrt(creal(q.a * conj(q.a) + q.b * conj(q.b)));
}
 
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;
}
 
 
static inline qpms_quat_t qpms_quat_standardise(qpms_quat_t p, double atol) {
    //assert(atol >= 0);
    double maxabs = 0;
    int maxi = 0;
    const double *arr = (double *) &(p.a);
    for(int i = 0; i < 4; ++i)
        if (fabs(arr[i]) > maxabs + atol) {
            maxi = i;
            maxabs = fabs(arr[i]);
        }
    if(arr[maxi] < 0) {
        p.a = -p.a;
        p.b = -p.b;
    }
    return p;
}       
 
static inline bool qpms_quat_isclose2(const qpms_quat_t p, const qpms_quat_t q, double atol) {
    return qpms_quat_norm(qpms_quat_sub(
                qpms_quat_standardise(p, atol),
                qpms_quat_standardise(q, atol))) <= atol;
}
 
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);
}
 
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);
}
 
static inline qpms_quat_t qpms_quat_log(const qpms_quat_t q) {
    const double n = qpms_quat_norm(q);
    const double imnorm = qpms_quat_imnorm(q);
    if (imnorm != 0.) {
        const double vc = acos(creal(q.a)/n) / imnorm;
        const qpms_quat_t r = {log(n) + cimag(q.a)*vc*I,
            q.b*vc};
        return r;
    }
    else {
        const qpms_quat_t r = {log(n), 0};
        return r;
    }
}
 
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);
}
 
 
static inline qpms_quat_t qpms_quat_inv(const qpms_quat_t q) {
    const double norm = qpms_quat_norm(q);
    return qpms_quat_rscale(1./(norm*norm),
            qpms_quat_conj(q));
}
            
static inline qpms_quat_t qpms_quat_from_cart3(const cart3_t c) {
    const qpms_quat4d_t q4 = {0, c.x, c.y, c.z};
    return qpms_quat_2c_from_4d(q4);
}
 
static inline cart3_t qpms_quat_to_cart3(const qpms_quat_t q) {
    const qpms_quat4d_t q4 = qpms_quat_4d_from_2c(q);
    const cart3_t c = {q4.ci, q4.cj, q4.ck};
    return c;
}
 
static inline cart3_t qpms_quat_rot_cart3(qpms_quat_t q, const cart3_t v) {
    q = qpms_quat_normalise(q);
    //const qpms_quat_t qc = qpms_quat_normalise(qpms_quat_pow(q, -1)); // implementation of _pow wrong!
    const qpms_quat_t qc = qpms_quat_conj(q);
    const qpms_quat_t vv = qpms_quat_from_cart3(v);
    return qpms_quat_to_cart3(qpms_quat_mult(q,
                qpms_quat_mult(vv, qc)));
}
 
static inline qpms_quat_t qpms_quat_from_rotvector(cart3_t v) {
        return qpms_quat_exp(qpms_quat_rscale(0.5,
                qpms_quat_from_cart3(v)));
}
 
 
complex double qpms_wignerD_elem(qpms_quat_t q, qpms_l_t l,
           qpms_m_t mp, qpms_m_t m);
 
 
complex double qpms_vswf_irot_elem_from_irot3(
        const qpms_irot3_t q, 
        qpms_l_t l, qpms_m_t mp, qpms_m_t m,
        bool pseudo 
        );
 
 
static inline int qpms_irot3_checkdet(const qpms_irot3_t p) {
    if (p.det != 1 && p.det != -1) abort();
    return 0;
}
 
static inline qpms_irot3_t qpms_irot3_mult(const qpms_irot3_t p, const qpms_irot3_t q) {
#ifndef NDEBUG
    qpms_irot3_checkdet(p);
    qpms_irot3_checkdet(q);
#endif
    const qpms_irot3_t r = {qpms_quat_normalise(qpms_quat_mult(p.rot, q.rot)), p.det*q.det};
    return r;
}
 
static inline qpms_irot3_t qpms_irot3_inv(qpms_irot3_t p) {
#ifndef NDEBUG
    qpms_irot3_checkdet(p);
#endif
    p.rot = qpms_quat_inv(p.rot);
    return p;
}
 
static inline qpms_irot3_t qpms_irot3_pow(const qpms_irot3_t p, int n) {
#ifndef NDEBUG
    qpms_irot3_checkdet(p);
#endif
    const qpms_irot3_t r = {qpms_quat_normalise(qpms_quat_pow(p.rot, n)),
        p.det == -1 ? min1pow(n) : 1};
    return r;
}
 
static inline bool qpms_irot3_isclose(const qpms_irot3_t p, const qpms_irot3_t q, double atol) {
    return qpms_quat_isclose2(p.rot, q.rot, atol) && p.det == q.det;
}
 
static inline cart3_t qpms_irot3_apply_cart3(const qpms_irot3_t p, const cart3_t v) {
#ifndef NDEBUG
    qpms_irot3_checkdet(p);
#endif
    return cart3_scale(p.det, qpms_quat_rot_cart3(p.rot, v));
}
 
// Some basic transformations with irot3 type
static const qpms_irot3_t QPMS_IROT3_IDENTITY = {{1, 0}, 1};
static const qpms_irot3_t QPMS_IROT3_XROT_PI =  {{0, I}, 1};
static const qpms_irot3_t QPMS_IROT3_YROT_PI =  {{0, 1}, 1};
static const qpms_irot3_t QPMS_IROT3_ZROT_PI =  {{I, 0}, 1};
static const qpms_irot3_t QPMS_IROT3_INVERSION = {{1, 0}, -1};
static const qpms_irot3_t QPMS_IROT3_XFLIP = {{0, I}, -1};
static const qpms_irot3_t QPMS_IROT3_YFLIP = {{0, 1}, -1};
static const qpms_irot3_t QPMS_IROT3_ZFLIP = {{I, 0}, -1};
 
static inline qpms_quat_t qpms_quat_zrot_angle(double angle) {
    qpms_quat_t q =  {cexp(I*(angle/2)), 0};
    return q;
}
 
static inline qpms_quat_t qpms_quat_zrot_Nk(double N, double k) {
    return qpms_quat_zrot_angle(2 * M_PI * k / N);
}
 
static inline qpms_irot3_t qpms_irot3_zrot_angle(double angle) {
    qpms_irot3_t q = {qpms_quat_zrot_angle(angle), 1};
    return q;
}
 
static inline qpms_irot3_t qpms_irot3_zrot_Nk(double N, double k) {
    return qpms_irot3_zrot_angle(2 * M_PI * k / N);
}
 
#endif //QPMS_WIGNER_H