diff --git a/examples/6field-simple/elm_6f.cxx b/examples/6field-simple/elm_6f.cxx
index 043aa80f44..6a51fc4221 100644
--- a/examples/6field-simple/elm_6f.cxx
+++ b/examples/6field-simple/elm_6f.cxx
@@ -366,7 +366,7 @@ class Elm_6f : public PhysicsModel {
       result = Grad_par(f, loc);
 
       if (nonlinear) {
-        result -= bracket(Psi, f, bm_mag) * tokamak_coordinates_factory.get_Bxy();
+        result -= bracket(Psi, f, bm_mag) * tokamak_coordinates_factory.Bxy();
       }
     }
 
@@ -678,7 +678,7 @@ class Elm_6f : public PhysicsModel {
 
     if (noshear) {
       if (include_curvature) {
-        b0xcv.z += tokamak_coordinates_factory.get_ShearFactor() * b0xcv.x;
+        b0xcv.z += tokamak_coordinates_factory.ShearFactor() * b0xcv.x;
       }
     }
 
@@ -688,7 +688,7 @@ class Elm_6f : public PhysicsModel {
     if (not mesh->IncIntShear) {
       // Dimits style, using local coordinate system
       if (include_curvature) {
-        b0xcv.z += tokamak_coordinates_factory.get_ShearFactor() * b0xcv.x;
+        b0xcv.z += tokamak_coordinates_factory.ShearFactor() * b0xcv.x;
       }
     }
 
@@ -832,7 +832,7 @@ class Elm_6f : public PhysicsModel {
       dump.add(sp_length, "sp_length", 1);
     }
 
-    J0 = SI::mu0 * Lbar * J0 / tokamak_coordinates_factory.get_Bxy();
+    J0 = SI::mu0 * Lbar * J0 / tokamak_coordinates_factory.Bxy();
     P0 = P0 / (SI::kb * (Tibar + Tebar) * eV_K / 2. * Nbar * density);
 
     b0xcv.x /= Bbar;
@@ -904,9 +904,9 @@ class Elm_6f : public PhysicsModel {
         q95 = q95_input; // use a constant for test
       } else {
         if (local_q) {
-          q95 = abs(tokamak_coordinates_factory.get_hthe()
-                    * tokamak_coordinates_factory.get_Btxy()
-                    / tokamak_coordinates_factory.get_Bpxy())
+          q95 = abs(tokamak_coordinates_factory.hthe()
+                    * tokamak_coordinates_factory.Btxy()
+                    / tokamak_coordinates_factory.Bpxy())
                 * q_alpha;
         } else {
           output.write("\tUsing q profile from grid.\n");
@@ -1037,29 +1037,29 @@ class Elm_6f : public PhysicsModel {
 
     if (mesh->IncIntShear) {
       // BOUT-06 style, using d/dx = d/dpsi + I * d/dz
-      coord->setIntShiftTorsion(tokamak_coordinates_factory.get_ShearFactor());
+      coord->setIntShiftTorsion(tokamak_coordinates_factory.ShearFactor());
     }
 
     // Set B field vector
 
     B0vec.covariant = false;
     B0vec.x = 0.;
-    B0vec.y = tokamak_coordinates_factory.get_Bpxy() / tokamak_coordinates_factory.get_hthe();
+    B0vec.y = tokamak_coordinates_factory.Bpxy() / tokamak_coordinates_factory.hthe();
     B0vec.z = 0.;
 
     // Set V0vec field vector
 
     V0vec.covariant = false;
     V0vec.x = 0.;
-    V0vec.y = Vp0 / tokamak_coordinates_factory.get_hthe();
-    V0vec.z = Vt0 / tokamak_coordinates_factory.get_Rxy();
+    V0vec.y = Vp0 / tokamak_coordinates_factory.hthe();
+    V0vec.z = Vt0 / tokamak_coordinates_factory.Rxy();
 
     // Set V0eff field vector
 
     V0eff.covariant = false;
     V0eff.x = 0.;
-    V0eff.y = -(tokamak_coordinates_factory.get_Btxy() / (tokamak_coordinates_factory.get_Bxy() * tokamak_coordinates_factory.get_Bxy())) * (Vp0 * tokamak_coordinates_factory.get_Btxy() - Vt0 * tokamak_coordinates_factory.get_Bpxy()) / tokamak_coordinates_factory.get_hthe();
-    V0eff.z = (tokamak_coordinates_factory.get_Bpxy() / (tokamak_coordinates_factory.get_Bxy() * tokamak_coordinates_factory.get_Bxy())) * (Vp0 * tokamak_coordinates_factory.get_Btxy() - Vt0 * tokamak_coordinates_factory.get_Bpxy()) / tokamak_coordinates_factory.get_Rxy();
+    V0eff.y = -(tokamak_coordinates_factory.Btxy() / (tokamak_coordinates_factory.Bxy() * tokamak_coordinates_factory.Bxy())) * (Vp0 * tokamak_coordinates_factory.Btxy() - Vt0 * tokamak_coordinates_factory.Bpxy()) / tokamak_coordinates_factory.hthe();
+    V0eff.z = (tokamak_coordinates_factory.Bpxy() / (tokamak_coordinates_factory.Bxy() * tokamak_coordinates_factory.Bxy())) * (Vp0 * tokamak_coordinates_factory.Btxy() - Vt0 * tokamak_coordinates_factory.Bpxy()) / tokamak_coordinates_factory.Rxy();
 
     Pe.setBoundary("P");
     Pi.setBoundary("P");
@@ -1101,10 +1101,10 @@ class Elm_6f : public PhysicsModel {
     if (diamag && diamag_phi0) {
       if (experiment_Er) { // get phi0 from grid file
         mesh->get(phi0, "Phi_0");
-        phi0 /= tokamak_coordinates_factory.get_Bxy() * Lbar * Va;
+        phi0 /= tokamak_coordinates_factory.Bxy() * Lbar * Va;
       } else {
         // Stationary equilibrium plasma. ExB velocity balances diamagnetic drift
-        phi0 = -Upara0 * Pi0 / tokamak_coordinates_factory.get_Bxy() / N0;
+        phi0 = -Upara0 * Pi0 / tokamak_coordinates_factory.Bxy() / N0;
       }
       SAVE_ONCE(phi0);
     }
@@ -1114,7 +1114,7 @@ class Elm_6f : public PhysicsModel {
     SAVE_ONCE(J0, P0);
     SAVE_ONCE(density, Lbar, Bbar, Tbar);
     SAVE_ONCE(Tibar, Tebar, Nbar);
-    Field2D tmp = tokamak_coordinates_factory.get_Bxy();
+    Field2D tmp = tokamak_coordinates_factory.Bxy();
     SAVE_ONCE(Va, tmp);
     SAVE_ONCE(Ti0, Te0, N0);
 
@@ -1138,13 +1138,13 @@ class Elm_6f : public PhysicsModel {
       Field2D logn0 = laplace_alpha * N0;
       Field3D Ntemp;
       Ntemp = N0;
-      ubyn = U * tokamak_coordinates_factory.get_Bxy() / Ntemp;
+      ubyn = U * tokamak_coordinates_factory.Bxy() / Ntemp;
       // Phi should be consistent with U
       if (laplace_alpha <= 0.0) {
-        phi = phiSolver->solve(ubyn) / tokamak_coordinates_factory.get_Bxy();
+        phi = phiSolver->solve(ubyn) / tokamak_coordinates_factory.Bxy();
       } else {
         phiSolver->setCoefC(logn0);
-        phi = phiSolver->solve(ubyn) / tokamak_coordinates_factory.get_Bxy();
+        phi = phiSolver->solve(ubyn) / tokamak_coordinates_factory.Bxy();
       }
     }
 
@@ -1219,9 +1219,9 @@ class Elm_6f : public PhysicsModel {
 
     //  Field2D lap_temp=0.0;
     Field2D logn0 = laplace_alpha * N0;
-    ubyn = U * tokamak_coordinates_factory.get_Bxy() / N0;
+    ubyn = U * tokamak_coordinates_factory.Bxy() / N0;
     if (diamag) {
-      ubyn -= Upara0 / N0 * Delp2(Pi) / tokamak_coordinates_factory.get_Bxy();
+      ubyn -= Upara0 / N0 * Delp2(Pi) / tokamak_coordinates_factory.Bxy();
       mesh->communicate(ubyn);
       ubyn.applyBoundary();
     }
@@ -1229,7 +1229,7 @@ class Elm_6f : public PhysicsModel {
     if (laplace_alpha > 0.0) {
       phiSolver->setCoefC(logn0);
     }
-    phi = phiSolver->solve(ubyn) / tokamak_coordinates_factory.get_Bxy();
+    phi = phiSolver->solve(ubyn) / tokamak_coordinates_factory.Bxy();
 
     mesh->communicate(phi);
 
@@ -1325,9 +1325,9 @@ class Elm_6f : public PhysicsModel {
 
     if (compress0) {
       if (nonlinear) {
-        Vepar = Vipar - tokamak_coordinates_factory.get_Bxy() * (Jpar) / N_tmp * Vepara;
+        Vepar = Vipar - tokamak_coordinates_factory.Bxy() * (Jpar) / N_tmp * Vepara;
       } else {
-        Vepar = Vipar - tokamak_coordinates_factory.get_Bxy() * (Jpar) / N0 * Vepara;
+        Vepar = Vipar - tokamak_coordinates_factory.Bxy() * (Jpar) / N0 * Vepara;
         Vepar.applyBoundary();
         mesh->communicate(Vepar);
       }
@@ -1365,13 +1365,13 @@ class Elm_6f : public PhysicsModel {
       ddt(Psi) = 0.0;
 
       if (spitzer_resist) {
-        ddt(Psi) = -Grad_parP(tokamak_coordinates_factory.get_Bxy() * phi) / tokamak_coordinates_factory.get_Bxy() - eta_spitzer * Jpar;
+        ddt(Psi) = -Grad_parP(tokamak_coordinates_factory.Bxy() * phi) / tokamak_coordinates_factory.Bxy() - eta_spitzer * Jpar;
       } else {
-        ddt(Psi) = -Grad_parP(tokamak_coordinates_factory.get_Bxy() * phi) / tokamak_coordinates_factory.get_Bxy() - eta * Jpar;
+        ddt(Psi) = -Grad_parP(tokamak_coordinates_factory.Bxy() * phi) / tokamak_coordinates_factory.Bxy() - eta * Jpar;
       }
 
       if (diamag) {
-        ddt(Psi) -= bracket(tokamak_coordinates_factory.get_Bxy() * phi0, Psi, bm_exb); // Equilibrium flow
+        ddt(Psi) -= bracket(tokamak_coordinates_factory.Bxy() * phi0, Psi, bm_exb); // Equilibrium flow
       }
 
       // Hyper-resistivity
@@ -1388,18 +1388,18 @@ class Elm_6f : public PhysicsModel {
 
       ddt(U) = 0.0;
 
-      ddt(U) = -SQ(tokamak_coordinates_factory.get_Bxy()) * bracket(Psi, J0, bm_mag) * tokamak_coordinates_factory.get_Bxy(); // Grad j term
+      ddt(U) = -SQ(tokamak_coordinates_factory.Bxy()) * bracket(Psi, J0, bm_mag) * tokamak_coordinates_factory.Bxy(); // Grad j term
 
       ddt(U) += 2.0 * Upara1 * b0xcv * Grad(P); // curvature term
 
-      ddt(U) += SQ(tokamak_coordinates_factory.get_Bxy()) * Grad_parP(Jpar); // b dot grad j
+      ddt(U) += SQ(tokamak_coordinates_factory.Bxy()) * Grad_parP(Jpar); // b dot grad j
 
       if (diamag) {
-        ddt(U) -= bracket(tokamak_coordinates_factory.get_Bxy() * phi0, U, bm_exb); // Equilibrium flow
+        ddt(U) -= bracket(tokamak_coordinates_factory.Bxy() * phi0, U, bm_exb); // Equilibrium flow
       }
 
       if (nonlinear) {
-        ddt(U) -= bracket(tokamak_coordinates_factory.get_Bxy() * phi, U, bm_exb); // Advection
+        ddt(U) -= bracket(tokamak_coordinates_factory.Bxy() * phi, U, bm_exb); // Advection
       }
 
       // parallel hyper-viscous diffusion for vector potential
@@ -1451,18 +1451,18 @@ class Elm_6f : public PhysicsModel {
 
       ddt(Ni) = 0.0;
 
-      ddt(Ni) -= bracket(tokamak_coordinates_factory.get_Bxy() * phi, N0, bm_exb);
+      ddt(Ni) -= bracket(tokamak_coordinates_factory.Bxy() * phi, N0, bm_exb);
 
       if (diamag) {
-        ddt(Ni) -= bracket(tokamak_coordinates_factory.get_Bxy() * phi0, Ni, bm_exb); // Equilibrium flow
+        ddt(Ni) -= bracket(tokamak_coordinates_factory.Bxy() * phi0, Ni, bm_exb); // Equilibrium flow
       }
 
       if (nonlinear) {
-        ddt(Ni) -= bracket(tokamak_coordinates_factory.get_Bxy() * phi, Ni, bm_exb); // Advection
+        ddt(Ni) -= bracket(tokamak_coordinates_factory.Bxy() * phi, Ni, bm_exb); // Advection
       }
 
       if (compress0) {
-        ddt(Ni) -= N0 * tokamak_coordinates_factory.get_Bxy() * Grad_parP(Vipar / tokamak_coordinates_factory.get_Bxy());
+        ddt(Ni) -= N0 * tokamak_coordinates_factory.Bxy() * Grad_parP(Vipar / tokamak_coordinates_factory.Bxy());
       }
 
       // 4th order Parallel diffusion terms
@@ -1481,18 +1481,18 @@ class Elm_6f : public PhysicsModel {
 
       ddt(Ti) = 0.0;
 
-      ddt(Ti) -= bracket(tokamak_coordinates_factory.get_Bxy() * phi, Ti0, bm_exb);
+      ddt(Ti) -= bracket(tokamak_coordinates_factory.Bxy() * phi, Ti0, bm_exb);
 
       if (diamag) {
-        ddt(Ti) -= bracket(phi0 * tokamak_coordinates_factory.get_Bxy(), Ti, bm_exb); // Equilibrium flow
+        ddt(Ti) -= bracket(phi0 * tokamak_coordinates_factory.Bxy(), Ti, bm_exb); // Equilibrium flow
       }
 
       if (nonlinear) {
-        ddt(Ti) -= bracket(phi * tokamak_coordinates_factory.get_Bxy(), Ti, bm_exb); // Advection
+        ddt(Ti) -= bracket(phi * tokamak_coordinates_factory.Bxy(), Ti, bm_exb); // Advection
       }
 
       if (compress0) {
-        ddt(Ti) -= 2.0 / 3.0 * Ti0 * tokamak_coordinates_factory.get_Bxy() * Grad_parP(Vipar / tokamak_coordinates_factory.get_Bxy());
+        ddt(Ti) -= 2.0 / 3.0 * Ti0 * tokamak_coordinates_factory.Bxy() * Grad_parP(Vipar / tokamak_coordinates_factory.Bxy());
       }
 
       if (diffusion_par > 0.0) {
@@ -1517,18 +1517,18 @@ class Elm_6f : public PhysicsModel {
 
       ddt(Te) = 0.0;
 
-      ddt(Te) -= bracket(tokamak_coordinates_factory.get_Bxy() * phi, Te0, bm_exb);
+      ddt(Te) -= bracket(tokamak_coordinates_factory.Bxy() * phi, Te0, bm_exb);
 
       if (diamag) {
-        ddt(Te) -= bracket(tokamak_coordinates_factory.get_Bxy() * phi0, Te, bm_exb); // Equilibrium flow
+        ddt(Te) -= bracket(tokamak_coordinates_factory.Bxy() * phi0, Te, bm_exb); // Equilibrium flow
       }
 
       if (nonlinear) {
-        ddt(Te) -= bracket(tokamak_coordinates_factory.get_Bxy() * phi, Te, bm_exb); // Advection
+        ddt(Te) -= bracket(tokamak_coordinates_factory.Bxy() * phi, Te, bm_exb); // Advection
       }
 
       if (compress0) {
-        ddt(Te) -= 2.0 / 3.0 * Te0 * tokamak_coordinates_factory.get_Bxy() * Grad_parP(Vepar / tokamak_coordinates_factory.get_Bxy());
+        ddt(Te) -= 2.0 / 3.0 * Te0 * tokamak_coordinates_factory.Bxy() * Grad_parP(Vepar / tokamak_coordinates_factory.Bxy());
       }
 
       if (diffusion_par > 0.0) {
@@ -1551,14 +1551,14 @@ class Elm_6f : public PhysicsModel {
       ddt(Vipar) = 0.0;
 
       ddt(Vipar) -= Vipara * Grad_parP(P) / N0;
-      ddt(Vipar) += Vipara * bracket(Psi, P0, bm_mag) * tokamak_coordinates_factory.get_Bxy() / N0;
+      ddt(Vipar) += Vipara * bracket(Psi, P0, bm_mag) * tokamak_coordinates_factory.Bxy() / N0;
 
       if (diamag) {
-        ddt(Vipar) -= bracket(tokamak_coordinates_factory.get_Bxy() * phi0, Vipar, bm_exb);
+        ddt(Vipar) -= bracket(tokamak_coordinates_factory.Bxy() * phi0, Vipar, bm_exb);
       }
 
       if (nonlinear) {
-        ddt(Vipar) -= bracket(tokamak_coordinates_factory.get_Bxy() * phi, Vipar, bm_exb);
+        ddt(Vipar) -= bracket(tokamak_coordinates_factory.Bxy() * phi, Vipar, bm_exb);
       }
 
       // parallel hyper-viscous diffusion for vector potential
diff --git a/examples/conducting-wall-mode/cwm.cxx b/examples/conducting-wall-mode/cwm.cxx
index 1fc55597fe..c869ea69ec 100644
--- a/examples/conducting-wall-mode/cwm.cxx
+++ b/examples/conducting-wall-mode/cwm.cxx
@@ -141,10 +141,10 @@ class CWM : public PhysicsModel {
     // add evolving variables to the communication object
     SOLVE_FOR(rho, te);
 
-    Field2D Rxy = tokamak_coordinates_factory.get_Rxy();
-    Field2D Bpxy = tokamak_coordinates_factory.get_Bpxy();
-    Field2D Btxy = tokamak_coordinates_factory.get_Btxy();
-    Field2D hthe = tokamak_coordinates_factory.get_hthe();
+    Field2D Rxy = tokamak_coordinates_factory.Rxy();
+    Field2D Bpxy = tokamak_coordinates_factory.Bpxy();
+    Field2D Btxy = tokamak_coordinates_factory.Btxy();
+    Field2D hthe = tokamak_coordinates_factory.hthe();
     SAVE_ONCE(Rxy, Bpxy, Btxy, Zxy, hthe);
     SAVE_ONCE(nu_hat, hthe0);
 
diff --git a/examples/dalf3/dalf3.cxx b/examples/dalf3/dalf3.cxx
index cb27d905b3..fdfc48db2a 100644
--- a/examples/dalf3/dalf3.cxx
+++ b/examples/dalf3/dalf3.cxx
@@ -162,7 +162,7 @@ class DALF3 : public PhysicsModel {
 
     if (lowercase(ptstr) == "shifted") {
       // Dimits style, using local coordinate system
-      b0xcv.z += tokamak_coordinates_factory.get_ShearFactor() * b0xcv.x;
+      b0xcv.z += tokamak_coordinates_factory.ShearFactor() * b0xcv.x;
       noshear = true;
     }
 
diff --git a/examples/elm-pb-outerloop/elm_pb_outerloop.cxx b/examples/elm-pb-outerloop/elm_pb_outerloop.cxx
index 5d6f1dca4f..4ac9774ac2 100644
--- a/examples/elm-pb-outerloop/elm_pb_outerloop.cxx
+++ b/examples/elm-pb-outerloop/elm_pb_outerloop.cxx
@@ -718,7 +718,7 @@ class ELMpb : public PhysicsModel {
     }
     const auto& metric = tokamak_coordinates_factory.make_tokamak_coordinates(noshear, Lbar, Bbar);
 
-    V0 = -tokamak_coordinates_factory.get_Rxy() * tokamak_coordinates_factory.get_Bpxy() * Dphi0 / tokamak_coordinates_factory.get_Bxy();
+    V0 = -tokamak_coordinates_factory.Rxy() * tokamak_coordinates_factory.Bpxy() * Dphi0 / tokamak_coordinates_factory.Bxy();
 
     if (simple_rmp) {
       include_rmp = true;
@@ -813,7 +813,7 @@ class ELMpb : public PhysicsModel {
 
     if (noshear) {
       if (include_curvature) {
-        b0xcv.z += tokamak_coordinates_factory.get_ShearFactor() * b0xcv.x;
+        b0xcv.z += tokamak_coordinates_factory.ShearFactor() * b0xcv.x;
       }
     }
 
@@ -823,7 +823,7 @@ class ELMpb : public PhysicsModel {
     if (not mesh->IncIntShear) {
       // Dimits style, using local coordinate system
       if (include_curvature) {
-        b0xcv.z += tokamak_coordinates_factory.get_ShearFactor() * b0xcv.x;
+        b0xcv.z += tokamak_coordinates_factory.ShearFactor() * b0xcv.x;
       }
     }
 
@@ -933,7 +933,7 @@ class ELMpb : public PhysicsModel {
     Field2D Te;
     Te = P0 / (2.0 * density * 1.602e-19); // Temperature in eV
 
-    J0 = -MU0 * Lbar * J0 / tokamak_coordinates_factory.get_Bxy();
+    J0 = -MU0 * Lbar * J0 / tokamak_coordinates_factory.Bxy();
     P0 = 2.0 * MU0 * P0 / (Bbar * Bbar);
     V0 = V0 / Va;
     Dphi0 *= Tbar;
@@ -1072,15 +1072,17 @@ class ELMpb : public PhysicsModel {
 
     B0vec.covariant = false;
     B0vec.x = 0.;
-    B0vec.y = tokamak_coordinates_factory.get_Bpxy() / tokamak_coordinates_factory.get_hthe();
+    B0vec.y = tokamak_coordinates_factory.Bpxy() / tokamak_coordinates_factory.hthe();
     B0vec.z = 0.;
 
     V0net.covariant = false; // presentation for net flow
     V0net.x = 0.;
-    V0net.y = tokamak_coordinates_factory.get_Rxy() * tokamak_coordinates_factory.get_Btxy() * tokamak_coordinates_factory.get_Bpxy() / (tokamak_coordinates_factory.get_hthe() * tokamak_coordinates_factory.get_Bxy() * tokamak_coordinates_factory.get_Bxy()) * Dphi0;
+    V0net.y =
+        tokamak_coordinates_factory.Rxy() * tokamak_coordinates_factory.Btxy() * tokamak_coordinates_factory.Bpxy() / (tokamak_coordinates_factory.hthe() * tokamak_coordinates_factory.Bxy() * tokamak_coordinates_factory.Bxy()) * Dphi0;
     V0net.z = -Dphi0;
 
-    U0 = B0vec * Curl(V0net) / tokamak_coordinates_factory.get_Bxy(); // get 0th vorticity for Kelvin-Holmholtz term
+    U0 = B0vec * Curl(V0net) / tokamak_coordinates_factory
+               .Bxy(); // get 0th vorticity for Kelvin-Holmholtz term
 
     /**************** SET EVOLVING VARIABLES *************/
 
@@ -1103,8 +1105,8 @@ class ELMpb : public PhysicsModel {
       SOLVE_FOR(Vpar);
       comms.add(Vpar);
 
-      beta = tokamak_coordinates_factory.get_Bxy() * tokamak_coordinates_factory.get_Bxy() / (0.5 + (tokamak_coordinates_factory.get_Bxy() * tokamak_coordinates_factory.get_Bxy() / (g * P0)));
-      gradparB = Grad_par(tokamak_coordinates_factory.get_Bxy()) / tokamak_coordinates_factory.get_Bxy();
+      beta = tokamak_coordinates_factory.Bxy() * tokamak_coordinates_factory.Bxy() / (0.5 + (tokamak_coordinates_factory.Bxy() * tokamak_coordinates_factory.Bxy() / (g * P0)));
+      gradparB = Grad_par(tokamak_coordinates_factory.Bxy()) / tokamak_coordinates_factory.Bxy();
 
       output.write("Beta in range {:e} -> {:e}\n", min(beta), max(beta));
     } else {
@@ -1137,10 +1139,10 @@ class ELMpb : public PhysicsModel {
     // Diamagnetic phi0
     if (diamag_phi0) {
       if (constn0) {
-        phi0 = -0.5 * dnorm * P0 / tokamak_coordinates_factory.get_Bxy();
+        phi0 = -0.5 * dnorm * P0 / tokamak_coordinates_factory.Bxy();
       } else {
         // Stationary equilibrium plasma. ExB velocity balances diamagnetic drift
-        phi0 = -0.5 * dnorm * P0 / tokamak_coordinates_factory.get_Bxy() / N0;
+        phi0 = -0.5 * dnorm * P0 / tokamak_coordinates_factory.Bxy() / N0;
       }
       SAVE_ONCE(phi0);
     } else {
@@ -1151,7 +1153,7 @@ class ELMpb : public PhysicsModel {
     // everything needed to recover physical units
     SAVE_ONCE(J0, P0);
     SAVE_ONCE(density, Lbar, Bbar, Tbar);
-    Field2D Bxy = tokamak_coordinates_factory.get_Bxy();
+    Field2D Bxy = tokamak_coordinates_factory.Bxy();
     SAVE_ONCE(Va, Bxy);
     SAVE_ONCE(Dphi0, U0);
     SAVE_ONCE(V0);
@@ -1211,7 +1213,7 @@ class ELMpb : public PhysicsModel {
 
     if (mesh->IncIntShear) {
       // BOUT-06 style, using d/dx = d/dpsi + I * d/dz
-      metric->setIntShiftTorsion(tokamak_coordinates_factory.get_ShearFactor());
+      metric->setIntShiftTorsion(tokamak_coordinates_factory.ShearFactor());
 
     }
 
@@ -1525,7 +1527,7 @@ class ELMpb : public PhysicsModel {
     auto P0_acc = Field2DAccessor<>(P0);
     auto J0_acc = Field2DAccessor<>(J0);
     auto phi0_acc = Field2DAccessor<>(phi0);
-    Field2D Bxy = tokamak_coordinates_factory.get_Bxy();
+    Field2D Bxy = tokamak_coordinates_factory.Bxy();
     auto B0_acc = Field2DAccessor<>(Bxy);
 
     // Evolving fields
@@ -1713,16 +1715,16 @@ class ELMpb : public PhysicsModel {
     // Vorticity equation
 
     // Grad j term
-    // ddt(U) = SQ(tokamak_coordinates_factory.get_Bxy()) * b0xGrad_dot_Grad(Psi, J0, CELL_CENTRE);
+    // ddt(U) = SQ(tokamak_coordinates_factory.Bxy()) * b0xGrad_dot_Grad(Psi, J0, CELL_CENTRE);
 
     // if (include_rmp) {
-    //   ddt(U) += SQ(tokamak_coordinates_factory.get_Bxy()) * b0xGrad_dot_Grad(rmp_Psi, J0, CELL_CENTRE);
+    //   ddt(U) += SQ(tokamak_coordinates_factory.Bxy()) * b0xGrad_dot_Grad(rmp_Psi, J0, CELL_CENTRE);
     // }
 
     ddt(U) += b0xcv * Grad(P); // curvature term
 
     // if (!nogradparj) { // Parallel current term
-    //   ddt(U) -= SQ(tokamak_coordinates_factory.get_Bxy()) * Grad_parP(Jpar, CELL_CENTRE); // b dot grad j
+    //   ddt(U) -= SQ(tokamak_coordinates_factory.Bxy()) * Grad_parP(Jpar, CELL_CENTRE); // b dot grad j
     // }
 
     if (withflow && K_H_term) { // K_H_term
diff --git a/examples/elm-pb/elm_pb.cxx b/examples/elm-pb/elm_pb.cxx
index 03bac5c569..3a3ac0fe9c 100644
--- a/examples/elm-pb/elm_pb.cxx
+++ b/examples/elm-pb/elm_pb.cxx
@@ -645,7 +645,7 @@ class ELMpb : public PhysicsModel {
       Dphi0 *= -1;
     }
 
-    V0 = -tokamak_coordinates_factory.get_Rxy() * tokamak_coordinates_factory.get_Bpxy() * Dphi0 / tokamak_coordinates_factory.get_Bxy();
+    V0 = -tokamak_coordinates_factory.Rxy() * tokamak_coordinates_factory.Bpxy() * Dphi0 / tokamak_coordinates_factory.Bxy();
 
     if (simple_rmp) {
       include_rmp = true;
@@ -740,7 +740,7 @@ class ELMpb : public PhysicsModel {
 
     if (noshear) {
       if (include_curvature) {
-        b0xcv.z += tokamak_coordinates_factory.get_ShearFactor() * b0xcv.x;
+        b0xcv.z += tokamak_coordinates_factory.ShearFactor() * b0xcv.x;
       }
     }
 
@@ -750,7 +750,7 @@ class ELMpb : public PhysicsModel {
     if (not mesh->IncIntShear) {
       // Dimits style, using local coordinate system
       if (include_curvature) {
-        b0xcv.z += tokamak_coordinates_factory.get_ShearFactor() * b0xcv.x;
+        b0xcv.z += tokamak_coordinates_factory.ShearFactor() * b0xcv.x;
       }
     }
 
@@ -867,7 +867,7 @@ class ELMpb : public PhysicsModel {
     Field2D Te;
     Te = P0 / (2.0 * density * 1.602e-19); // Temperature in eV
 
-    J0 = -MU0 * Lbar * J0 / tokamak_coordinates_factory.get_Bxy();
+    J0 = -MU0 * Lbar * J0 / tokamak_coordinates_factory.Bxy();
     P0 = 2.0 * MU0 * P0 / (Bbar * Bbar);
     V0 = V0 / Va;
     Dphi0 *= Tbar;
@@ -884,7 +884,7 @@ class ELMpb : public PhysicsModel {
 
     if (mesh->IncIntShear) {
       // BOUT-06 style, using d/dx = d/dpsi + I * d/dz
-      metric->setIntShiftTorsion(tokamak_coordinates_factory.get_ShearFactor());
+      metric->setIntShiftTorsion(tokamak_coordinates_factory.ShearFactor());
     }
 
     if (constn0) {
@@ -1016,15 +1016,17 @@ class ELMpb : public PhysicsModel {
 
     B0vec.covariant = false;
     B0vec.x = 0.;
-    B0vec.y = tokamak_coordinates_factory.get_Bpxy() / tokamak_coordinates_factory.get_hthe();
+    B0vec.y = tokamak_coordinates_factory.Bpxy() / tokamak_coordinates_factory.hthe();
     B0vec.z = 0.;
 
     V0net.covariant = false; // presentation for net flow
     V0net.x = 0.;
-    V0net.y = tokamak_coordinates_factory.get_Rxy() * tokamak_coordinates_factory.get_Btxy() * tokamak_coordinates_factory.get_Bpxy() / (tokamak_coordinates_factory.get_hthe() * tokamak_coordinates_factory.get_Bxy() * tokamak_coordinates_factory.get_Bxy()) * Dphi0;
+    V0net.y =
+        tokamak_coordinates_factory.Rxy() * tokamak_coordinates_factory.Btxy() * tokamak_coordinates_factory.Bpxy() / (tokamak_coordinates_factory.hthe() * tokamak_coordinates_factory.Bxy() * tokamak_coordinates_factory.Bxy()) * Dphi0;
     V0net.z = -Dphi0;
 
-    U0 = B0vec * Curl(V0net) / tokamak_coordinates_factory.get_Bxy(); // get 0th vorticity for Kelvin-Holmholtz term
+    U0 = B0vec * Curl(V0net) / tokamak_coordinates_factory
+               .Bxy(); // get 0th vorticity for Kelvin-Holmholtz term
 
     /**************** SET EVOLVING VARIABLES *************/
 
@@ -1047,8 +1049,8 @@ class ELMpb : public PhysicsModel {
       SOLVE_FOR(Vpar);
       comms.add(Vpar);
 
-      beta = tokamak_coordinates_factory.get_Bxy() * tokamak_coordinates_factory.get_Bxy() / (0.5 + (tokamak_coordinates_factory.get_Bxy() * tokamak_coordinates_factory.get_Bxy() / (g * P0)));
-      gradparB = Grad_par(tokamak_coordinates_factory.get_Bxy()) / tokamak_coordinates_factory.get_Bxy();
+      beta = tokamak_coordinates_factory.Bxy() * tokamak_coordinates_factory.Bxy() / (0.5 + (tokamak_coordinates_factory.Bxy() * tokamak_coordinates_factory.Bxy() / (g * P0)));
+      gradparB = Grad_par(tokamak_coordinates_factory.Bxy()) / tokamak_coordinates_factory.Bxy();
 
       output.write("Beta in range {:e} -> {:e}\n", min(beta), max(beta));
     } else {
@@ -1081,10 +1083,10 @@ class ELMpb : public PhysicsModel {
     // Diamagnetic phi0
     if (diamag_phi0) {
       if (constn0) {
-        phi0 = -0.5 * dnorm * P0 / tokamak_coordinates_factory.get_Bxy();
+        phi0 = -0.5 * dnorm * P0 / tokamak_coordinates_factory.Bxy();
       } else {
         // Stationary equilibrium plasma. ExB velocity balances diamagnetic drift
-        phi0 = -0.5 * dnorm * P0 / tokamak_coordinates_factory.get_Bxy() / N0;
+        phi0 = -0.5 * dnorm * P0 / tokamak_coordinates_factory.Bxy() / N0;
       }
       SAVE_ONCE(phi0);
     }
@@ -1093,7 +1095,7 @@ class ELMpb : public PhysicsModel {
     // everything needed to recover physical units
     SAVE_ONCE(J0, P0);
     SAVE_ONCE(density, Lbar, Bbar, Tbar);
-    Field2D tmp = tokamak_coordinates_factory.get_Bxy();
+    Field2D tmp = tokamak_coordinates_factory.Bxy();
     SAVE_ONCE(Va, tmp);
     SAVE_ONCE(Dphi0, U0);
     SAVE_ONCE(V0);
@@ -1161,7 +1163,7 @@ class ELMpb : public PhysicsModel {
     Field3D result = Grad_par(f, loc);
 
     if (nonlinear) {
-      const auto Bxy = tokamak_coordinates_factory.get_Bxy();
+      const auto Bxy = tokamak_coordinates_factory.Bxy();
       result -= bracket(interp_to(Psi, loc), f, bm_mag) * Bxy;
 
       if (include_rmp) {
@@ -1301,12 +1303,12 @@ class ELMpb : public PhysicsModel {
         }
 
         if (diamag) {
-          phi -= 0.5 * dnorm * P / tokamak_coordinates_factory.get_Bxy();
+          phi -= 0.5 * dnorm * P / tokamak_coordinates_factory.Bxy();
         }
       } else {
         ubyn = U / N0;
         if (diamag) {
-          ubyn -= 0.5 * dnorm / (N0 * tokamak_coordinates_factory.get_Bxy()) * Delp2(P);
+          ubyn -= 0.5 * dnorm / (N0 * tokamak_coordinates_factory.Bxy()) * Delp2(P);
           mesh->communicate(ubyn);
         }
         // Invert laplacian for phi
@@ -1451,9 +1453,9 @@ class ELMpb : public PhysicsModel {
 
     if (evolve_jpar) {
       // Jpar
-      Field3D B0U = tokamak_coordinates_factory.get_Bxy() * U;
+      Field3D B0U = tokamak_coordinates_factory.Bxy() * U;
       mesh->communicate(B0U);
-      ddt(Jpar) = -Grad_parP(B0U, loc) / tokamak_coordinates_factory.get_Bxy() + eta * Delp2(Jpar);
+      ddt(Jpar) = -Grad_parP(B0U, loc) / tokamak_coordinates_factory.Bxy() + eta * Delp2(Jpar);
 
       if (relax_j_vac) {
         // Make ddt(Jpar) relax to zero.
@@ -1462,16 +1464,16 @@ class ELMpb : public PhysicsModel {
       }
     } else {
       // Vector potential
-      ddt(Psi) = -Grad_parP(phi * tokamak_coordinates_factory.get_Bxy(), loc) / tokamak_coordinates_factory.get_Bxy() + eta * Jpar;
+      ddt(Psi) = -Grad_parP(phi * tokamak_coordinates_factory.Bxy(), loc) / tokamak_coordinates_factory.Bxy() + eta * Jpar;
 
       if (eHall) { // electron parallel pressure
         ddt(Psi) += 0.25 * delta_i
-                    * (Grad_parP(tokamak_coordinates_factory.get_Bxy() * P, loc) / tokamak_coordinates_factory.get_Bxy()
-                       + bracket(interp_to(P0, loc), Psi, bm_mag) * tokamak_coordinates_factory.get_Bxy());
+                    * (Grad_parP(tokamak_coordinates_factory.Bxy() * P, loc) / tokamak_coordinates_factory.Bxy()
+                       + bracket(interp_to(P0, loc), Psi, bm_mag) * tokamak_coordinates_factory.Bxy());
       }
 
       if (diamag_phi0) { // Equilibrium flow
-        ddt(Psi) -= bracket(interp_to(phi0, loc), Psi, bm_exb) * tokamak_coordinates_factory.get_Bxy();
+        ddt(Psi) -= bracket(interp_to(phi0, loc), Psi, bm_exb) * tokamak_coordinates_factory.Bxy();
       }
 
       if (withflow) { // net flow
@@ -1479,7 +1481,7 @@ class ELMpb : public PhysicsModel {
       }
 
       if (diamag_grad_t) { // grad_par(T_e) correction
-        ddt(Psi) += 1.71 * dnorm * 0.5 * Grad_parP(P, loc) / tokamak_coordinates_factory.get_Bxy();
+        ddt(Psi) += 1.71 * dnorm * 0.5 * Grad_parP(P, loc) / tokamak_coordinates_factory.Bxy();
       }
 
       if (hyperresist > 0.0) { // Hyper-resistivity
@@ -1517,15 +1519,15 @@ class ELMpb : public PhysicsModel {
     Psi_loc = interp_to(Psi, CELL_CENTRE, "RGN_ALL");
     Psi_loc.applyBoundary();
     // Grad j term
-    ddt(U) = SQ(tokamak_coordinates_factory.get_Bxy()) * b0xGrad_dot_Grad(Psi_loc, J0, CELL_CENTRE);
+    ddt(U) = SQ(tokamak_coordinates_factory.Bxy()) * b0xGrad_dot_Grad(Psi_loc, J0, CELL_CENTRE);
     if (include_rmp) {
-      ddt(U) += SQ(tokamak_coordinates_factory.get_Bxy()) * b0xGrad_dot_Grad(rmp_Psi, J0, CELL_CENTRE);
+      ddt(U) += SQ(tokamak_coordinates_factory.Bxy()) * b0xGrad_dot_Grad(rmp_Psi, J0, CELL_CENTRE);
     }
 
     ddt(U) += b0xcv * Grad(P); // curvature term
 
     if (!nogradparj) {                                 // Parallel current term
-      ddt(U) -= SQ(tokamak_coordinates_factory.get_Bxy()) * Grad_parP(Jpar, CELL_CENTRE); // b dot grad j
+      ddt(U) -= SQ(tokamak_coordinates_factory.Bxy()) * Grad_parP(Jpar, CELL_CENTRE); // b dot grad j
     }
 
     if (withflow && K_H_term) { // K_H_term
@@ -1541,7 +1543,7 @@ class ELMpb : public PhysicsModel {
     }
 
     if (nonlinear) { // Advection
-      ddt(U) -= bracket(phi, U, bm_exb) * tokamak_coordinates_factory.get_Bxy();
+      ddt(U) -= bracket(phi, U, bm_exb) * tokamak_coordinates_factory.Bxy();
     }
 
     // Viscosity terms
@@ -1585,11 +1587,11 @@ class ELMpb : public PhysicsModel {
       Pi = 0.5 * P;
       Pi0 = 0.5 * P0;
 
-      Dperp2Phi0 = Field3D(Delp2(tokamak_coordinates_factory.get_Bxy() * phi0));
+      Dperp2Phi0 = Field3D(Delp2(tokamak_coordinates_factory.Bxy() * phi0));
       Dperp2Phi0.applyBoundary();
       mesh->communicate(Dperp2Phi0);
 
-      Dperp2Phi = Delp2(tokamak_coordinates_factory.get_Bxy() * phi);
+      Dperp2Phi = Delp2(tokamak_coordinates_factory.Bxy() * phi);
       Dperp2Phi.applyBoundary();
       mesh->communicate(Dperp2Phi);
 
@@ -1601,35 +1603,35 @@ class ELMpb : public PhysicsModel {
       Dperp2Pi.applyBoundary();
       mesh->communicate(Dperp2Pi);
 
-      bracketPhi0P = bracket(tokamak_coordinates_factory.get_Bxy() * phi0, Pi, bm_exb);
+      bracketPhi0P = bracket(tokamak_coordinates_factory.Bxy() * phi0, Pi, bm_exb);
       bracketPhi0P.applyBoundary();
       mesh->communicate(bracketPhi0P);
 
-      bracketPhiP0 = bracket(tokamak_coordinates_factory.get_Bxy() * phi, Pi0, bm_exb);
+      bracketPhiP0 = bracket(tokamak_coordinates_factory.Bxy() * phi, Pi0, bm_exb);
       bracketPhiP0.applyBoundary();
       mesh->communicate(bracketPhiP0);
 
-      ddt(U) -= 0.5 * Upara2 * bracket(Pi, Dperp2Phi0, bm_exb) / tokamak_coordinates_factory.get_Bxy();
-      ddt(U) -= 0.5 * Upara2 * bracket(Pi0, Dperp2Phi, bm_exb) / tokamak_coordinates_factory.get_Bxy();
-      Field3D B0phi = tokamak_coordinates_factory.get_Bxy() * phi;
+      ddt(U) -= 0.5 * Upara2 * bracket(Pi, Dperp2Phi0, bm_exb) / tokamak_coordinates_factory.Bxy();
+      ddt(U) -= 0.5 * Upara2 * bracket(Pi0, Dperp2Phi, bm_exb) / tokamak_coordinates_factory.Bxy();
+      Field3D B0phi = tokamak_coordinates_factory.Bxy() * phi;
       mesh->communicate(B0phi);
-      Field3D B0phi0 = tokamak_coordinates_factory.get_Bxy() * phi0;
+      Field3D B0phi0 = tokamak_coordinates_factory.Bxy() * phi0;
       mesh->communicate(B0phi0);
-      ddt(U) += 0.5 * Upara2 * bracket(B0phi, Dperp2Pi0, bm_exb) / tokamak_coordinates_factory.get_Bxy();
-      ddt(U) += 0.5 * Upara2 * bracket(B0phi0, Dperp2Pi, bm_exb) / tokamak_coordinates_factory.get_Bxy();
-      ddt(U) -= 0.5 * Upara2 * Delp2(bracketPhi0P) / tokamak_coordinates_factory.get_Bxy();
-      ddt(U) -= 0.5 * Upara2 * Delp2(bracketPhiP0) / tokamak_coordinates_factory.get_Bxy();
+      ddt(U) += 0.5 * Upara2 * bracket(B0phi, Dperp2Pi0, bm_exb) / tokamak_coordinates_factory.Bxy();
+      ddt(U) += 0.5 * Upara2 * bracket(B0phi0, Dperp2Pi, bm_exb) / tokamak_coordinates_factory.Bxy();
+      ddt(U) -= 0.5 * Upara2 * Delp2(bracketPhi0P) / tokamak_coordinates_factory.Bxy();
+      ddt(U) -= 0.5 * Upara2 * Delp2(bracketPhiP0) / tokamak_coordinates_factory.Bxy();
 
       if (nonlinear) {
-        Field3D B0phi = tokamak_coordinates_factory.get_Bxy() * phi;
+        Field3D B0phi = tokamak_coordinates_factory.Bxy() * phi;
         mesh->communicate(B0phi);
         bracketPhiP = bracket(B0phi, Pi, bm_exb);
         bracketPhiP.applyBoundary();
         mesh->communicate(bracketPhiP);
 
-        ddt(U) -= 0.5 * Upara2 * bracket(Pi, Dperp2Phi, bm_exb) / tokamak_coordinates_factory.get_Bxy();
-        ddt(U) += 0.5 * Upara2 * bracket(B0phi, Dperp2Pi, bm_exb) / tokamak_coordinates_factory.get_Bxy();
-        ddt(U) -= 0.5 * Upara2 * Delp2(bracketPhiP) / tokamak_coordinates_factory.get_Bxy();
+        ddt(U) -= 0.5 * Upara2 * bracket(Pi, Dperp2Phi, bm_exb) / tokamak_coordinates_factory.Bxy();
+        ddt(U) += 0.5 * Upara2 * bracket(B0phi, Dperp2Pi, bm_exb) / tokamak_coordinates_factory.Bxy();
+        ddt(U) -= 0.5 * Upara2 * Delp2(bracketPhiP) / tokamak_coordinates_factory.Bxy();
       }
     }
 
@@ -1659,7 +1661,7 @@ class ELMpb : public PhysicsModel {
       }
 
       if (nonlinear) { // Advection
-        ddt(P) -= bracket(phi, P, bm_exb) * tokamak_coordinates_factory.get_Bxy();
+        ddt(P) -= bracket(phi, P, bm_exb) * tokamak_coordinates_factory.Bxy();
       }
     }
 
@@ -1704,7 +1706,7 @@ class ELMpb : public PhysicsModel {
       ddt(Vpar) = -0.5 * (Grad_par(P, loc) + Grad_par(P0, loc));
 
       if (nonlinear) {
-        ddt(Vpar) -= bracket(interp_to(phi, loc), Vpar, bm_exb) * tokamak_coordinates_factory.get_Bxy(); // Advection
+        ddt(Vpar) -= bracket(interp_to(phi, loc), Vpar, bm_exb) * tokamak_coordinates_factory.Bxy(); // Advection
       }
     }
 
@@ -1795,7 +1797,7 @@ class ELMpb : public PhysicsModel {
     mesh->communicate(Jrhs, ddt(P));
 
     Field3D U1 = ddt(U);
-    U1 += (gamma * tokamak_coordinates_factory.get_Bxy() * tokamak_coordinates_factory.get_Bxy()) * Grad_par(Jrhs, CELL_CENTRE) + (gamma * b0xcv) * Grad(P);
+    U1 += (gamma * tokamak_coordinates_factory.Bxy() * tokamak_coordinates_factory.Bxy()) * Grad_par(Jrhs, CELL_CENTRE) + (gamma * b0xcv) * Grad(P);
 
     // Second matrix, solving Alfven wave dynamics
     static std::unique_ptr<InvertPar> invU{nullptr};
@@ -1804,7 +1806,7 @@ class ELMpb : public PhysicsModel {
     }
 
     invU->setCoefA(1.);
-    invU->setCoefB(-SQ(gamma) * tokamak_coordinates_factory.get_Bxy() * tokamak_coordinates_factory.get_Bxy());
+    invU->setCoefB(-SQ(gamma) * tokamak_coordinates_factory.Bxy() * tokamak_coordinates_factory.Bxy());
     ddt(U) = invU->solve(U1);
     ddt(U).applyBoundary();
 
@@ -1812,9 +1814,9 @@ class ELMpb : public PhysicsModel {
     Field3D phi3 = phiSolver->solve(ddt(U));
     mesh->communicate(phi3);
     phi3.applyBoundary("neumann");
-    Field3D B0phi3 = tokamak_coordinates_factory.get_Bxy() * phi3;
+    Field3D B0phi3 = tokamak_coordinates_factory.Bxy() * phi3;
     mesh->communicate(B0phi3);
-    ddt(Psi) = ddt(Psi) - gamma * Grad_par(B0phi3, loc) / tokamak_coordinates_factory.get_Bxy();
+    ddt(Psi) = ddt(Psi) - gamma * Grad_par(B0phi3, loc) / tokamak_coordinates_factory.Bxy();
     ddt(Psi).applyBoundary();
 
     return 0;
@@ -1846,14 +1848,14 @@ class ELMpb : public PhysicsModel {
     Field3D JP = -b0xGrad_dot_Grad(phi, P0);
     JP.setBoundary("P");
     JP.applyBoundary();
-    Field3D B0phi = tokamak_coordinates_factory.get_Bxy() * phi;
+    Field3D B0phi = tokamak_coordinates_factory.Bxy() * phi;
     mesh->communicate(B0phi);
-    Field3D JPsi = -Grad_par(B0phi, loc) / tokamak_coordinates_factory.get_Bxy();
+    Field3D JPsi = -Grad_par(B0phi, loc) / tokamak_coordinates_factory.Bxy();
     JPsi.setBoundary("Psi");
     JPsi.applyBoundary();
 
-    Field3D JU = b0xcv * Grad(ddt(P)) - SQ(tokamak_coordinates_factory.get_Bxy()) * Grad_par(Jpar, CELL_CENTRE)
-                 + SQ(tokamak_coordinates_factory.get_Bxy()) * b0xGrad_dot_Grad(ddt(Psi), J0, CELL_CENTRE);
+    Field3D JU = b0xcv * Grad(ddt(P)) - SQ(tokamak_coordinates_factory.Bxy()) * Grad_par(Jpar, CELL_CENTRE)
+                 + SQ(tokamak_coordinates_factory.Bxy()) * b0xGrad_dot_Grad(ddt(Psi), J0, CELL_CENTRE);
     JU.setBoundary("U");
     JU.applyBoundary();
 
diff --git a/examples/gyro-gem/gem.cxx b/examples/gyro-gem/gem.cxx
index eee64f8473..be92ddea0a 100644
--- a/examples/gyro-gem/gem.cxx
+++ b/examples/gyro-gem/gem.cxx
@@ -246,7 +246,7 @@ class GEM : public PhysicsModel {
 
     if (mesh->get(Bbar, "Bbar")) {
       if (mesh->get(Bbar, "bmag")) {
-        Bbar = max(tokamak_coordinates_factory.get_Bxy(), true);
+        Bbar = max(tokamak_coordinates_factory.Bxy(), true);
       }
     }
     Bbar = options["Bbar"].withDefault(Bbar); // Override in options file
@@ -349,7 +349,7 @@ class GEM : public PhysicsModel {
 
     B0vec.covariant = false;
     B0vec.x = 0.;
-    B0vec.y = tokamak_coordinates_factory.get_Bpxy() / tokamak_coordinates_factory.get_hthe();
+    B0vec.y = tokamak_coordinates_factory.Bpxy() / tokamak_coordinates_factory.hthe();
     B0vec.z = 0.;
 
     // Precompute this for use in RHS
diff --git a/include/bout/tokamak_coordinates_factory.hxx b/include/bout/tokamak_coordinates_factory.hxx
index def2be950e..5690d4fdab 100644
--- a/include/bout/tokamak_coordinates_factory.hxx
+++ b/include/bout/tokamak_coordinates_factory.hxx
@@ -90,12 +90,12 @@ public:
     return coord;
   }
 
-  const Field2D& get_Rxy() const { return Rxy_m; }
-  const Field2D& get_Bpxy() const { return Bpxy_m; }
-  const Field2D& get_Btxy() const { return Btxy_m; }
-  const Field2D& get_Bxy() const { return Bxy_m; }
-  const Field2D& get_hthe() const { return hthe_m; }
-  const FieldMetric& get_ShearFactor() const { return ShearFactor_m; }
+  const Field2D& Rxy() const { return Rxy_m; }
+  const Field2D& Bpxy() const { return Bpxy_m; }
+  const Field2D& Btxy() const { return Btxy_m; }
+  const Field2D& Bxy() const { return Bxy_m; }
+  const Field2D& hthe() const { return hthe_m; }
+  const FieldMetric& ShearFactor() const { return ShearFactor_m; }
 };
 
 #endif //BOUT_TOKAMAK_COORDINATES_FACTORY_HXX
diff --git a/tests/integrated/test-drift-instability/2fluid.cxx b/tests/integrated/test-drift-instability/2fluid.cxx
index 022ae55314..23dbb7e837 100644
--- a/tests/integrated/test-drift-instability/2fluid.cxx
+++ b/tests/integrated/test-drift-instability/2fluid.cxx
@@ -134,7 +134,7 @@ class TwoFluid : public PhysicsModel {
     const bool ShiftXderivs = (*globalOptions)["ShiftXderivs"].withDefault(false);
     if (ShiftXderivs) {
       ShearFactor = 0.0; // I disappears from metric
-      b0xcv.z += tokamak_coordinates_factory.get_ShearFactor() * b0xcv.x;
+      b0xcv.z += tokamak_coordinates_factory.ShearFactor() * b0xcv.x;
       noshear = true;
     }
 
diff --git a/tests/integrated/test-interchange-instability/2fluid.cxx b/tests/integrated/test-interchange-instability/2fluid.cxx
index 385d98f235..843797e365 100644
--- a/tests/integrated/test-interchange-instability/2fluid.cxx
+++ b/tests/integrated/test-interchange-instability/2fluid.cxx
@@ -104,7 +104,7 @@ class Interchange : public PhysicsModel {
                                                                  bmag / 1e4, ShearFactor);
 
     if (ShiftXderivs) {
-      b0xcv.z += tokamak_coordinates_factory.get_ShearFactor() * b0xcv.x;
+      b0xcv.z += tokamak_coordinates_factory.ShearFactor() * b0xcv.x;
     }
 
     /************** NORMALISE QUANTITIES *****************/