qpms/quaternions.h

   1 /*! \file quaternions.h
   2  * \brief Quaternions and Wigner matrices
   3  */
   4 #ifndef QPMS_WIGNER_H
   5 #define QPMS_WIGNER_H
   6
   7 #include "qpms_types.h"
   8 #include "vectors.h"
   9 #include "tiny_inlines.h"
  10
  11 /// Just some arbitrarily chosen "default" value for quaternion comparison tolerance.
  12 #define QPMS_QUAT_ATOL (1e-10)
  13
  14 /// Conversion from the 4*double to the 2*complex quaternion.
  15 // TODO is this really correct?
  16 // I.e. do the axis from moble's text match this convention?
  17 static inline qpms_quat_t qpms_quat_2c_from_4d (qpms_quat4d_t q) {
  18         qpms_quat_t q2c = {q.c1 + I * q.ck, q.cj + I * q.ci};
  19         return q2c;
  20 }
  21
  22 /// Conversion from the 2*complex to the 4*double quaternion.
  23 // TODO is this really correct?
  24 // I.e. do the axis from moble's text match this convention?
  25 static inline qpms_quat4d_t qpms_quat_4d_from_2c (qpms_quat_t q) {
  26         qpms_quat4d_t q4d = {creal(q.a), cimag(q.b), creal(q.b), cimag(q.a)};
  27         return q4d;
  28 }
  29
  30 /// Quaternion multiplication.
  31 /**
  32  * \f[ (P Q)_a = P_a Q_a - \bar P_b Q_b, \f]
  33  * \f[ (P Q)_b = P_b Q_a + \bar P_a Q_b. \f]
  34  */
  35 static inline qpms_quat_t qpms_quat_mult(qpms_quat_t p, qpms_quat_t q) {
  36         qpms_quat_t r;
  37         r.a = p.a * q.a - conj(p.b) * q.b;
  38         r.b = p.b * q.a + conj(p.a) * q.b;
  39         return r;
  40 }
  41
  42 /// Quaternion addition.
  43 static inline qpms_quat_t qpms_quat_add(qpms_quat_t p, qpms_quat_t q) {
  44         qpms_quat_t r;
  45         r.a = p.a+q.a;
  46         r.b = p.b+q.b;
  47         return r;
  48 }
  49
  50 /// Quaternion substraction.
  51 static inline qpms_quat_t qpms_quat_sub(qpms_quat_t p, qpms_quat_t q) {
  52         qpms_quat_t r;
  53         r.a = p.a-q.a;
  54         r.b = p.b-q.b;
  55         return r;
  56 }
  57
  58 /// Exponential function of a quaternion \f$e^Q$\f.
  59 static inline qpms_quat_t qpms_quat_exp(const qpms_quat_t q) {
  60         const qpms_quat4d_t q4 = qpms_quat_4d_from_2c(q);
  61         const double vn = sqrt(q4.ci*q4.ci + q4.cj*q4.cj + q4.ck *q4.ck);
  62         const double ea = exp(q4.c1);
  63         const double cv = vn ? (ea*sin(vn)/vn) : ea; // "vector" part common prefactor
  64         const qpms_quat4d_t r4 = {ea * cos(vn), cv*q4.ci, cv*q4.cj, cv*q4.ck};
  65         return qpms_quat_2c_from_4d(r4);
  66 }
  67
  68 /// Quaternion scaling with a real number.
  69 static inline qpms_quat_t qpms_quat_rscale(double s, qpms_quat_t q) {
  70         qpms_quat_t r = {s * q.a, s * q.b};
  71         return r;
  72 }
  73
  74 // quaternion "basis"
  75 /// Quaternion real unit.
  76 static const qpms_quat_t QPMS_QUAT_1 = {1, 0};
  77 /// Quaternion imaginary unit i.
  78 static const qpms_quat_t QPMS_QUAT_I = {0, I};
  79 /// Quaternion imaginury unik j.
  80 static const qpms_quat_t QPMS_QUAT_J = {0, 1};
  81 /// Quaternion imaginary unit k.
  82 static const qpms_quat_t QPMS_QUAT_K = {I, 0};
  83
  84 /// Quaternion conjugation.
  85 static inline qpms_quat_t qpms_quat_conj(const qpms_quat_t q) {
  86         qpms_quat_t r = {conj(q.a), -q.b};
  87         return r;
  88 }
  89
  90 /// Quaternion norm.
  91 static inline double qpms_quat_norm(const qpms_quat_t q) {
  92         return sqrt(creal(q.a * conj(q.a) + q.b * conj(q.b)));
  93 }
  94
  95 /// Test approximate equality of quaternions.
  96 static inline bool qpms_quat_isclose(const qpms_quat_t p, const qpms_quat_t q, double atol) {
  97         return qpms_quat_norm(qpms_quat_sub(p,q)) <= atol;
  98 }
  99
 100 /// "Standardises" a quaternion to have the largest component "positive".
 101 /**
 102  * This is to remove the ambiguity stemming from the double cover of SO(3).
 103  */
 104 static inline qpms_quat_t qpms_quat_standardise(qpms_quat_t p, double atol) {
 105         //assert(atol >= 0);
 106         double maxabs = 0;
 107         int maxi = 0;
 108         const double *arr = (double *) &(p.a);
 109         for(int i = 0; i < 4; ++i)
 110                 if (fabs(arr[i]) > maxabs + atol) {
 111                         maxi = i;
 112                         maxabs = fabs(arr[i]);
 113                 }
 114         if(arr[maxi] < 0) {
 115                 p.a = -p.a;
 116                 p.b = -p.b;
 117         }
 118         return p;
 119 }
 120
 121 /// Test approximate equality of "standardised" quaternions, so that \f$-q\f$ is considered equal to \f$q\f$.
 122 static inline bool qpms_quat_isclose2(const qpms_quat_t p, const qpms_quat_t q, double atol) {
 123         return qpms_quat_norm(qpms_quat_sub(
 124                                 qpms_quat_standardise(p, atol),
 125                                 qpms_quat_standardise(q, atol))) <= atol;
 126 }
 127
 128 /// Norm of the quaternion imaginary (vector) part.
 129 static inline double qpms_quat_imnorm(const qpms_quat_t q) {
 130         const double z = cimag(q.a), x = cimag(q.b), y = creal(q.b);
 131         return sqrt(z*z + x*x + y*y);
 132 }
 133
 134 /// Quaternion normalisation to unit norm.
 135 static inline qpms_quat_t qpms_quat_normalise(qpms_quat_t q) {
 136         double n = qpms_quat_norm(q);
 137         return qpms_quat_rscale(1/n, q);
 138 }
 139
 140 /// Logarithm of a quaternion.
 141 static inline qpms_quat_t qpms_quat_log(const qpms_quat_t q) {
 142         const double n = qpms_quat_norm(q);
 143         const double imnorm = qpms_quat_imnorm(q);
 144         if (imnorm != 0.) {
 145                 const double vc = acos(creal(q.a)/n) / imnorm;
 146                 const qpms_quat_t r = {log(n) + cimag(q.a)*vc*I,
 147                         q.b*vc};
 148                 return r;
 149         }
 150         else {
 151                 const qpms_quat_t r = {log(n), 0};
 152                 return r;
 153         }
 154 }
 155
 156 /// Quaternion power to a real exponent.
 157 static inline qpms_quat_t qpms_quat_pow(const qpms_quat_t q, const double exponent) {
 158         const qpms_quat_t qe = qpms_quat_rscale(exponent,
 159                         qpms_quat_log(q));
 160         return qpms_quat_exp(qe);
 161 }
 162
 163 /// Quaternion inversion.
 164 /** \f[ q^{-1} = \frac{q*}{|q|}. \f] */
 165 static inline qpms_quat_t qpms_quat_inv(const qpms_quat_t q) {
 166         const double norm = qpms_quat_norm(q);
 167         return qpms_quat_rscale(1./(norm*norm),
 168                         qpms_quat_conj(q));
 169 }
 170
 171 /// Make a pure imaginary quaternion from a 3d cartesian vector.
 172 static inline qpms_quat_t qpms_quat_from_cart3(const cart3_t c) {
 173         const qpms_quat4d_t q4 = {0, c.x, c.y, c.z};
 174         return qpms_quat_2c_from_4d(q4);
 175 }
 176
 177 /// Make a 3d cartesian vector from the imaginary part of a quaternion.
 178 static inline cart3_t qpms_quat_to_cart3(const qpms_quat_t q) {
 179         const qpms_quat4d_t q4 = qpms_quat_4d_from_2c(q);
 180         const cart3_t c = {q4.ci, q4.cj, q4.ck};
 181         return c;
 182 }
 183
 184 /// Rotate a 3-dimensional cartesian vector using the quaternion/versor representation.
 185 static inline cart3_t qpms_quat_rot_cart3(qpms_quat_t q, const cart3_t v) {
 186         q = qpms_quat_normalise(q);
 187         //const qpms_quat_t qc = qpms_quat_normalise(qpms_quat_pow(q, -1)); // implementation of _pow wrong!
 188         const qpms_quat_t qc = qpms_quat_conj(q);
 189         const qpms_quat_t vv = qpms_quat_from_cart3(v);
 190         return qpms_quat_to_cart3(qpms_quat_mult(q,
 191                                 qpms_quat_mult(vv, qc)));
 192 }
 193
 194 /// Versor quaternion from rotation vector (norm of the vector is the rotation angle).
 195 static inline qpms_quat_t qpms_quat_from_rotvector(cart3_t v) {
 196         return qpms_quat_exp(qpms_quat_rscale(0.5,
 197                                 qpms_quat_from_cart3(v)));
 198 }
 199
 200 /// Wigner D matrix element from a rotator quaternion for integer \a l.
 201 /**
 202  * The D matrix are calculated using formulae (3), (4), (6), (7) from
 203  * http://moble.github.io/spherical_functions/WignerDMatrices.html
 204  */
 205 complex double qpms_wignerD_elem(qpms_quat_t q, qpms_l_t l,
 206                qpms_m_t mp, qpms_m_t m);
 207
 208 /// A VSWF representation element of the O(3) group.
 209 /**
 210  * TODO more doc.
 211  */
 212 complex double qpms_vswf_irot_elem_from_irot3(
 213                 const qpms_irot3_t q, ///< The O(3) element in the quaternion representation.
 214                 qpms_l_t l, qpms_m_t mp, qpms_m_t m,
 215                 bool pseudo ///< Determines the sign of improper rotations. True for magnetic waves, false otherwise.
 216                 );
 217
 218
 219 static inline int qpms_irot3_checkdet(const qpms_irot3_t p) {
 220         if (p.det != 1 && p.det != -1) abort();
 221         return 0;
 222 }
 223
 224 /// Improper rotation multiplication.
 225 static inline qpms_irot3_t qpms_irot3_mult(const qpms_irot3_t p, const qpms_irot3_t q) {
 226 #ifndef NDEBUG
 227         qpms_irot3_checkdet(p);
 228         qpms_irot3_checkdet(q);
 229 #endif
 230         const qpms_irot3_t r = {qpms_quat_normalise(qpms_quat_mult(p.rot, q.rot)), p.det*q.det};
 231         return r;
 232 }
 233
 234 /// Improper rotation inverse operation.
 235 static inline qpms_irot3_t qpms_irot3_inv(qpms_irot3_t p) {
 236 #ifndef NDEBUG
 237         qpms_irot3_checkdet(p);
 238 #endif
 239         p.rot = qpms_quat_inv(p.rot);
 240         return p;
 241 }
 242
 243 /// Improper rotation power \f$ p^n \f$.
 244 static inline qpms_irot3_t qpms_irot3_pow(const qpms_irot3_t p, int n) {
 245 #ifndef NDEBUG
 246         qpms_irot3_checkdet(p);
 247 #endif
 248         const qpms_irot3_t r = {qpms_quat_normalise(qpms_quat_pow(p.rot, n)),
 249                 p.det == -1 ? min1pow(n) : 1};
 250         return r;
 251 }
 252
 253 /// Test approximate equality of irot3.
 254 static inline bool qpms_irot3_isclose(const qpms_irot3_t p, const qpms_irot3_t q, double atol) {
 255         return qpms_quat_isclose2(p.rot, q.rot, atol) && p.det == q.det;
 256 }
 257
 258 /// Apply an improper rotation onto a 3d cartesian vector.
 259 static inline cart3_t qpms_irot3_apply_cart3(const qpms_irot3_t p, const cart3_t v) {
 260 #ifndef NDEBUG
 261         qpms_irot3_checkdet(p);
 262 #endif
 263         return cart3_scale(p.det, qpms_quat_rot_cart3(p.rot, v));
 264 }
 265
 266 // Some basic transformations with irot3 type
 267 /// Identity
 268 static const qpms_irot3_t QPMS_IROT3_IDENTITY = {{1, 0}, 1};
 269 /// \f$ \pi \f$ rotation around x axis.
 270 static const qpms_irot3_t QPMS_IROT3_XROT_PI =  {{0, I}, 1};
 271 /// \f$ \pi \f$ rotation around y axis.
 272 static const qpms_irot3_t QPMS_IROT3_YROT_PI =  {{0, 1}, 1};
 273 /// \f$ \pi \f$ rotation around z axis.
 274 static const qpms_irot3_t QPMS_IROT3_ZROT_PI =  {{I, 0}, 1};
 275 /// Spatial inversion.
 276 static const qpms_irot3_t QPMS_IROT3_INVERSION = {{1, 0}, -1};
 277 /// yz-plane mirror symmetry
 278 static const qpms_irot3_t QPMS_IROT3_XFLIP = {{0, I}, -1};
 279 /// xz-plane mirror symmetry
 280 static const qpms_irot3_t QPMS_IROT3_YFLIP = {{0, 1}, -1};
 281 /// xy-plane mirror symmetry
 282 static const qpms_irot3_t QPMS_IROT3_ZFLIP = {{I, 0}, -1};
 283
 284 /// versor representing rotation around z-axis.
 285 static inline qpms_quat_t qpms_quat_zrot_angle(double angle) {
 286         qpms_quat_t q =  {cexp(I*(angle/2)), 0};
 287         return q;
 288 }
 289
 290 /// versor representing rotation \f$ C_N^k \f$, i.e. of angle \f$ 2\pi k / N\f$ around z axis.
 291 static inline qpms_quat_t qpms_quat_zrot_Nk(double N, double k) {
 292         return qpms_quat_zrot_angle(2 * M_PI * k / N);
 293 }
 294
 295 /// Rotation around z-axis.
 296 static inline qpms_irot3_t qpms_irot3_zrot_angle(double angle) {
 297         qpms_irot3_t q = {qpms_quat_zrot_angle(angle), 1};
 298         return q;
 299 }
 300
 301 /// Rotation \f$ C_N^k \f$, i.e. of angle \f$ 2\pi k / N\f$ around z axis.
 302 static inline qpms_irot3_t qpms_irot3_zrot_Nk(double N, double k) {
 303         return qpms_irot3_zrot_angle(2 * M_PI * k / N);
 304 }
 305
 306 #endif //QPMS_WIGNER_H