diff --git a/CMakeLists.txt b/CMakeLists.txt
index c561e91fa5..3d431d2dcc 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -361,6 +361,7 @@ set(BOUT_SOURCES
   ${CMAKE_CURRENT_BINARY_DIR}/include/bout/revision.hxx
   ${CMAKE_CURRENT_BINARY_DIR}/include/bout/version.hxx
   ./include/bout/tokamak_coordinates.hxx
+  ./src/mesh/tokamak_coordinates.cxx
 )
 
 
diff --git a/include/bout/tokamak_coordinates.hxx b/include/bout/tokamak_coordinates.hxx
index 7a9833c025..cb84549424 100644
--- a/include/bout/tokamak_coordinates.hxx
+++ b/include/bout/tokamak_coordinates.hxx
@@ -10,9 +10,11 @@
 #include "bout/utils.hxx"
 #include "bout/vector2d.hxx"
 
+
 class TokamakCoordinates {
 
 private:
+
   Mesh& mesh_m;
   Field2D Rxy_m;
   Field2D Bpxy_m;
@@ -22,73 +24,17 @@ private:
   FieldMetric ShearFactor_m;
   FieldMetric dx_m;
 
-  void normalise(BoutReal Lbar, BoutReal Bbar, BoutReal ShearFactor) {
+  void normalise(BoutReal Lbar, BoutReal Bbar, BoutReal ShearFactor);
 
-    Rxy_m /= Lbar;
-    Bpxy_m /= Bbar;
-    Btxy_m /= Bbar;
-    Bxy_m /= Bbar;
-    hthe_m /= Lbar;
-    ShearFactor_m *= Lbar * Lbar * Bbar * ShearFactor;
-    dx_m /= Lbar * Lbar * Bbar;
-  }
+  BoutReal get_sign_of_bp();
 
-public:
-  TokamakCoordinates(Mesh& mesh) : mesh_m(mesh) {
 
-    mesh.get(Rxy_m, "Rxy");
-    //    mesh->get(Zxy, "Zxy");
-    mesh.get(Bpxy_m, "Bpxy");
-    mesh.get(Btxy_m, "Btxy");
-    mesh.get(Bxy_m, "Bxy");
-    mesh.get(hthe_m, "hthe");
-    mesh.get(ShearFactor_m, "sinty");
-    mesh.get(dx_m, "dpsi");
-  }
+public:
 
-  BoutReal get_sign_of_bp() {
-    if (min(Bpxy_m, true) < 0.0) {
-      return -1.0;
-    }
-    return 1.0;
-  }
+  TokamakCoordinates(Mesh& mesh);
 
   Coordinates* make_coordinates(const bool noshear, BoutReal Lbar, BoutReal Bbar,
-                                        BoutReal ShearFactor = 1.0) {
-
-    normalise(Lbar, Bbar, ShearFactor);
-
-    const BoutReal sign_of_bp = get_sign_of_bp();
-
-    if (noshear) {
-      ShearFactor_m = 0.0;
-    }
-
-    auto* coord = mesh_m.getCoordinates();
-
-    const auto g11 = SQ(Rxy_m * Bpxy_m);
-    const auto g22 = 1.0 / SQ(hthe_m);
-    const auto g33 = SQ(ShearFactor_m) * g11 + SQ(Bxy_m) / g11;
-    const auto g12 = 0.0;
-    const auto g13 = -ShearFactor_m * g11;
-    const auto g23 = -sign_of_bp * Btxy_m / (hthe_m * Bpxy_m * Rxy_m);
-
-    const auto g_11 = 1.0 / g11 + SQ(ShearFactor_m * Rxy_m);
-    const auto g_22 = SQ(Bxy_m * hthe_m / Bpxy_m);
-    const auto g_33 = Rxy_m * Rxy_m;
-    const auto g_12 = sign_of_bp * Btxy_m * hthe_m * ShearFactor_m * Rxy_m / Bpxy_m;
-    const auto g_13 = ShearFactor_m * Rxy_m * Rxy_m;
-    const auto g_23 = sign_of_bp * Btxy_m * hthe_m * Rxy_m / Bpxy_m;
-
-    coord->setMetricTensor(ContravariantMetricTensor(g11, g22, g33, g12, g13, g23),
-                           CovariantMetricTensor(g_11, g_22, g_33, g_12, g_13, g_23));
-
-    coord->setJ(hthe_m / Bpxy_m);
-    coord->setBxy(Bxy_m);
-    coord->setDx(dx_m);
-
-    return coord;
-  }
+                                        BoutReal ShearFactor = 1.0);
 
   const Field2D& Rxy() const { return Rxy_m; }
   const Field2D& Bpxy() const { return Bpxy_m; }
diff --git a/src/mesh/tokamak_coordinates.cxx b/src/mesh/tokamak_coordinates.cxx
new file mode 100644
index 0000000000..6613d61cb1
--- /dev/null
+++ b/src/mesh/tokamak_coordinates.cxx
@@ -0,0 +1,69 @@
+#include <bout/tokamak_coordinates.hxx>
+
+
+void TokamakCoordinates::normalise(BoutReal Lbar, BoutReal Bbar, BoutReal ShearFactor) {
+
+    Rxy_m /= Lbar;
+    Bpxy_m /= Bbar;
+    Btxy_m /= Bbar;
+    Bxy_m /= Bbar;
+    hthe_m /= Lbar;
+    ShearFactor_m *= Lbar * Lbar * Bbar * ShearFactor;
+    dx_m /= Lbar * Lbar * Bbar;
+  }
+
+  TokamakCoordinates::TokamakCoordinates(Mesh& mesh) : mesh_m(mesh) {
+
+    mesh.get(Rxy_m, "Rxy");
+    //    mesh->get(Zxy, "Zxy");
+    mesh.get(Bpxy_m, "Bpxy");
+    mesh.get(Btxy_m, "Btxy");
+    mesh.get(Bxy_m, "Bxy");
+    mesh.get(hthe_m, "hthe");
+    mesh.get(ShearFactor_m, "sinty");
+    mesh.get(dx_m, "dpsi");
+  }
+
+  BoutReal TokamakCoordinates::get_sign_of_bp() {
+    if (min(Bpxy_m, true) < 0.0) {
+      return -1.0;
+    }
+    return 1.0;
+  }
+
+  Coordinates* TokamakCoordinates::make_coordinates(const bool noshear, BoutReal Lbar,
+                                                    BoutReal Bbar, BoutReal ShearFactor) {
+
+    normalise(Lbar, Bbar, ShearFactor);
+
+    const BoutReal sign_of_bp = get_sign_of_bp();
+
+    if (noshear) {
+      ShearFactor_m = 0.0;
+    }
+
+    auto* coord = mesh_m.getCoordinates();
+
+    const auto g11 = SQ(Rxy_m * Bpxy_m);
+    const auto g22 = 1.0 / SQ(hthe_m);
+    const auto g33 = SQ(ShearFactor_m) * g11 + SQ(Bxy_m) / g11;
+    const auto g12 = 0.0;
+    const auto g13 = -ShearFactor_m * g11;
+    const auto g23 = -sign_of_bp * Btxy_m / (hthe_m * Bpxy_m * Rxy_m);
+
+    const auto g_11 = 1.0 / g11 + SQ(ShearFactor_m * Rxy_m);
+    const auto g_22 = SQ(Bxy_m * hthe_m / Bpxy_m);
+    const auto g_33 = Rxy_m * Rxy_m;
+    const auto g_12 = sign_of_bp * Btxy_m * hthe_m * ShearFactor_m * Rxy_m / Bpxy_m;
+    const auto g_13 = ShearFactor_m * Rxy_m * Rxy_m;
+    const auto g_23 = sign_of_bp * Btxy_m * hthe_m * Rxy_m / Bpxy_m;
+
+    coord->setMetricTensor(ContravariantMetricTensor(g11, g22, g33, g12, g13, g23),
+                           CovariantMetricTensor(g_11, g_22, g_33, g_12, g_13, g_23));
+
+    coord->setJ(hthe_m / Bpxy_m);
+    coord->setBxy(Bxy_m);
+    coord->setDx(dx_m);
+
+    return coord;
+  }