diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt
index eff61951..e049ebe7 100644
--- a/examples/CMakeLists.txt
+++ b/examples/CMakeLists.txt
@@ -336,6 +336,9 @@ target_link_libraries(threewindows OpenGL::GL glfw Freetype::Freetype)
 add_executable(rod rod.cpp)
 target_link_libraries(rod OpenGL::GL glfw Freetype::Freetype)
 
+add_executable(rod_with_normals rod_with_normals.cpp)
+target_link_libraries(rod_with_normals OpenGL::GL glfw Freetype::Freetype)
+
 add_executable(rod_pan rod_pan.cpp)
 target_link_libraries(rod_pan OpenGL::GL glfw Freetype::Freetype)
 
@@ -367,6 +370,9 @@ if(NOT APPLE)
   add_executable(geodesic geodesic.cpp)
   target_link_libraries(geodesic OpenGL::GL glfw Freetype::Freetype)
 
+  add_executable(geodesic_with_normals geodesic_with_normals.cpp)
+  target_link_libraries(geodesic_with_normals OpenGL::GL glfw Freetype::Freetype)
+
   add_executable(geodesic_ce geodesic_ce.cpp)
   target_link_libraries(geodesic_ce OpenGL::GL glfw Freetype::Freetype)
 endif()
diff --git a/examples/geodesic_with_normals.cpp b/examples/geodesic_with_normals.cpp
new file mode 100644
index 00000000..3bca7161
--- /dev/null
+++ b/examples/geodesic_with_normals.cpp
@@ -0,0 +1,72 @@
+/*
+ * Visualize a sequence of icosahedral geodesics, showing their normals
+ */
+
+#include <iostream>
+#include <fstream>
+#include <cmath>
+#include <array>
+#include <stdexcept>
+#include <string>
+
+#include <sm/vec>
+
+#include <mplot/Visual.h>
+#include <mplot/ColourMap.h>
+#include <mplot/GeodesicVisual.h>
+#include <mplot/NormalsVisual.h>
+
+int main()
+{
+    int rtn = -1;
+
+    mplot::Visual v(1024, 768, "Geodesic Polyhedra with normals");
+    v.showCoordArrows (true);
+    // Set the Visual to rotate about the nearest VisualModel (Change at runtime with Ctrl-k)
+    v.rotateAboutNearest (true);
+    // In this example, use the 'rotate about a scene vertical axis' mode
+    v.rotateAboutVertical (true);
+
+    try {
+        sm::vec<float, 3> offset = {};
+        sm::vec<float, 3> step = { 2.2f };
+
+        mplot::ColourMap<float> cm (mplot::ColourMapType::Jet);
+        int imax = 4;
+        for (int i = 0; i < imax; ++i) {
+            auto cl = cm.convert (i / static_cast<float>(imax - 1));
+            auto gv1 = std::make_unique<mplot::GeodesicVisual<float>> (offset + step * i, 0.9f);
+            v.bindmodel (gv1);
+            gv1->iterations = i;
+            std::string lbl = std::string("iterations = ") + std::to_string(i);
+            gv1->addLabel (lbl, {0, -1, 0}, mplot::TextFeatures(0.06f));
+            gv1->cm.setType (mplot::ColourMapType::Jet);
+            gv1->colour_bb = cl;
+            gv1->finalize();
+
+            // re-colour after construction
+            auto gv1p = v.addVisualModel (gv1);
+            float imax_mult = 1.0f / static_cast<float>(imax);
+            // sequential colouring:
+            size_t sz1 = gv1p->data.size();
+            gv1p->data.linspace (0.0f, 1+i * imax_mult, sz1);
+            gv1p->reinitColours();
+
+            gv1p->vertex_postprocess(); // creates the triangles and normals required for NormalsVisual
+
+            // Create an associate normals model
+            auto nrm = std::make_unique<mplot::NormalsVisual<>> (gv1p);
+            v.bindmodel (nrm);
+            nrm->finalize();
+            v.addVisualModel (nrm);
+        }
+
+        v.keepOpen();
+
+    } catch (const std::exception& e) {
+        std::cerr << "Caught exception: " << e.what() << std::endl;
+        rtn = -1;
+    }
+
+    return rtn;
+}
diff --git a/examples/rod_with_normals.cpp b/examples/rod_with_normals.cpp
new file mode 100644
index 00000000..aff1d57b
--- /dev/null
+++ b/examples/rod_with_normals.cpp
@@ -0,0 +1,72 @@
+/*
+ * Visualize a Rod
+ */
+#include <iostream>
+#include <fstream>
+#include <cmath>
+#include <array>
+#include <stdexcept>
+#include <string>
+
+#include <sm/vec>
+
+#include <mplot/Visual.h>
+#include <mplot/ColourMap.h>
+#include <mplot/RodVisual.h>
+#include <mplot/NormalsVisual.h>
+
+int main()
+{
+    int rtn = -1;
+
+    mplot::Visual v(1024, 768, "Visualization");
+    v.zNear = 0.001;
+    //v.showCoordArrows (true);
+    //v.coordArrowsInScene (true);
+    // For a white background:
+    v.backgroundWhite();
+    // Switch on a mix of diffuse/ambient lighting
+    v.lightingEffects(true);
+
+    try {
+        constexpr sm::vec<float, 3> colour1 = { 1.0, 0.0, 0.0 };
+        sm::vec<float, 3> offset = { 0.0, 0.0, 0.0 };
+        sm::vec<float, 3> start = { 0, 0, 0 };
+        sm::vec<float, 3> end = { 0.25, 0, 0 };
+        auto rvm = std::make_unique<mplot::RodVisual<>> (offset, start, end, 0.1f, colour1, colour1);
+        v.bindmodel (rvm);
+        rvm->use_oriented_tube = false;
+        rvm->finalize();
+        auto rvmp1 = v.addVisualModel (rvm);
+        rvmp1->vertex_postprocess();
+
+        auto nrm = std::make_unique<mplot::NormalsVisual<>> (rvmp1);
+        v.bindmodel (nrm);
+        nrm->finalize();
+        v.addVisualModel (nrm);
+
+        constexpr sm::vec<float, 3> colour2 = { 0.0, 0.9, 0.4 };
+        sm::vec<float, 3> start2 = { -0.1, 0.2, 0.6 };
+        sm::vec<float, 3> end2 = { 0.2, 0.4, 0.6 };
+        // You can reuse the unique_ptr rvm, once you've transferred ownership with v.addVisualModel (rvm)
+        rvm = std::make_unique<mplot::RodVisual<>>(offset, start2, end2, 0.05f, colour2);
+        v.bindmodel (rvm);
+        rvm->finalize();
+        auto rvmp2 = v.addVisualModel (rvm);
+        rvmp2->vertex_postprocess();
+
+        // Create an associate normals model
+        nrm = std::make_unique<mplot::NormalsVisual<>> (rvmp2);
+        v.bindmodel (nrm);
+        nrm->finalize();
+        v.addVisualModel (nrm);
+
+        v.keepOpen();
+
+    } catch (const std::exception& e) {
+        std::cerr << "Caught exception: " << e.what() << std::endl;
+        rtn = -1;
+    }
+
+    return rtn;
+}
diff --git a/maths b/maths
index 548e7a03..de5a7954 160000
--- a/maths
+++ b/maths
@@ -1 +1 @@
-Subproject commit 548e7a03193e24f875756a7957e4efbb3beaef30
+Subproject commit de5a7954bbbb0cfcb885c7f82a16693ff44eb51c
diff --git a/mplot/CMakeLists.txt b/mplot/CMakeLists.txt
index 83286a28..c5388ab1 100644
--- a/mplot/CMakeLists.txt
+++ b/mplot/CMakeLists.txt
@@ -76,6 +76,7 @@ install(
   HSVWheelVisual.h
   IcosaVisual.h
   LengthscaleVisual.h
+  NormalsVisual.h
   PointRowsMeshVisual.h
   PointRowsVisual.h
   PolarVisual.h
diff --git a/mplot/NormalsVisual.h b/mplot/NormalsVisual.h
new file mode 100644
index 00000000..cdac03e2
--- /dev/null
+++ b/mplot/NormalsVisual.h
@@ -0,0 +1,90 @@
+#pragma once
+
+/*!
+ * \file Declares NormalsVisual to visualize the normals of another VisualModel
+ */
+
+#include <array>
+#include <sm/vec>
+#include <mplot/colour.h>
+#include <mplot/VisualModel.h>
+
+namespace mplot {
+
+    //! A class to visualize normals for another model
+    template <int glver = mplot::gl::version_4_1>
+    class NormalsVisual : public VisualModel<glver>
+    {
+    public:
+        NormalsVisual(mplot::VisualModel<glver>* _mymodel)
+        {
+            this->mymodel = _mymodel;
+            this->viewmatrix = _mymodel->getViewMatrix();
+        }
+
+        void initializeVertices()
+        {
+            if (mymodel == nullptr) {
+                std::cout << "NormalsVisual: I have no model; returning\n";
+                return;
+            }
+
+            const float cone_r = this->thickness * this->scale_factor * 2.0f;
+            const float tube_r = this->thickness * this->scale_factor;
+            // Copy data out of my model...
+            std::vector<float> mymodelPositions = mymodel->getVertexPositions();
+            std::vector<float> mymodelNormals = mymodel->getVertexNormals();
+            std::vector<float> mymodelColors = {};
+            if (!singlecolour) { mymodelColors = mymodel->getVertexColors(); }
+            // And interpret it
+            auto vp = reinterpret_cast<const std::vector<sm::vec<float, 3>>*>(&mymodelPositions);
+            auto vn = reinterpret_cast<const std::vector<sm::vec<float, 3>>*>(&mymodelNormals);
+            auto vc = reinterpret_cast<const std::vector<std::array<float, 3>>*>(&mymodelColors);
+
+            for (uint32_t ii = 0; ii < vn->size(); ++ii) {
+                // (*vp)[ii] is position, (*vn)[ii] is normal
+                std::array<float, 3> _clr = clr;
+                if (!singlecolour) { _clr = (*vc)[ii]; }
+                this->computeArrow ((*vp)[ii], ((*vp)[ii] + (*vn)[ii] * this->scale_factor),
+                                    _clr, tube_r, this->arrowhead_prop, cone_r, this->shapesides);
+            }
+
+            std::array<uint32_t, 3> ti = {};
+            sm::vec<float, 3> nv = {};
+            sm::vec<float, 3> nvc = {};
+            sm::vec<float, 3> nvd = {};
+            sm::vec<float, 3> pos = {};
+            // We also have vp1 (public) and triangles (also public)
+            for (auto t : mymodel->triangles) {
+                std::tie(ti, nv, nvc, nvd) = t;
+                // Plot tn at mean location of ti
+                pos = (mymodel->vp1[ti[0]] + mymodel->vp1[ti[1]] + mymodel->vp1[ti[2]]) / 3.0f;
+                // Mesh triangle normals
+                this->computeArrow (pos, (pos + nv * this->scale_factor),
+                                    clr, tube_r, this->arrowhead_prop, cone_r, this->shapesides);
+                // Computed triangle normals
+                this->computeArrow (pos, (pos + nvc * this->scale_factor),
+                                    clrnc, tube_r, this->arrowhead_prop, cone_r, this->shapesides);
+                this->computeArrow (pos, (pos + nvd * scale_factor),
+                                    clrnd, tube_r, this->arrowhead_prop, cone_r, this->shapesides);
+            }
+        };
+
+        // The model for which we will plot normal vectors
+        mplot::VisualModel<glver>* mymodel = nullptr;
+        // How many sides to each normal vector
+        int shapesides = 12;
+        // thickness for the normal vectors
+        float thickness = 0.025f;
+        // What proportion of the arrow length should the arrowhead length be?
+        float arrowhead_prop = 0.25f;
+        // How much to linearly scale the size of the vector
+        float scale_factor = 0.1f;
+        // Vector single colour
+        bool singlecolour = false;
+        std::array<float, 3> clr = mplot::colour::grey20;
+        std::array<float, 3> clrnc = mplot::colour::grey60; // computed norm
+        std::array<float, 3> clrnd = mplot::colour::grey90; // computed norm
+    };
+
+} // namespace mplot
diff --git a/mplot/RodVisual.h b/mplot/RodVisual.h
index 510ec9f8..8bf95b7d 100644
--- a/mplot/RodVisual.h
+++ b/mplot/RodVisual.h
@@ -66,10 +66,13 @@ namespace mplot {
             }
         }
 
-        void update (const sm::vec<float, 3>& s, const sm::vec<float, 3>& e)
+        template <std::size_t N = 3> requires (N == 3) || (N == 4)
+        void update (const sm::vec<float, N>& s, const sm::vec<float, N>& e)
         {
-            this->start_coord = s;
-            this->end_coord = e;
+            for (std::size_t i = 0; i < 3; ++i) {
+                this->start_coord[i] = s[i];
+                this->end_coord[i] = e[i];
+            }
             this->reinit();
         }
 
@@ -81,7 +84,8 @@ namespace mplot {
         float radius = 1.0f;
         //! If true, use face_uz and face_uy to draw the tube, else get a square-ended tube
         bool use_oriented_tube = true;
-        //! Face directions for the oriented tube
+        //! Face directions for the oriented tube's *end* cap. Choose carefully so that face_uy ^
+        //! face_uz gives the normal for the end cap.
         sm::vec<float, 3> face_uy = sm::vec<float, 3>::uy();
         sm::vec<float, 3> face_uz = sm::vec<float, 3>::uz();
 
diff --git a/mplot/VisualModelBase.h b/mplot/VisualModelBase.h
index 6d69191f..1287342f 100644
--- a/mplot/VisualModelBase.h
+++ b/mplot/VisualModelBase.h
@@ -227,8 +227,8 @@ namespace mplot {
         // The edges that make up the same triangles as are shown with this->indices, but in terms of vp1.
         // Each edge must be two indices in *ascending numerical order*
         std::set<std::array<uint32_t, 2>> edges;
-        // Triangles too. Might be more useful than edges.
-        sm::vvec<std::tuple<std::array<uint32_t, 3>, sm::vec<float, 3>>> triangles;
+        // Triangles too. Might be more useful than edges. Triangle given as indices into vp1
+        sm::vvec<std::tuple<std::array<uint32_t, 3>, sm::vec<float, 3>, sm::vec<float, 3>, sm::vec<float, 3>>> triangles;
 
         // Return index of vp1 that is closest to scene_coord. Can use vp1_to_indices to find the indices
         // into vertexPositions and vertexNormals that this index in the topographic mesh relates to.
@@ -259,16 +259,16 @@ namespace mplot {
             return trivert;
         }
 
-        sm::vvec<uint32_t> neighbours (const uint32_t idx) const
+        sm::vvec<uint32_t> neighbours (const uint32_t _idx) const
         {
             sm::vvec<uint32_t> rtn;
-            // Search edges to find those that include idx and then pack up the other ends in a return object
+            // Search edges to find those that include _idx and then pack up the other ends in a return object
             for (auto e : this->edges) {
                 // we have e[0] and e[1]
-                if (e[0] == idx) {
+                if (e[0] == _idx) {
                     // neighb is e[1]
                     rtn.push_back (e[1]);
-                } else if (e[1] == idx) {
+                } else if (e[1] == _idx) {
                     // neighb is e[0]
                     rtn.push_back (e[0]);
                 }
@@ -276,13 +276,13 @@ namespace mplot {
             return rtn;
         }
 
-        sm::vvec<std::array<uint32_t, 3>> neighbour_triangles (const uint32_t idx) const
+        sm::vvec<std::array<uint32_t, 3>> neighbour_triangles (const uint32_t _idx) const
         {
             sm::vvec<std::array<uint32_t, 3>> rtn;
             for (auto t: this->triangles) {
-                auto [ti, tn] = t;
-                // If it includes idx, add it to rtn
-                if (ti[0] == idx || ti[1] == idx || ti[2] == idx) {
+                auto [ti, tn, tnc, tnd] = t;
+                // If it includes _idx, add it to rtn
+                if (ti[0] == _idx || ti[1] == _idx || ti[2] == _idx) {
                     rtn.push_back (ti);
                 }
             }
@@ -310,8 +310,8 @@ namespace mplot {
         find_triangle_crossing (const sm::vec<float, 3>& coord, const sm::vec<float, 3>& vdir) const
         {
             for (auto tri : triangles) {
-                auto [ti, tn] = tri;
-                auto [isect, p] = sm::algo::ray_tri_intersection<float> (this->vp1[ti[0]], this->vp1[ti[1]], this->vp1[ti[2]], coord, vdir);
+                auto [ti, tn, tnc, tnd] = tri;
+                auto [isect, p] = sm::algo::ray_tri_intersection<float> (this->vp1[ti[0]], this->vp1[ti[1]], this->vp1[ti[2]], coord - (vdir / 2.0f), vdir);
                 if (isect) { return {p, ti, tn}; }
             }
 
@@ -335,7 +335,7 @@ namespace mplot {
             sm::vec<float> other_n = {fmax, fmax, fmax};
             sm::vec<float> my_n = {fmax, fmax, fmax}; // debug
             for (auto tri : triangles) {
-                auto [ti, tn] = tri;
+                auto [ti, tn, tnc, tnd] = tri;
                 if (ti == not_this) {
                     if constexpr (debug_normals) { my_n = tn; }
                     continue;
@@ -360,7 +360,9 @@ namespace mplot {
         std::tuple<sm::vec<float, 3>, std::array<uint32_t, 3>, sm::vec<float, 3>>
         find_triangle_crossing (const sm::vec<float, 3>& coord) const
         {
-            return this->find_triangle_crossing (coord, (this->bb.mid() - coord));
+            sm::vec<float, 3> vdir = this->bb.mid() - coord;
+            vdir.renormalize();
+            return this->find_triangle_crossing (coord, vdir);
         }
 
         /*!
@@ -374,7 +376,8 @@ namespace mplot {
          */
         void vertex_postprocess() // make_neighbour_mesh() ?
         {
-            constexpr bool debug = true;
+            constexpr bool debug = false;
+            constexpr bool debug_reorder = false;
 
             if constexpr (debug) { std::cout << __func__ << " called\n"; }
             // For each vertex, search for other vertices that have the same or almost the same location
@@ -457,23 +460,34 @@ namespace mplot {
                 // Direct population of triangles
                 std::array<uint32_t, 3> t = { equiv_top[indices[i]], equiv_top[indices[i+1]], equiv_top[indices[i+2]] };
 
-                // The normal vector for this triangle is easy to get, as we're dealing with indices already
+                // The normal vector for this triangle could be obtained from the mesh normals, but
+                // we can't trust them (though they're easy to get, as we're dealing with indices
+                // already). However, use this to ensure that our triangle indices order is in
+                // agreement with mesh normal as far as direction goes.
                 sm::vec<float> trinorm = this->get_normal (indices[i]) + this->get_normal (indices[i+1]) + this->get_normal (indices[i+2]) ;
                 trinorm.renormalize();
-                if constexpr (debug) { std::cout << "Triangle normal: " << trinorm << std::endl; }
 
-                // Check rotational sense of triangles
-                sm::vec<float, 3> n = (vp1[t[1]] - vp1[t[0]]).cross (vp1[t[2]] - vp1[t[0]]);
+                // Compute trinorm as well and compare with the one from the mesh - perhaps it's
+                // different? We really want the right normal.
+                const sm::vec<float>& tv0 = vp1[t[0]];
+                const sm::vec<float>& tv1 = vp1[t[1]];
+                const sm::vec<float>& tv2 = vp1[t[2]];
+                sm::vec<float>  nx = (tv1 - tv0);
+                sm::vec<float>  ny = (tv2 - tv0);
+                sm::vec<float, 3> n = nx.cross (ny);
+                n.renormalize();
+
+                // Check rotational sense of triangles?
                 if (n.dot (trinorm) < 0.0f) {
                     // need to swap order in t:
-                    if constexpr (debug) { std::cout << "Triangle reordered (corners 1 and 2 switched)\n"; }
+                    if constexpr (debug_reorder) { std::cout << "Triangle reordered (corners 1 and 2 switched)\n"; }
                     uint32_t ti = t[2];
                     t[2] = t[1];
                     t[1] = ti;
+                    n = -n; // Also reverse n
                 }
 
-                this->triangles.push_back ({t, trinorm});
-                // If required, can store trinorm into a container like triangle_normals with a push_back.
+                this->triangles.push_back ({t, n, nx, ny}); // n is computed normal
             }
 
             if constexpr (debug) {
@@ -481,12 +495,11 @@ namespace mplot {
                     std::cout << "Edge: " << e[0] << "," << e[1] << std::endl;
                 }
                 for (auto t : this->triangles) {
-                    auto [ti, tn] = t;
+                    auto [ti, tn, tnc, tnd] = t;
                     std::cout << "Tri: " << ti[0] << "," << ti[1] << "," << ti[2] << ", norm " << tn << std::endl;
                 }
+                std::cout << this->edges.size() << " edges and " << this->triangles.size() << " triangles in model '" << this->name << "'\n";
             }
-
-            std::cout << this->edges.size() << " edges and " << this->triangles.size() << " triangles in " << this->name << "\n";
         }
 
         /**
@@ -869,6 +882,12 @@ namespace mplot {
         void twodimensional (const bool val) { this->flags.set (vm_bools::twodimensional, val); }
         bool twodimensional() const { return this->flags.test (vm_bools::twodimensional); }
 
+        //! Getter for vertex positions (for mplot::NormalsVisual)
+        std::vector<float> getVertexPositions() { return this->vertexPositions; }
+        //! Getter for vertex normals (for mplot::NormalsVisual)
+        std::vector<float> getVertexNormals() { return this->vertexNormals; }
+        std::vector<float> getVertexColors() { return this->vertexColors; }
+
     protected:
 
         //! The model-specific view matrix. Used to transform the pose of the model in the scene.
@@ -996,7 +1015,7 @@ namespace mplot {
          * \param start The start of the tube
          * \param end The end of the tube
          * \param _ux a vector in the x axis direction for the end face
-         * \param _uy a vector in the y axis direction
+         * \param _uy a vector in the y axis direction. _ux ^ _uy gives the end cap normal
          * \param colStart The tube starting colour
          * \param colEnd The tube's ending colour
          * \param r Radius of the tube
@@ -1015,7 +1034,7 @@ namespace mplot {
             sm::vec<float> vend = end;
 
             // v is a face normal
-            sm::vec<float> v = _uy.cross(_ux);
+            sm::vec<float> v = _ux.cross(_uy);
             v.renormalize();
 
             // If bounding box, populate different buffers:
@@ -1755,7 +1774,7 @@ namespace mplot {
             for (int i = 0; i < n_faces; ++i) { // For each face in the geodesic...
                 sm::vec<F, 3> norm = { F{0}, F{0}, F{0} };
                 for (auto vtx : geo.poly.faces[i]) { // For each vertex in face...
-                    norm += vtx; // Add to the face norm
+                    norm += geo.poly.vertices[vtx]; // Add to the face norm
                     this->vertex_push (geo.poly.vertices[vtx].as_float() * r + so, this->vertexPositions);
                 }
                 sm::vec<float, 3> nf = (norm / F{3}).as_float();