diff --git a/qpms/wigner.c b/qpms/wigner.c
new file mode 100644
index 0000000..779e405
--- /dev/null
+++ b/qpms/wigner.c
@@ -0,0 +1,57 @@
+#include "wigner.h"
+#include "tiny_inlines.h"
+#include "kahansum.h"
+#define WIGNER_ATOL (1e-15)
+
+#define MIN(a,b) ((a)<(b)?(a):b)
+#define MAX(a,b) ((a)>(b)?(a):b)
+
+complex double qpms_wignerD_elem(const qpms_quat_t R,
+    const qpms_l_t l, const qpms_m_t mp, const qpms_m_t m) {
+  // TODO do some optimisations... The combinatoric coeffs could be precomputed.
+  if (abs(m) > l || abs(mp) > l)  abort();//return 0; // It's safer to crash. :D
+  const double mRa = cabs(R.a), mRb = cabs(R.b);
+  if (mRa < WIGNER_ATOL) {
+    if (m != -mp) return 0;
+    else return min1pow(l+m) * cpow(R.b, 2*m);
+  } else if (mRb < WIGNER_ATOL) {
+    if (m != mp) return 0;
+    else return cpow(R.a, 2*m);
+  } else if (mRa >= mRb) {
+    complex double prefac = exp(.5*(lgamma(l+m+1) + lgamma(l-m+1)
+        - lgamma(l+mp+1) - lgamma(l-mp+1))) 
+      * pow(mRa, 2*l-2*m)
+      * cpow(R.a, mp+m)
+      * cpow(R.b, -mp+m);
+    double sum, sumc; kahaninit( &sum, &sumc);
+    // FIXME this is probably quite stupid way to calculate the sum
+    for(int rho=MAX(0, -m+mp); rho <= MIN(l-m, l+mp); ++rho)
+      kahanadd(&sum, &sumc,
+          exp( 
+            lgamma(l+mp+1) - lgamma(rho+1) - lgamma(l+mp-rho+1)
+            + lgamma(l-mp+1) - lgamma(l-rho-m+1) - lgamma(-mp+rho+m+1)
+          )
+          * pow(-mRb*mRb/(mRa*mRa), rho)
+      );
+    return prefac * sum;
+  } else { // (mRa < mRb)
+    complex double prefac = min1pow(l-m) * exp(.5*(lgamma(l+m+1) + lgamma(l-m+1)
+        - lgamma(l+mp+1) - lgamma(l-mp+1))) 
+      * pow(mRb, 2*l-2*m)
+      * cpow(R.a, mp+m)
+      * cpow(R.b, -mp+m);
+    double sum, sumc; kahaninit( &sum, &sumc);
+    // FIXME this is probably quite stupid way to calculate the sum
+    for(int rho=MAX(0, -m-mp); rho <= MIN(l-mp, l-m); ++rho)
+      kahanadd(&sum, &sumc,
+          exp( 
+            lgamma(l+mp+1) - lgamma(l-m-rho+1) - lgamma(mp+m+rho+1)
+            + lgamma(l-mp+1) - lgamma(rho+1) - lgamma(l-mp-rho+1)
+          )
+          * pow(-mRa*mRa/(mRb*mRb), rho)
+      );
+    return prefac * sum;
+  }
+}
+
+
diff --git a/qpms/wigner.h b/qpms/wigner.h
new file mode 100644
index 0000000..efcf0d9
--- /dev/null
+++ b/qpms/wigner.h
@@ -0,0 +1,88 @@
+/*! \file wigner.h
+ * \brief Quaternions and Wigner matrices
+ */
+#ifndef QPMS_WIGNER_H
+#define QPMS_WIGNER_H
+
+#include "qpms_types.h"
+
+/// Quaternion type.
+/**
+ * Internaly represented as a pair of complex numbers,
+ * \f$ Q_a = Q_1 + iQ_z, Q_b = Q_y + i Q_x\f$.
+ */
+typedef struct qpms_quat_t {
+	complex double a, b;
+} qpms_quat_t;
+
+/// Quaternion type as four doubles.
+typedef struct qpms_quat4d_t {
+	double c1, ci, cj, ck;
+} qpms_quat4d_t;
+
+/// Conversion from the 4*double to the 2*complex quaternion.
+// 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;
+}
+
+/// Quaternion multiplication.
+/**
+ * \f[ (P Q)_a = P_a Q_a - \bar P_b Q_b, \f]
+ * \f[ (P Q)_b = P_b Q_a + \bar P_a Q_b. \f]
+ */
+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;
+}
+
+/// Quaternion addition.
+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;
+}
+
+/// 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};
+	return r;
+}
+
+// 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 conjugation.
+static inline qpms_quat_t qpms_quat_conj(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) {
+	return sqrt(creal(q.a * conj(q.a) + q.b * conj(q.b)));
+}
+
+/// 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);
+}
+
+/// Wigner D matrix element from a rotator quaternion for integer l.
+/**
+ * The D matrix are calculated using formulae (3), (4), (6), (7) from
+ * http://moble.github.io/spherical_functions/WignerDMatrices.html
+ */
+complex double qpms_wignerD_elem(qpms_quat_t q, qpms_l_t l,
+	       qpms_m_t	mp, qpms_m_t m);
+
+#endif //QPMS_WIGNER_H