+/*
+ * Compute a rotation quaternion between two vectors
+ *
+ * cos(θ) + u * sin(θ)
+ *
+ * where θ is the angle between the two vectors and u
+ * is a unit vector axis of rotation
+ */
+
+static inline void ao_quaternion_vectors_to_rotation(struct ao_quaternion *r,
+ const struct ao_quaternion *a,
+ const struct ao_quaternion *b)
+{
+ /*
+ * The cross product will point orthogonally to the two
+ * vectors, forming our rotation axis. The length will be
+ * sin(θ), so these values are already multiplied by that.
+ */
+
+ float x = a->y * b->z - a->z * b->y;
+ float y = a->z * b->x - a->x * b->z;
+ float z = a->x * b->y - a->y * b->x;
+
+ float s_2 = x*x + y*y + z*z;
+ float s = sqrtf(s_2);
+
+ /* cos(θ) = a · b / (|a| |b|).
+ *
+ * a and b are both unit vectors, so the divisor is one
+ */
+ float c = a->x*b->x + a->y*b->y + a->z*b->z;
+
+ float c_half = sqrtf ((1 + c) / 2);
+ float s_half = sqrtf ((1 - c) / 2);
+
+ /*
+ * Divide out the sine factor from the
+ * cross product, then multiply in the
+ * half sine factor needed for the quaternion
+ */
+ float s_scale = s_half / s;
+
+ r->x = x * s_scale;
+ r->y = y * s_scale;
+ r->z = z * s_scale;
+
+ r->r = c_half;
+
+ ao_quaternion_normalize(r, r);
+}
+