diff --git a/vector.hpp b/vector.hpp
index 03198605..8ea6d205 100644
--- a/vector.hpp
+++ b/vector.hpp
@@ -17,11 +17,13 @@
 #ifndef CRUFT_UTIL_VECTOR_HPP
 #define CRUFT_UTIL_VECTOR_HPP
 
-#include "./coord/fwd.hpp"
-#include "./coord.hpp"
+#include "coord/fwd.hpp"
+#include "coord/ops.hpp"
+#include "coord.hpp"
 
-#include "maths.hpp"
+#include "debug.hpp"
 #include "json/fwd.hpp"
+#include "maths.hpp"
 
 #include <cstddef>
 #include <cmath>
@@ -80,8 +82,6 @@ namespace util {
     //-------------------------------------------------------------------------
     // given a vector find two vectors which produce an orthonormal basis.
     //
-    // we use frisvad's method, avoids explicit normalisation. a good
-    // alternative is hughes-moeller, but the paper is hard to find.
     template <typename T>
     std::pair<
         util::vector<3,T>,
@@ -89,6 +89,10 @@ namespace util {
     >
     make_basis (const util::vector<3,T> n)
     {
+#if 0
+        // frisvad's method avoids explicit normalisation. a good alternative
+        // is hughes-moeller, but the paper is hard to find.
+        CHECK (is_normalised (n));
 
         // avoid a singularity
         if (n.z < -T(0.9999999)) {
@@ -101,10 +105,37 @@ namespace util {
         const T a = 1 / (1 + n.z);
         const T b = -n.x * n.y * a;
 
-        return {
-            { 1 - n.x * n.x * a, b, -n.x },
-            { b, 1 - n.y * n.y * a, -n.y }
-        };
+        const util::vector<3,T> v0 { 1 - n.x * n.x * a, b, -n.x };
+        const util::vector<3,T> v1 { b, 1 - n.y * n.y * a, -n.y };
+
+        CHECK (is_normalised (v0));
+        CHECK (is_normalised (v1));
+
+        return { v0, v1 };
+#else
+        // huges-moeller isn't as fast, but is more accurate
+        if(util::abs (n.x) > util::abs (n.z))
+        {
+            // Normalization factor for b2
+            auto const a = rsqrt (n.x * n.x + n.y * n.y);
+            util::vector<3,T> b1 { -n.y * a, n.x * a, 0 };
+
+             // Cross product using that b2 has a zero component
+            util::vector<3,T> b0 { b1.y * n.z, -b1.x * n.z, b1.x * n.y - b1.y * n.x };
+
+            return { b0, b1 };
+        }
+        else
+        {
+            // Normalization factor for b2
+            auto const a = rsqrt (n.y * n.y + n.z * n.z);
+            util::vector<3,T> b1 { 0.0f, -n.z * a, n.y * a };
+            // Cross product using that b2 has a zero component
+            util::vector<3,T> b0 { b1.y * n.z - b1.z * n.y, b1.z * n.x, -b1.y * n.x };
+
+            return { b0, b1 };
+        }
+#endif
     }