+/******************************************************************************
+ *
+ * Calculate the dot product of two vectors u and v
+ * Dot product u.v = |u||v|cos(theta)
+ * Where theta = angle between u and v
+ */
+static inline double
+dot (const Vector3D u, const Vector3D v)
+{
+ return (u.x * v.x) + (u.y * v.y) + (u.z * v.z);
+}
+
+/******************************************************************************
+ *
+ * Calculate the cross product of two vectors.
+ * Gives a vector perpendicular to u and v with magnitude |u||v|sin(theta)
+ * Where theta = angle between u and v
+ */
+static inline Vector3D
+cross (const Vector3D u, const Vector3D v)
+{
+ Vector3D result;
+
+ result.x = (u.y * v.z - u.z * v.y);
+ result.y = (u.z * v.x - u.x * v.z);
+ result.z = (u.x * v.y - u.y * v.x);
+
+ return result;
+}
+
+/******************************************************************************
+ *
+ * Add vector v to vector u
+ */
+static inline void
+add (Vector3D *u, const Vector3D v)
+{
+ u->x = u->x + v.x;
+ u->y = u->y + v.y;
+ u->z = u->z + v.z;
+}
+
+/******************************************************************************
+ *
+ * Subtract vector v from vector u
+ */
+static inline Vector3D
+subtract (const Vector3D u, const Vector3D v)
+{
+ Vector3D result;
+
+ result.x = u.x - v.x;
+ result.y = u.y - v.y;
+ result.z = u.z - v.z;
+
+ return result;
+}
+
+/******************************************************************************
+ *
+ * multiply vector v by scalar s
+ */
+static inline Vector3D
+scale (const Vector3D v, const double s)
+{
+ Vector3D result;
+
+ result.x = v.x * s;
+ result.y = v.y * s;
+ result.z = v.z * s;
+ return result;
+}
+
+/******************************************************************************
+ *
+ * normalise vector v
+ */
+static inline Vector3D
+normalise (const Vector3D v)
+{
+ Vector3D result;
+ double magnitude;
+
+ magnitude = sqrt (dot(v, v));
+
+ if (magnitude > 1e-300)
+ {
+ result = scale (v, 1.0 / magnitude);
+ }
+ else
+ {
+ printf("zero\n");
+ result = zero_vector;
+ }
+ return result;
+}
+
+/******************************************************************************
+ *
+ * Calculate the transform matrix for the given quaternion
+ */
+static void
+quaternion_transform (Quaternion q, GLdouble * transform)
+{
+ GLdouble x, y, z, w;
+ x = q.x;
+ y = q.y;
+ z = q.z;
+ w = q.w;
+
+ transform[0] = (w * w) + (x * x) - (y * y) - (z * z);
+ transform[1] = (2.0 * x * y) + (2.0 * w * z);
+ transform[2] = (2.0 * x * z) - (2.0 * w * y);
+ transform[3] = 0.0;
+
+ transform[4] = (2.0 * x * y) - (2.0 * w * z);
+ transform[5] = (w * w) - (x * x) + (y * y) - (z * z);
+ transform[6] = (2.0 * y * z) + (2.0 * w * x);
+ transform[7] = 0.0;
+
+ transform[8] = (2.0 * x * z) + (2.0 * w * y);
+ transform[9] = (2.0 * y * z) - (2.0 * w * x);
+ transform[10] = (w * w) - (x * x) - (y * y) + (z * z);
+ transform[11] = 0.0;
+
+ transform[12] = 0.0;
+ transform[13] = 0.0;
+ transform[14] = 0.0;
+ transform[15] = (w * w) + (x * x) + (y * y) + (z * z);
+}
+
+/******************************************************************************
+ *
+ * Apply a matrix transform to the given vector
+ */
+static inline Vector3D
+vector_transform (Vector3D u, GLdouble * t)
+{
+ Vector3D result;
+
+ result.x = (u.x * t[0] + u.y * t[4] + u.z * t[8] + 1.0 * t[12]);
+ result.y = (u.x * t[1] + u.y * t[5] + u.z * t[9] + 1.0 * t[13]);
+ result.z = (u.x * t[2] + u.y * t[6] + u.z * t[10] + 1.0 * t[14]);
+
+ return result;
+}
+
+/******************************************************************************
+ *
+ * Return a node that is on an arc between node1 and node2, where distance
+ * is the proportion of the distance from node1 to the total arc.
+ */
+static Vector3D
+partial (Vector3D node1, Vector3D node2, double distance)
+{
+ Vector3D result;
+ Vector3D rotation_axis;
+ GLdouble transformation[16];
+ double angle;
+ Quaternion rotation;
+
+ rotation_axis = normalise (cross (node1, node2));
+ angle = acos (dot (node1, node2)) * distance;
+
+ rotation.x = rotation_axis.x * sin (angle / 2.0);
+ rotation.y = rotation_axis.y * sin (angle / 2.0);
+ rotation.z = rotation_axis.z * sin (angle / 2.0);
+ rotation.w = cos (angle / 2.0);
+
+ quaternion_transform (rotation, transformation);
+
+ result = vector_transform (node1, transformation);
+
+ return result;
+}
+