diff --git a/manual/sphinx/user_docs/coordinates.rst b/manual/sphinx/user_docs/coordinates.rst
index 19ebeca5b5..f5aa37962d 100644
--- a/manual/sphinx/user_docs/coordinates.rst
+++ b/manual/sphinx/user_docs/coordinates.rst
@@ -20,36 +20,39 @@ BOUT++ tokamak models, and useful derivations and expressions.
 Orthogonal toroidal coordinates
 ===============================
 
-Starting with an orthogonal toroidal coordinate system
-`\left(\psi, \theta, \zeta\right)`, where `\psi` is the
-poloidal flux, `\theta` the poloidal angle (from `0` to
-`2\pi`), and `\zeta` the toroidal angle (also `0` to
-`2\pi`). We have that the magnetic field
-`{\boldsymbol{B}}` can be expressed as
+Starting with an orthogonal, right-handed toroidal coordinate system
+`\left(r, \theta, \phi\right)`. `\theta` is the poloidal angle (from
+`0` to `2\pi`) in the clockwise direction in the right R-Z
+plane. `\phi` is the toroidal angle (also `0` to `2\pi`) going
+anti-clockwise from the top of the tokamak.
 
-.. math::
+We define the poloidal magnetic field `B_{pol}` as the component of
+the magnetic field in the `\theta` direction, and the toroidal field `B_\text{tor}`
+as the component of the magnetic field in the `\phi` direction.
 
-   \begin{aligned}
-    {\boldsymbol{B}}=& B_\theta \nabla \theta + B_\zeta \nabla \zeta\\ =& B_\theta
-       {\boldsymbol{e}}_\theta + B_\zeta {\boldsymbol{e}}_\zeta\\ =& {B_{\text{pol}}}h_\theta {\boldsymbol{e}}_\theta + {B_{\text{tor}}}
-       R {\boldsymbol{e}}_\zeta\\ =& {B_{\text{pol}}}\hat{{\boldsymbol{e}}}_\theta + {B_{\text{tor}}}\hat{{\boldsymbol{e}}}_\zeta\end{aligned}
+We now introduce the poloidal flux `\psi` as the new radial
+coordinate.  If the poloidal magnetic field `B_\text{pol}` is positive
+then `\psi` increases with radius; if `B_\text{pol}` is negative then
+`\psi` decreases with radius. To keep the coordinate system
+right-handed, we define a new toroidal coordinate `\zeta` which is
+defined as `\zeta = \sigma_{B\text{pol}}\phi`, where the sign of the
+poloidal magnetic field is `\sigma_{B\text{pol}} \equiv {B_{\text{pol}}}/
+\left|{B_{\text{pol}}}\right|`. If `B_\text{pol} > 0` then `\zeta` is
+anti-clockwise looking down from above the tokamak, and if `B_\text{pol} <
+0` then `\zeta` is clockwise. This coordinate system `\left(\psi,
+\theta, \zeta\right)` is orthogonal and right-handed.
 
-The magnitudes of the unit vectors are
+The magnitudes of the basis vectors are
 
 .. math::
 
    \begin{aligned}
-   \left|\hat{{\boldsymbol{e}}}_\psi\right| = \frac{1}{R\left|{B_{\text{pol}}}\right|} \qquad \left|\hat{{\boldsymbol{e}}}_\theta\right| = {h_\theta}
-   \qquad \left|\hat{{\boldsymbol{e}}}_\zeta\right| = R\end{aligned}
+   \left|{\boldsymbol{e}}_\psi\right| = \frac{1}{R\left|{B_{\text{pol}}}\right|} \qquad \left|\boldsymbol{e}_\theta\right| = {h_\theta}
+   \qquad \left|\boldsymbol{e}_\zeta\right| = R
+   \end{aligned}
 
-where `{h_\theta}` is the poloidal arc length per radian. The
-coordinate system is right handed, so
-`\hat{{\boldsymbol{e}}}_\psi\times\hat{{\boldsymbol{e}}}_\theta
-= \hat{{\boldsymbol{e}}}_\zeta`,
-`\hat{{\boldsymbol{e}}}_\psi\times\hat{{\boldsymbol{e}}}_\zeta =
--\hat{{\boldsymbol{e}}}_\theta` and
-`\hat{{\boldsymbol{e}}}_\theta\times\hat{{\boldsymbol{e}}}_\zeta
-= \hat{{\boldsymbol{e}}}_\psi`.  The covariant metric coefficients are
+where `{h_\theta}` is the poloidal arc length per radian.
+The non-zero covariant metric coefficients are
 
 .. math::
 
@@ -65,69 +68,87 @@ and the magnitudes of the reciprocal vectors are therefore
    \left|\nabla\psi\right| = R\left|{B_{\text{pol}}}\right| \qquad \left|\nabla\theta\right| = \frac{1}{h_\theta}
    \qquad \left|\nabla\zeta\right| = \frac{1}{R}\end{aligned}
 
-Because the coordinate system is orthogonal, `g^{ii} = 1/g_{ii}`
-and so the cross-products can be calculated as
+The cross products are:
+
+.. math::
+
+   \boldsymbol{e}_\psi\times\boldsymbol{e}_\theta = J_{\psi\theta\zeta} \nabla\zeta \qquad 
+   \boldsymbol{e}_\psi\times\boldsymbol{e}_\zeta = -J_{\psi\theta\zeta} \nabla\theta \qquad
+   \boldsymbol{e}_\theta\times\boldsymbol{e}_\zeta = J_{\psi\theta\zeta} \nabla\psi
+
+where `J_{\psi\theta\zeta} = h_\theta / \left|{B_{\text{pol}}}\right|` is the Jacobian, which is
+always positive. Similarly,
 
 .. math::
 
    \begin{aligned}
-   \nabla\psi\times\nabla\theta = &\hat{{\boldsymbol{e}}}^\psi\times \hat{{\boldsymbol{e}}}^\theta =
-       g^{\psi\psi}{\boldsymbol{e}}_\psi\times g^{\theta\theta}{\boldsymbol{e}}_\theta \nonumber \\ =
-       & g^{\psi\psi}g^{\theta\theta}h_\psi h_\theta
-       \hat{e}_\psi\times\hat{e}_\theta \nonumber \\ = &\frac{1}{h_\psi
-   h_\theta}\hat{{\boldsymbol{e}}}_\zeta = \frac{R\left|{B_{\text{pol}}}\right|}{h_\theta}\hat{e}_\zeta\end{aligned}
+   \nabla\psi \times \nabla\theta = \frac{1}{J_{\psi\theta\zeta}} \boldsymbol{e}_\zeta \qquad
+   \nabla\psi \times \nabla\zeta = - \frac{1}{J_{\psi\theta\zeta}} \boldsymbol{e}_\theta \qquad
+   \nabla\theta \times \nabla\zeta = \frac{1}{J_{\psi\theta\zeta}} \boldsymbol{e}_\psi
+   \end{aligned}
 
-Similarly,
+The magnetic field `{\boldsymbol{B}}` can be expressed as
 
 .. math::
 
    \begin{aligned}
-   \nabla\psi\times\nabla\zeta = -\left|{B_{\text{pol}}}\right|\hat{{\boldsymbol{e}}}_\theta \qquad
-   \nabla\theta\times\nabla\zeta = \frac{1}{Rh_\theta}\hat{{\boldsymbol{e}}}_\psi =
-   \frac{1}{h_\theta R^2\left|{B_{\text{pol}}}\right|}\nabla \psi\end{aligned}
+    {\boldsymbol{B}}=& B_{\text{pol}} \frac{\boldsymbol{e}_\theta}{h_\theta} + B_{\text{tor}} \frac{\boldsymbol{e}_\phi}{R} \\
+    =& {B_{\text{pol}}}\hat{{\boldsymbol{e}}}_\theta + {B_{\text{tor}}}\hat{{\boldsymbol{e}}}_\phi\end{aligned}
+
+where the hats on the basis vectors indicate unit directions e.g. `\hat{{\boldsymbol{e}}}_\theta = {\boldsymbol{e}}_\theta / \left|{\boldsymbol{e}}_\theta\right|`.
 
 Field-aligned coordinates
 =========================
 
 In order to efficiently simulate (predominantly) field-aligned
-structures, grid-points are placed in a field-aligned coordinate system.
-We define
-`\sigma_{B\theta} \equiv {B_{\text{pol}}}/ \left|{B_{\text{pol}}}\right|`
-i.e. the sign of the poloidal field. The new coordinates
-`\left(x,y,z\right)` are defined by:
+structures, the standard coordinate system used by BOUT++ models is a
+Clebsch system where grid-points are aligned to the magnetic field
+along the `y` coordinate.
+
+To align to the magnetic field we define a local field line pitch `\nu`:
 
 .. math::
-   :label: eq:coordtransform
 
    \begin{aligned}
-   x = {\sigma_{B\theta}}\left(\psi - \psi_0\right) \qquad y = \theta \qquad z = \sigma_{B\theta}
-   \left(\zeta - \int_{\theta_0}^{\theta}\nu\left(\psi,\theta\right)d\theta\right)
+   \nu\left(\psi, \theta\right) = \frac{{\boldsymbol{B}}\cdot\nabla\phi}{{\boldsymbol{B}}\cdot\nabla\theta} =
+   \frac{{B_{\text{tor}}}{h_\theta}}{{B_{\text{pol}}}R}
    \end{aligned}
+The sign of the poloidal field `{B_{\text{pol}}}` and toroidal field 
+`{B_{\text{tor}}}` can be either + or -.
 
-Where `\nu` is the local field-line pitch given by
+The field-aligned coordinates `\left(x,y,z\right)` are defined by:
 
 .. math::
+   :label: eq:coordtransform
 
    \begin{aligned}
-   \nu\left(\psi, \theta\right) = \frac{{\boldsymbol{B}}\cdot\nabla\zeta}{{\boldsymbol{B}}\cdot\nabla\theta} =
-   \frac{{B_{\text{tor}}}{h_\theta}}{{B_{\text{pol}}}R} = \frac{\left(F/R\right)h_\theta}{{B_{\text{pol}}}R} = FJ/R^2\end{aligned}
-
-where `F={B_{\text{tor}}}R` is a function only of `\psi`
-(sometimes called the poloidal current function).
+   x = {\sigma_{B\text{pol}}}\left(\psi - \psi_0\right) \qquad y = \theta \qquad z = \sigma_{B\text{pol}}
+   \left(\phi - \int_{\theta_0}^{\theta}\nu\left(\psi,\theta\right)d\theta\right)
+   \end{aligned}
 
 The coordinate system is chosen so that `x` increases radially
-outwards, from plasma to the wall. The sign of the toroidal field
-`{B_{\text{tor}}}` can then be either + or -.
+outwards, from plasma to the wall. The `y` coordinate increases in the
+same direction as `\theta` i.e. clockwise in the right-hand poloidal
+plane. The `z` coordinate increases in the same direction as `\zeta`
+i.e.  anti-clockwise looking from the top if `B_{pol}>0` and clockwise
+if `B_{pol} < 0`.
+
+This coordinate system is right-handed if `B_{pol}>0`, and left-handed if `B_{pol}<0`.
+The Jacobian of this coordinate system, `J_{xyz} = {h_\theta} / {B_{\text{pol}}}`, can
+therefore be positive or negative. This therefore differs from the Jacobian for the
+orthogonal system above: `J_{xyz} = \sigma_{B\text{pol}} J_{\psi\theta\zeta}`.
 
-The contravariant basis vectors are therefore
+The reciprocal basis vectors are
 
 .. math::
 
    \begin{aligned}
-   \nabla x = {\sigma_{B\theta}}\nabla \psi \qquad \nabla y = \nabla \theta \qquad \nabla z =
-   {\sigma_{B\theta}}\left(\nabla\zeta - \left[\int_{\theta_0}^\theta{\frac{\partial \nu\left(\psi,
-   \theta\right)}{\partial \psi}} d\theta\right] \nabla\psi - \nu\left(\psi, \theta\right)\nabla\theta\right)\end{aligned}
-
+   \nabla x = {\sigma_{B\text{pol}}}\nabla \psi \qquad
+   \nabla y = \nabla \theta \qquad
+   \nabla z = \nabla\zeta - \sigma_{B\text{pol}}\left[\int_{\theta_0}^\theta{\frac{\partial \nu\left(\psi,\theta\right)}{\partial \psi}} d\theta\right] \nabla\psi
+   - \sigma_{B\text{pol}}\nu\left(\psi, \theta\right)\nabla\theta
+   \end{aligned}
+  
 The term in square brackets is the integrated local shear:
 
 .. math::
@@ -135,6 +156,24 @@ The term in square brackets is the integrated local shear:
    \begin{aligned}
    I = \int_{y_0}^y\frac{\partial\nu\left(x, y\right)}{\partial\psi}dy\end{aligned}
 
+  
+The basis vectors are:
+
+.. math::
+   
+   \begin{aligned}
+   \boldsymbol{e}_x =& J_{xyz}\left(\nabla y \times \nabla z\right) = {\sigma_{B\text{pol}}} {\boldsymbol{e}}_\psi + I{\boldsymbol{e}}_\zeta \\
+   \boldsymbol{e}_y =& J_{xyz}\left(\nabla z \times \nabla x\right) = {\boldsymbol{e}}_\theta + \nu{\boldsymbol{e}}_\phi \\
+   \boldsymbol{e}_z =& J_{xyz}\left(\nabla x \times \nabla y\right) = {\boldsymbol{e}}_\zeta
+   \end{aligned}
+ 
+where `{\boldsymbol{e}}_\phi =
+{\sigma_{B\text{pol}}}{\boldsymbol{e}}_\zeta` is always anticlockwise when
+seen from above the tokamak looking down. The direction of
+`{\boldsymbol{e}}_\zeta` depends on the sign of the poloidal field
+`\sigma_{B\text{pol}}`. Note that `J_{xyz} = \sigma_{B\text{pol}} J_{\psi\theta\zeta}`, and
+can be either positive or negative.
+
 Magnetic field
 --------------
 
@@ -143,7 +182,7 @@ Magnetic field is given in Clebsch form by:
 .. math::
 
    \begin{aligned}
-   {\boldsymbol{B}}= \nabla z\times \nabla x = \frac{1}{J}{\boldsymbol{e}}_y\end{aligned}
+   {\boldsymbol{B}}= \nabla z\times \nabla x = \frac{1}{J_{xyz}}{\boldsymbol{e}}_y\end{aligned}
 
 The contravariant components of this are then
 
@@ -164,7 +203,7 @@ and the covariant components calculated using `g_{ij}` as
 .. math::
 
    \begin{aligned}
-   B_x = {\sigma_{B\theta}}{B_{\text{tor}}}I R \qquad B_y = \frac{B^2 {h_\theta}}{{B_{\text{pol}}}} \qquad B_z = {\sigma_{B\theta}}{B_{\text{tor}}}R\end{aligned}
+   B_x = {\sigma_{B\text{pol}}}{B_{\text{tor}}}I R \qquad B_y = \frac{B^2 {h_\theta}}{{B_{\text{pol}}}} \qquad B_z = {\sigma_{B\text{pol}}}{B_{\text{tor}}}R\end{aligned}
 
 The unit vector in the direction of equilibrium `{\boldsymbol{B}}` is
 therefore
@@ -172,7 +211,7 @@ therefore
 .. math::
 
    \begin{aligned}
-   {\boldsymbol{b}} = \frac{1}{JB}{\boldsymbol{e}}_y = \frac{1}{JB}\left[g_{xy}\nabla x + g_{yy}\nabla y
+   {\boldsymbol{b}} = \frac{1}{J_{xyz}B}{\boldsymbol{e}}_y = \frac{1}{J_{xyz}B}\left[g_{xy}\nabla x + g_{yy}\nabla y
    + g_{yz}\nabla z\right]\end{aligned}
 
 Jacobian and metric tensors
@@ -183,7 +222,7 @@ The Jacobian of this coordinate system is
 .. math::
 
    \begin{aligned}
-   J^{-1} \equiv \left(\nabla x\times\nabla y\right)\cdot\nabla z = {B_{\text{pol}}}/ {h_\theta}\end{aligned}
+   J_{xyz}^{-1} \equiv \left(\nabla x\times\nabla y\right)\cdot\nabla z = {B_{\text{pol}}}/ {h_\theta}\end{aligned}
 
 which can be either positive or negative, depending on the sign of
 `{B_{\text{pol}}}`. The contravariant metric tensor is
@@ -195,8 +234,8 @@ given by:
    g^{ij} \equiv {\boldsymbol{e}}^i \cdot{\boldsymbol{e}}^j \equiv \nabla u^i \cdot \nabla u^j = \left(%
    \begin{array}{ccc}
    \left(R{B_{\text{pol}}}\right)^2 & 0 & -I\left(R{B_{\text{pol}}}\right)^2 \\
-   0 & 1 / {h_\theta}^2 & -{\sigma_{B\theta}}\nu / {h_\theta}^2 \\
-   -I\left(R{B_{\text{pol}}}\right)^2 & -{\sigma_{B\theta}}\nu / {h_\theta}^2 & I^2\left(R{B_{\text{pol}}}\right)^2 + B^2 /
+   0 & 1 / {h_\theta}^2 & -{\sigma_{B\text{pol}}}\nu / {h_\theta}^2 \\
+   -I\left(R{B_{\text{pol}}}\right)^2 & -{\sigma_{B\text{pol}}}\nu / {h_\theta}^2 & I^2\left(R{B_{\text{pol}}}\right)^2 + B^2 /
    \left(R{B_{\text{pol}}}\right)^2
    \end{array}
    %
@@ -209,15 +248,185 @@ and the covariant metric tensor:
    \begin{aligned}
    g_{ij} \equiv {\boldsymbol{e}}_i \cdot{\boldsymbol{e}}_j = \left(%
    \begin{array}{ccc}
-   I^2 R^2 + 1 / {\left({R{B_{\text{pol}}}}\right)^2}& {\sigma_{B\theta}}{B_{\text{tor}}}{h_\theta}I R / {B_{\text{pol}}}& I R^2 \\
-   {\sigma_{B\theta}}{B_{\text{tor}}}{h_\theta}I R / {B_{\text{pol}}}& B^2{h_\theta}^2 / {B_{\text{pol}}}^2 & {\sigma_{B\theta}}{B_{\text{tor}}}{h_\theta}R / {B_{\text{pol}}}\\
-   I R^2 & {\sigma_{B\theta}}{B_{\text{tor}}}{h_\theta}R / {B_{\text{pol}}}& R^2
+   I^2 R^2 + 1 / {\left({R{B_{\text{pol}}}}\right)^2}& {\sigma_{B\text{pol}}}{B_{\text{tor}}}{h_\theta}I R / {B_{\text{pol}}}& I R^2 \\
+   {\sigma_{B\text{pol}}}{B_{\text{tor}}}{h_\theta}I R / {B_{\text{pol}}}& B^2{h_\theta}^2 / {B_{\text{pol}}}^2 & {\sigma_{B\text{pol}}}{B_{\text{tor}}}{h_\theta}R / {B_{\text{pol}}}\\
+   I R^2 & {\sigma_{B\text{pol}}}{B_{\text{tor}}}{h_\theta}R / {B_{\text{pol}}}& R^2
    \end{array}
    %
     \right)\end{aligned}
 
-Differential operators
-----------------------
+or equivalently:
+
+.. math::
+
+   \begin{aligned}
+   g_{ij} = \left(%
+   \begin{array}{ccc}
+   I^2 R^2 + 1 / {\left({R{B_{\text{pol}}}}\right)^2}& {\sigma_{B\text{pol}}} I \nu R^2 & I R^2 \\
+   {\sigma_{B\text{pol}}} I \nu R^2 & J_{xyz}^2B^2 & {\sigma_{B\text{pol}}} \nu R^2 \\
+   I R^2 & {\sigma_{B\text{pol}}}\nu R^2 & R^2
+   \end{array}
+   %
+   \right)\end{aligned}
+
+zShift
+------
+
+The `\texttt{zShift}` is used to connect grid cells along the magnetic
+field. It is the `z` angle of a point on a field line relative to a
+reference location:
+
+.. math::
+
+   \begin{aligned}
+   \texttt{zShift}\left(x, y\right) &= \int_{y = 0}^{y}\frac{{\boldsymbol{B}}\cdot\nabla z}{{\boldsymbol{B}}\cdot\nabla y} dy \\
+   &= \int_{\theta = 0}^{\theta}\frac{{\sigma_{B\text{pol}}}{B_{\text{tor}}}{h_\theta}}{{B_{\text{pol}}}R} d\theta \\
+   &= {\sigma_{B\text{pol}}} \int_{\theta = 0}^{\theta} \nu d\theta
+   \end{aligned}
+
+The `\texttt{ShiftAngle}` is then defined as the change in
+`\texttt{zShift}` between `y=0` and `y=2\pi`: It is the change in the
+`z` coordinate after one poloidal circuit in `y`.
+
+Note that `\texttt{zShift}` can be related to the integrated shear `I`:
+
+.. math::
+
+   \begin{aligned}
+   I = \int_{y_0}^y\frac{\partial\nu\left(x, y\right)}{\partial\psi}dy = \frac{\partial}{\partial x} \texttt{zShift}
+   \end{aligned}
+
+
+Right-handed field-aligned coordinates
+======================================
+
+If the poloidal magnetic field is negative, i.e. anti-clockwise in the right-hand R-Z plane, then the above
+coordinate system is left-handed and the Jacobian `J_{xyz}` is negative.
+To obtain a consistently right-handed coordinate system, we have to reverse the direction of the `y` coordinate
+when `B_{pol} < 0`:
+
+This `\left(x,\eta,z\right)` coordinate system is defined by:
+
+.. math::
+   :label: eq:coordtransform2
+
+   \begin{aligned}
+   x = {\sigma_{B\text{pol}}}\left(\psi - \psi_0\right) \qquad \eta = {\sigma_{B\text{pol}}}\theta \qquad z = \sigma_{B\text{pol}}
+   \left(\phi - \int_{\theta_0}^{\theta}\nu\left(\psi,\theta\right)d\theta\right)
+   \end{aligned}
+
+The radial coordinate `x` always points outwards. The `\eta` coordinate
+increases in the direction of the poloidal magnetic field: clockwise
+in the right-hand poloidal plane if `B_{pol} > 0`, and anti-clockwise
+otherwise.  The `z` coordinate increases in the same direction as
+`\zeta` i.e.  anti-clockwise looking from the top if `B_{pol}>0` and
+clockwise if `B_{pol} < 0`.
+
+This is still a Clebsch coordinate system:
+
+.. math::
+
+   \begin{aligned}
+   {\boldsymbol{B}}= \nabla z\times \nabla x = \frac{1}{J_{x\eta z}}{\boldsymbol{e}}_\eta
+   \end{aligned}
+
+but the Jacobian is now always positive:
+
+.. math::
+
+   \begin{aligned}
+   J_{x\eta z} = h_\theta / \left|B_{\text{pol}}\right|
+   \end{aligned}
+
+
+The reciprocal basis vectors are
+
+.. math::
+   
+   \begin{aligned}
+   \nabla x =& {\sigma_{B\text{pol}}} \nabla \psi \\
+   \nabla \eta =& {\sigma_{B\text{pol}}} \nabla \theta \\
+   \nabla z =& \nabla \zeta - {\sigma_{B\text{pol}}} I \nabla \psi - {\sigma_{B\text{pol}}}\nu\nabla\theta
+   \end{aligned}
+
+and basis vectors
+
+.. math::
+   
+   \begin{aligned}
+   \boldsymbol{e}_x =& J_{x\eta z}\left(\nabla y \times \nabla z\right) = {\sigma_{B\text{pol}}} {\boldsymbol{e}}_\psi + I{\boldsymbol{e}}_\zeta \\
+   \boldsymbol{e}_\eta =& J_{x\eta z}\left(\nabla z \times \nabla x\right) = {\sigma_{B\text{pol}}} {\boldsymbol{e}}_\theta + \nu{\boldsymbol{e}}_\zeta \\
+   \boldsymbol{e}_z =& J_{x\eta z}\left(\nabla x \times \nabla y\right) = {\boldsymbol{e}}_\zeta
+   \end{aligned}
+
+
+The contravariant metric tensor is:
+
+.. math::
+
+   \begin{aligned}
+   g^{ij} \equiv {\boldsymbol{e}}^i \cdot{\boldsymbol{e}}^j \equiv \nabla u^i \cdot \nabla u^j = \left(%
+   \begin{array}{ccc}
+   \left(R{B_{\text{pol}}}\right)^2 & 0 & -I\left(R{B_{\text{pol}}}\right)^2 \\
+   0 & 1 / {h_\theta}^2 & -\nu / {h_\theta}^2 \\
+   -I\left(R{B_{\text{pol}}}\right)^2 & -\nu / {h_\theta}^2 & I^2\left(R{B_{\text{pol}}}\right)^2 + B^2 /
+   \left(R{B_{\text{pol}}}\right)^2
+   \end{array}
+   %
+   \right)\end{aligned}
+
+and the covariant metric tensor:
+
+.. math::
+
+   \begin{aligned}
+   g_{ij} = \left(%
+   \begin{array}{ccc}
+   I^2 R^2 + 1 / {\left({R{B_{\text{pol}}}}\right)^2}& I \nu R^2 & I R^2 \\
+   I \nu R^2 & J_{x\eta z}^2B^2 & \nu R^2 \\
+   I R^2 & \nu R^2 & R^2
+   \end{array}
+   %
+   \right)\end{aligned}
+
+The `\texttt{zShift}` quantity is the `z` angle of a point on a field
+line relative to a reference location. This is a scalar which doesn't
+change if the sign of the `\eta` coordinate is reversed:
+
+.. math::
+
+   \begin{aligned}
+   \texttt{zShift}\left(x, \eta\right) = \int_{\eta = 0}^{\eta}\frac{{\boldsymbol{B}}\cdot\nabla z}{{\boldsymbol{B}}\cdot\nabla \eta} d\eta =
+   \int_{\theta = 0}^{{\sigma_{B\text{pol}}}\theta}\frac{{\sigma_{B\text{pol}}}{B_{\text{tor}}}{h_\theta}}{{B_{\text{pol}}}R} d\theta
+   \end{aligned}
+
+The `\texttt{ShiftAngle}` quantity is related to `\texttt{zShift}`: It
+is the change in `\texttt{zShift}` from `\eta=0` to `\eta=2\pi`. It
+therefore does change sign if the `\eta` direction is reversed.
+
+The differences from the previous `\left(x,y,z\right)` coordinate
+system are that `g_{xy}`, `g_{yz}`, `g^{yz}`, `J` and
+`\texttt{ShiftAngle}` are multiplied by `{\sigma_{B\text{pol}}}` to obtain
+their equivalents in the `\left(x,\eta,z\right)` coordinate system. If
+`B_{pol} < 0` so the poloidal magnetic field is anticlockwise in the
+right-hand R-Z plane, then the `\eta` direction changes.
+
+
+Differential operators in field-aligned coordinates
+===================================================
+
+These operators are valid for either `\left(x,y,z\right)` or
+`\left(x,\eta,z\right)` field-aligned coordinates defined
+above. Unless explicitly stated, in the sections that follow `y` will
+be used to indicate the parallel coordinate (`y` or `\eta`). In a few
+places the sign of `B_\text{pol}` may appear, depending on whether `y`
+or `\eta` is used for the parallel coordinate, so we define
+
+.. math::
+    \sigma_y = \begin{cases}
+        \sigma_{B\text{pol}} & \text{if using }(x,y,z) \\
+        +1 & \text{if using }(x,\eta,z)
+    \end{cases}
 
 The derivative of a scalar field `f` along the *unperturbed*
 magnetic field `{\boldsymbol{b}}_0` is given by
@@ -226,9 +435,9 @@ magnetic field `{\boldsymbol{b}}_0` is given by
 
    \begin{aligned}
    \partial^0_{||}f \equiv {\boldsymbol{b}}_0 \cdot\nabla f =
-   \frac{1}{\sqrt{g_{yy}}}{\frac{\partial f}{\partial y}} = \frac{{B_{\text{pol}}}}{B{h_\theta}}{\frac{\partial f}{\partial y}}\end{aligned}
+   \frac{1}{JB}{\frac{\partial f}{\partial y}} = \frac{\sigma_y|{B_{\text{pol}}}|}{B{h_\theta}}{\frac{\partial f}{\partial y}}\end{aligned}
 
-whilst the parallel divergence is given by
+Note that J could be positive or negative. The parallel divergence is given by
 
 .. math::
 
@@ -245,7 +454,7 @@ the Laplacian operator is given by
        \frac{\partial^2}{\partial y^2}\left|\nabla y\right|^2 +
        \frac{\partial^2}{\partial z^2}\left|\nabla z\right|^2 \nonumber \\
        &-2\frac{\partial^2}{\partial x\partial z}I\left(R{B_{\text{pol}}}\right)^2 -
-       2\frac{\partial^2}{\partial y\partial z}\frac{\nu}{h_\theta^2}\\
+       2\frac{\partial^2}{\partial y\partial z}\frac{\sigma_y\nu}{h_\theta^2}\\
        &+\frac{\partial}{\partial x}\nabla^2x + \frac{\partial}{\partial
    y}\nabla^2y + \frac{\partial}{\partial z}\nabla^2z \nonumber\end{aligned}
 
@@ -263,7 +472,7 @@ Using equation :eq:`eq:laplace_expand` for
 
    \begin{aligned}
    \nabla^2z = -\frac{{B_{\text{pol}}}}{h_\theta}\left[\frac{\partial}{\partial x}\left(IR^2{B_{\text{pol}}}
-   h_\theta\right) + \frac{\partial}{\partial y}\left(\frac{\nu}{{B_{\text{pol}}}h_\theta}\right)\right]\end{aligned}
+   h_\theta\right) + \sigma_y \frac{\partial}{\partial y}\left(\frac{\nu}{{B_{\text{pol}}}h_\theta}\right)\right]\end{aligned}
 
 Neglecting some parallel derivative terms, the perpendicular Laplacian
 can be written:
@@ -281,7 +490,7 @@ The second derivative along the equilibrium field
 
    \begin{aligned}
    \partial^2_{||}\phi = \partial^0_{||}\left(\partial^0_{||}\phi\right) =
-   \frac{1}{\sqrt{g_{yy}}}{\frac{\partial }{\partial y}}\left(\frac{1}{\sqrt{g_{yy}}}\right){\frac{\partial  \phi}{\partial y}}
+   \frac{1}{JB}{\frac{\partial }{\partial y}}\left(\frac{1}{JB}\right){\frac{\partial \phi}{\partial y}}
    + \frac{1}{g_{yy}}\frac{\partial^2\phi}{\partial y^2}\end{aligned}
 
 A common expression (the Poisson bracket in reduced MHD) is (from
@@ -291,7 +500,7 @@ equation :eq:`eq:brackets`)):
 
    \begin{aligned}
    {\boldsymbol{b}}_0\cdot\nabla\phi\times\nabla A =
-   \frac{1}{J\sqrt{g_{yy}}}\left[\left(g_{yy}{\frac{\partial \phi}{\partial z}} -
+   \frac{1}{J^2B}\left[\left(g_{yy}{\frac{\partial \phi}{\partial z}} -
    g_{yz}{\frac{\partial \phi}{\partial y}}\right){\frac{\partial A}{\partial x}} + \left(g_{yz}{\frac{\partial \phi}{\partial x}} -
    g_{xy}{\frac{\partial \phi}{\partial z}}\right){\frac{\partial A}{\partial y}} + \left(g_{xy}{\frac{\partial \phi}{\partial y}} -
    g_{yy}{\frac{\partial \phi}{\partial x}}\right){\frac{\partial A}{\partial z}}\right]\end{aligned}
@@ -315,14 +524,14 @@ Components of the magnetic field in field-aligned coordinates:
 .. math::
 
    \begin{aligned}
-   B^y = \frac{{B_{\text{pol}}}}{{h_\theta}} \qquad B^x = B^z = 0\end{aligned}
+   B^y = \frac{\sigma_y{|B_{\text{pol}}|}}{{h_\theta}} \qquad B^x = B^z = 0\end{aligned}
 
 and
 
 .. math::
 
    \begin{aligned}
-   B_x = {\sigma_{B\theta}}{B_{\text{tor}}}I R \qquad B_y = \frac{B^2{h_\theta}}{{B_{\text{pol}}}} \qquad B_z = {\sigma_{B\theta}}{B_{\text{tor}}}R\end{aligned}
+   B_x = {\sigma_{B\text{pol}}}{B_{\text{tor}}}I R \qquad B_y = \sigma_y\frac{B^2{h_\theta}}{{|B_{\text{pol}}|}} \qquad B_z = {\sigma_{B\text{pol}}}{B_{\text{tor}}}R\end{aligned}
 
 Calculate current `{\boldsymbol{J}}= \frac{1}{\mu}{\nabla\times
 {\boldsymbol{B}} }`
@@ -338,10 +547,10 @@ since `{B_{\text{tor}}}R` is a flux-surface quantity, and
 .. math::
 
    \begin{aligned}
-   \left({\nabla\times {\boldsymbol{B}} }\right)^y =& -{\sigma_{B\theta}}\frac{{B_{\text{pol}}}}{{h_\theta}}{\frac{\partial }{\partial x}}\left({B_{\text{tor}}}R\right) \\
+   \left({\nabla\times {\boldsymbol{B}} }\right)^y =& -{\sigma_y\sigma_{B\text{pol}}}\frac{{B_{\text{pol}}}}{{h_\theta}}{\frac{\partial }{\partial x}}\left({B_{\text{tor}}}R\right) \\
        \left({\nabla\times {\boldsymbol{B}} }\right)^z =&
        \frac{{B_{\text{pol}}}}{{h_\theta}}\left[{\frac{\partial }{\partial x}}\left(\frac{B^2{h_\theta}}{{B_{\text{pol}}}}\right) -
-       {\sigma_{B\theta}}{\frac{\partial }{\partial y}}\left({B_{\text{tor}}}I R\right)\right]\end{aligned}
+       {\sigma_{B\text{pol}}}{\frac{\partial }{\partial y}}\left({B_{\text{tor}}}I R\right)\right]\end{aligned}
 
 The second term can be simplified, again using
 `{B_{\text{tor}}}R` constant on flux-surfaces:
@@ -349,7 +558,7 @@ The second term can be simplified, again using
 .. math::
 
    \begin{aligned}
-   {\frac{\partial }{\partial y}}\left({B_{\text{tor}}}I R\right) = {\sigma_{B\theta}}{B_{\text{tor}}}R{\frac{\partial \nu}{\partial x}} \qquad \nu =
+   {\frac{\partial }{\partial y}}\left({B_{\text{tor}}}I R\right) = {\sigma_{B\text{pol}}}{B_{\text{tor}}}R{\frac{\partial \nu}{\partial x}} \qquad \nu =
    \frac{{h_\theta}{B_{\text{tor}}}}{R{B_{\text{pol}}}}\end{aligned}
 
 From these, calculate covariant components:
@@ -362,8 +571,8 @@ From these, calculate covariant components:
        \frac{IR^2{B_{\text{pol}}}}{{h_\theta}}\left[{\frac{\partial }{\partial x}}\left(\frac{B^2{h_\theta}}{{B_{\text{pol}}}}\right) - {B_{\text{tor}}}
        R{\frac{\partial \nu}{\partial x}}\right] \nonumber\\
    %
-   \left({\nabla\times {\boldsymbol{B}} }\right)_y =& -{\sigma_{B\theta}}\frac{B^2{h_\theta}}{{B_{\text{pol}}}}{\frac{\partial }{\partial x}}\left({B_{\text{tor}}}R\right) +
-       {\sigma_{B\theta}}{B_{\text{tor}}}R\left[{\frac{\partial }{\partial x}}\left(\frac{B^2{h_\theta}}{{B_{\text{pol}}}}\right) - {B_{\text{tor}}}R{\frac{\partial \nu}{\partial x}}\right]
+   \left({\nabla\times {\boldsymbol{B}} }\right)_y =& -{\sigma_y\sigma_{B\text{pol}}}\frac{B^2{h_\theta}}{{B_{\text{pol}}}}{\frac{\partial }{\partial x}}\left({B_{\text{tor}}}R\right) +
+       {\sigma_y\sigma_{B\text{pol}}}{B_{\text{tor}}}R\left[{\frac{\partial }{\partial x}}\left(\frac{B^2{h_\theta}}{{B_{\text{pol}}}}\right) - {B_{\text{tor}}}R{\frac{\partial \nu}{\partial x}}\right]
        \\
    %
    \left({\nabla\times {\boldsymbol{B}} }\right)_z =& -{B_{\text{tor}}}R{\frac{\partial }{\partial x}}\left({B_{\text{tor}}}R\right) +
@@ -452,24 +661,28 @@ Parallel current
 .. math::
 
    \begin{aligned}
-   J_{||} = {\boldsymbol{b}}\cdot{\boldsymbol{J}}\qquad b^y = \frac{{B_{\text{pol}}}}{B{h_\theta}}\end{aligned}
+   J_{||} = {\boldsymbol{b}}\cdot{\boldsymbol{J}}\qquad b^y = \sigma_y\frac{{|B_{\text{pol}}|}}{B{h_\theta}}\end{aligned}
 
 and from equation :eq:`eq:curlb_y`:
 
 .. math::
 
    \begin{aligned}
-   J_y = \frac{{\sigma_{B\theta}}}{\mu_0}\left\{-\frac{B^2{h_\theta}}{{B_{\text{pol}}}}{\frac{\partial }{\partial x}}\left({B_{\text{tor}}}R\right) + {B_{\text{tor}}}
-   R\left[{\frac{\partial }{\partial x}}\left(\frac{B^2{h_\theta}}{{B_{\text{pol}}}}\right) - {\sigma_{B\theta}}{B_{\text{tor}}}R{\frac{\partial \nu}{\partial x}}\right]\right\}\end{aligned}
+   J_y = \frac{{\sigma_y\sigma_{B\text{pol}}}}{\mu_0}\left\{-\frac{B^2{h_\theta}}{{B_{\text{pol}}}}{\frac{\partial }{\partial x}}\left({B_{\text{tor}}}R\right) + {B_{\text{tor}}}
+   R\left[{\frac{\partial }{\partial x}}\left(\frac{B^2{h_\theta}}{{B_{\text{pol}}}}\right) - {B_{\text{tor}}}R{\frac{\partial \nu}{\partial x}}\right]\right\}\end{aligned}
 
 since `J_{||} = b^yJ_y`,
 
 .. math::
 
    \begin{aligned}
-   \mu_0 J_{||} ={\sigma_{B\theta}}\frac{{B_{\text{pol}}}{B_{\text{tor}}}
+   \mu_0 J_{||} =\frac{{B_{\text{pol}}}{B_{\text{tor}}}
    R}{B{h_\theta}}\left[{\frac{\partial }{\partial x}}\left(\frac{B^2{h_\theta}}{{B_{\text{pol}}}}\right) - {B_{\text{tor}}}R{\frac{\partial \nu}{\partial x}}\right] -
-   {\sigma_{B\theta}}B{\frac{\partial }{\partial x}}\left({B_{\text{tor}}}R\right)\end{aligned}
+   B{\frac{\partial }{\partial x}}\left({B_{\text{tor}}}R\right)\end{aligned}
+
+Note, this does not depend on our coordinate choices, so does not
+depend on `\sigma_y` or `\sigma_{B\text{pol}}`, as it should not
+since `\mu_0 J_\parallel` is a scalar quantity.
 
 Curvature
 ---------
@@ -487,18 +700,16 @@ Re-arranging, this becomes:
    {\boldsymbol{b}}\times{\boldsymbol{\kappa}} = \nabla\times{\boldsymbol{b}}-
    {\boldsymbol{b}}\left({\boldsymbol{b}}\cdot\left(\nabla\times{\boldsymbol{b}}\right)\right)\end{aligned}
 
-Components of `\nabla\times{\boldsymbol{b}}` are:
+Components of `\nabla\times{\boldsymbol{b}}` are [#curvature]_:
 
 .. math::
 
    \begin{aligned}
-   \left(\nabla\times{\boldsymbol{b}}\right)^x =& {\sigma_{B\theta}}\frac{{B_{\text{pol}}}}{{h_\theta}}{\frac{\partial }{\partial y}}\left(\frac{{B_{\text{tor}}}
-   R}{B}\right) \\ \left(\nabla\times{\boldsymbol{b}}\right)^y =&
-       -{\sigma_{B\theta}}\frac{{B_{\text{pol}}}}{{h_\theta}}{\frac{\partial }{\partial x}}\left(\frac{{B_{\text{tor}}}R}{B}\right) \\
-       \left(\nabla\times{\boldsymbol{b}}\right)^z =&
-       \frac{{B_{\text{pol}}}}{{h_\theta}}{\frac{\partial }{\partial x}}\left(\frac{B{h_\theta}}{{B_{\text{pol}}}}\right) - {\sigma_{B\theta}}\frac{{B_{\text{pol}}}{B_{\text{tor}}}
-       R}{{h_\theta}B}{\frac{\partial \nu}{\partial x}} - {\sigma_{B\theta}}I\frac{{B_{\text{pol}}}}{{h_\theta}}{\frac{\partial }{\partial y}}\left(\frac{{B_{\text{tor}}}
-       R}{B}\right) \\\end{aligned}
+   \left(\nabla\times{\boldsymbol{b}}\right)^x =& {\sigma_y}\frac{{B_{\text{pol}}}}{{h_\theta}}{\frac{\partial }{\partial y}}\left(\frac{{B_{\text{tor}}}
+   R}{B}\right) \\ 
+   \left(\nabla\times{\boldsymbol{b}}\right)^y =& -{\sigma_y}\frac{{B_{\text{pol}}}}{{h_\theta}}{\frac{\partial }{\partial x}}\left(\frac{{B_{\text{tor}}}R}{B}\right) \\
+   \left(\nabla\times{\boldsymbol{b}}\right)^z =& \frac{{B_{\text{pol}}}}{{h_\theta}}{\frac{\partial }{\partial x}}\left(\frac{B{h_\theta}}{{B_{\text{pol}}}}\right) - \frac{{B_{\text{pol}}}{B_{\text{tor}}} R}{{h_\theta}B}{\frac{\partial \nu}{\partial x}} - {\sigma_y}I\frac{{|B_{\text{pol}}|}}{{h_\theta}}{\frac{\partial }{\partial y}}\left(\frac{{B_{\text{tor}}} R}{B}\right) \\
+   \end{aligned}
 
 giving:
 
@@ -507,38 +718,38 @@ giving:
 
    \begin{aligned}
    {\boldsymbol{\kappa}} =& -\frac{{B_{\text{pol}}}}{B h_\theta}\left[{\frac{\partial }{\partial x}}\left(\frac{B
-   h_\theta}{{B_{\text{pol}}}}\right) - {\sigma_{B\theta}}{\frac{\partial }{\partial y}}\left(\frac{{B_{\text{tor}}}I R}{B}\right)\right]\nabla x \nonumber
-   \\ &+ {\sigma_{B\theta}}\frac{{B_{\text{pol}}}}{B h_\theta}{\frac{\partial }{\partial y}}\left(\frac{{B_{\text{tor}}}R}{B}\right)\nabla z
+   h_\theta}{{B_{\text{pol}}}}\right) - {\sigma_y\sigma_{B\text{pol}}}{\frac{\partial }{\partial y}}\left(\frac{{B_{\text{tor}}}I R}{B}\right)\right]\nabla x \nonumber
+   \\ &+ {\sigma_y}\frac{{B_{\text{pol}}}}{B h_\theta}{\frac{\partial }{\partial y}}\left(\frac{{B_{\text{tor}}}R}{B}\right)\nabla z
    \end{aligned}
 
 .. math::
 
    \begin{aligned}
-   {\boldsymbol{b}}\cdot\left(\nabla\times{\boldsymbol{b}}\right) = -{\sigma_{B\theta}}B{\frac{\partial }{\partial x}}\left(\frac{{B_{\text{tor}}}R}{B}\right) +
-   {\sigma_{B\theta}}\frac{{B_{\text{tor}}}{B_{\text{pol}}}R}{B{h_\theta}}{\frac{\partial }{\partial x}}\left(\frac{B{h_\theta}}{{B_{\text{pol}}}}\right) -
-   \frac{{B_{\text{pol}}}{B_{\text{tor}}}^2R^2}{{h_\theta}B^2}{\frac{\partial \nu}{\partial x}}\end{aligned}
+   {\boldsymbol{b}}\cdot\left(\nabla\times{\boldsymbol{b}}\right) = -{\sigma_{B\text{pol}}}B{\frac{\partial }{\partial x}}\left(\frac{{B_{\text{tor}}}R}{B}\right) +
+   {\sigma_{B\text{pol}}}\frac{{B_{\text{tor}}}{B_{\text{pol}}}R}{B{h_\theta}}{\frac{\partial }{\partial x}}\left(\frac{B{h_\theta}}{{B_{\text{pol}}}}\right) -
+   \sigma_{B\text{pol}}\frac{{B_{\text{pol}}}{B_{\text{tor}}}^2R^2}{{h_\theta}B^2}{\frac{\partial \nu}{\partial x}}\end{aligned}
 
 therefore,
 
 .. math::
 
    \begin{aligned}
-   \left({\boldsymbol{b}}\times{\boldsymbol{\kappa}}\right)^x =& {\sigma_{B\theta}}\frac{{B_{\text{pol}}}}{{h_\theta}}{\frac{\partial }{\partial y}}\left(\frac{{B_{\text{tor}}}
-   R}{B}\right) = -{\sigma_{B\theta}}\frac{{B_{\text{pol}}}{B_{\text{tor}}}R}{{h_\theta}B^2}{\frac{\partial B}{\partial y}} \\
-   \left({\boldsymbol{b}}\times{\boldsymbol{\kappa}}\right)^y =& \frac{{B_{\text{pol}}}^2{B_{\text{tor}}}^2
-   R^2}{B^3{h_\theta}^2}{\frac{\partial \nu}{\partial x}} - {\sigma_{B\theta}}\frac{{B_{\text{pol}}}^2{B_{\text{tor}}}
+   \left({\boldsymbol{b}}\times{\boldsymbol{\kappa}}\right)^x =& {\sigma_y}\frac{{B_{\text{pol}}}}{{h_\theta}}{\frac{\partial }{\partial y}}\left(\frac{{B_{\text{tor}}}
+   R}{B}\right) = -{\sigma_y}\frac{{B_{\text{pol}}}{B_{\text{tor}}}R}{{h_\theta}B^2}{\frac{\partial B}{\partial y}} \\
+   \left({\boldsymbol{b}}\times{\boldsymbol{\kappa}}\right)^y =& \sigma_y\frac{{B_{\text{pol}}}^2{B_{\text{tor}}}^2
+   R^2}{B^3{h_\theta}^2}{\frac{\partial \nu}{\partial x}} - {\sigma_y}\frac{{B_{\text{pol}}}^2{B_{\text{tor}}}
    R}{B^2{h_\theta}^2}{\frac{\partial }{\partial x}}\left(\frac{B{h_\theta}}{{B_{\text{pol}}}}\right) \\
    \left({\boldsymbol{b}}\times{\boldsymbol{\kappa}}\right)^z =&
-   \frac{{B_{\text{pol}}}}{{h_\theta}}{\frac{\partial }{\partial x}}\left(\frac{B{h_\theta}}{{B_{\text{pol}}}}\right) - {\sigma_{B\theta}}\frac{{B_{\text{pol}}}{B_{\text{tor}}}
-   R}{{h_\theta}B}{\frac{\partial \nu}{\partial x}} - I\left({\boldsymbol{b}}\times{\boldsymbol{\kappa}}\right)^x\end{aligned}
+   \frac{{B_{\text{pol}}}}{{h_\theta}}{\frac{\partial }{\partial x}}\left(\frac{B{h_\theta}}{{B_{\text{pol}}}}\right) - \frac{{B_{\text{pol}}}{B_{\text{tor}}}
+   R}{{h_\theta}B}{\frac{\partial \nu}{\partial x}} - \sigma_{B\text{pol}} I\left({\boldsymbol{b}}\times{\boldsymbol{\kappa}}\right)^x\end{aligned}
 
 Using equation :eq:`eq:xbalance`:
 
 .. math::
 
    \begin{aligned}
-   B{\frac{\partial }{\partial x}}\left(\frac{B{h_\theta}}{{B_{\text{pol}}}}\right) + \frac{B{h_\theta}}{{B_{\text{pol}}}}{\frac{\partial B}{\partial x}} - {\sigma_{B\theta}}{B_{\text{tor}}}
-   R{\frac{\partial }{\partial x}}\left(\frac{{B_{\text{tor}}}{h_\theta}}{R{B_{\text{pol}}}}\right) + \frac{\mu_0{h_\theta}}{{B_{\text{pol}}}}{\frac{\partial P}{\partial x}} =
+   \sigma_{B\text{pol}}B{\frac{\partial }{\partial x}}\left(\frac{B{h_\theta}}{{B_{\text{pol}}}}\right) + \sigma_{B\text{pol}}\frac{B{h_\theta}}{{B_{\text{pol}}}}{\frac{\partial B}{\partial x}} - {\sigma_{B\text{pol}}}{B_{\text{tor}}}
+   R{\frac{\partial }{\partial x}}\left(\frac{{B_{\text{tor}}}{h_\theta}}{R{B_{\text{pol}}}}\right) + \sigma_{B\text{pol}}\frac{\mu_0{h_\theta}}{{B_{\text{pol}}}}{\frac{\partial P}{\partial x}} =
    0\end{aligned}
 
 we can re-write the above components as:
@@ -546,12 +757,26 @@ we can re-write the above components as:
 .. math::
 
    \begin{aligned}
-   \left({\boldsymbol{b}}\times{\boldsymbol{\kappa}}\right)^y =& {\sigma_{B\theta}}\frac{{B_{\text{pol}}}{B_{\text{tor}}}
+   \left({\boldsymbol{b}}\times{\boldsymbol{\kappa}}\right)^y =& {\sigma_y}\frac{{B_{\text{pol}}}{B_{\text{tor}}}
    R}{B^2{h_\theta}}\left[\frac{\mu_0}{B}{\frac{\partial P}{\partial x}} + {\frac{\partial B}{\partial x}}\right] \\
    \left({\boldsymbol{b}}\times{\boldsymbol{\kappa}}\right)^z =& -\frac{\mu_0}{B}{\frac{\partial P}{\partial x}} - {\frac{\partial B}{\partial x}} -
-   I\left({\boldsymbol{b}}\times{\boldsymbol{\kappa}}\right)^x\end{aligned}
-
-Curvature from `{\boldsymbol{b}}\times{\boldsymbol{\kappa}}`
+   \sigma_{B\text{pol}}I\left({\boldsymbol{b}}\times{\boldsymbol{\kappa}}\right)^x\end{aligned}
+
+.. [#curvature] Note on signs: `\nabla\times\boldsymbol{b}` should
+                flip sign if we flip the magnetic field direction (i.e.
+                `B_\text{pol}\rightarrow -B_\text{pol}` and
+                `B_\text{tor} \rightarrow -B_\text{tor}`). Under this
+                flip, the `x`-coordinate stays the same and the
+                `z`-coordinate flips sign. The `y`-coordinate stays
+                the same, or the `\eta`-coordinate flips sign.
+                Therefore the x-component of `\nabla\times\boldsymbol{b}` 
+                should flip sign, the `z`-component should not flip
+                sign (product of two sign flips), and the
+                `y`-component should flip sign if '`y`' is `y` and not
+                flip sign if '`y`' is `\eta`.
+
+
+Curvature from `{\nabla\times\left(\frac{\boldsymbol{b}}{B}\right)}`
 ------------------------
 
 The vector `{\boldsymbol{b}}\times{\boldsymbol{\kappa}}` is an
@@ -571,39 +796,41 @@ operator in curvilinear coordinates (see appendix):
    \begin{aligned}
    \nabla\times\left(\frac{{\boldsymbol{b}}}{B}\right) =&
        \frac{{B_{\text{pol}}}}{{h_\theta}}\left[\left({\frac{\partial }{\partial x}}\left(\frac{{h_\theta}}{{B_{\text{pol}}}}\right) -
-       {\frac{\partial }{\partial y}}\left(\frac{{\sigma_{B\theta}}{B_{\text{tor}}}IR}{B^2}\right)\right){\boldsymbol{e}}_z \right.  \\ &+
-       {\frac{\partial }{\partial y}}\left(\frac{{\sigma_{B\theta}}{B_{\text{tor}}}R}{B^2}\right){\boldsymbol{e}}_x \\ &+
-       \left.{\frac{\partial }{\partial x}}\left(\frac{{\sigma_{B\theta}}{B_{\text{tor}}}R}{B^2}\right){\boldsymbol{e}}_y\right]\end{aligned}
+       \sigma_y{\frac{\partial }{\partial y}}\left(\frac{{\sigma_{B\text{pol}}}{B_{\text{tor}}}IR}{B^2}\right)\right){\boldsymbol{e}}_z \right.  \\ &+
+       \sigma_y{\frac{\partial }{\partial y}}\left(\frac{{B_{\text{tor}}}R}{B^2}\right){\boldsymbol{e}}_x \\ &-
+       \sigma_y\left.{\frac{\partial }{\partial x}}\left(\frac{{B_{\text{tor}}}R}{B^2}\right){\boldsymbol{e}}_y\right]\end{aligned}
 
 This can be simplified using
 
 .. math::
 
    \begin{aligned}
-   {\frac{\partial }{\partial y}}\left(\frac{{\sigma_{B\theta}}{B_{\text{tor}}}IR}{B^2}\right) = I{\sigma_{B\theta}}{B_{\text{tor}}}
-   R{\frac{\partial }{\partial y}}\left(\frac{1}{B^2}\right) + \frac{{B_{\text{tor}}}R}{B^2}{\frac{\partial \nu}{\partial x}}\end{aligned}
+   {\sigma_y\frac{\partial }{\partial y}}\left(\frac{{\sigma_{B\text{pol}}}{B_{\text{tor}}}IR}{B^2}\right) = I{\sigma_{B\text{pol}}}{B_{\text{tor}}}
+   R{\sigma_y\frac{\partial }{\partial y}}\left(\frac{1}{B^2}\right) + \frac{{B_{\text{tor}}}R}{B^2}{\frac{\partial \nu}{\partial x}}\end{aligned}
 
 to give
 
 .. math::
 
    \begin{aligned}
-     \left({\boldsymbol{b}}\times{\boldsymbol{\kappa}}\right)^x =& -{\sigma_{B\theta}}\frac{{B_{\text{pol}}}{B_{\text{tor}}}R}{{h_\theta}B^2}{\frac{\partial B}{\partial y}} \\
-       \left({\boldsymbol{b}}\times{\boldsymbol{\kappa}}\right)^y =& -{\sigma_{B\theta}}\frac{B{B_{\text{pol}}}}{2{h_\theta}}{\frac{\partial }{\partial x}}\left(\frac{{B_{\text{tor}}}
-   R}{B^2}\right) \\ \left({\boldsymbol{b}}\times{\boldsymbol{\kappa}}\right)^z =&
-       \frac{B{B_{\text{pol}}}}{2{h_\theta}}{\frac{\partial }{\partial x}}\left(\frac{{h_\theta}}{{B_{\text{pol}}}}\right) - \frac{{B_{\text{pol}}}{B_{\text{tor}}}
-       R}{2{h_\theta}B}{\frac{\partial \nu}{\partial x}} - I\left({\boldsymbol{b}}\times{\boldsymbol{\kappa}}\cdot\nabla\right)^x\end{aligned}
+     {\frac{B}{2}\left(\nabla\times\frac{\boldsymbol{b}}{B}\right)^x} =& {-{\sigma_y}\frac{{B_{\text{pol}}}{B_{\text{tor}}}R}{{h_\theta}B^2}{\frac{\partial B}{\partial y}}} \\
+     {\frac{B}{2}\left(\nabla\times\frac{\boldsymbol{b}}{B}\right)^y} =& {-{\sigma_y}\frac{B{B_{\text{pol}}}}{2{h_\theta}}{\frac{\partial }{\partial x}}\left(\frac{{B_{\text{tor}}}
+   R}{B^2}\right)} \\
+     {\frac{B}{2}\left(\nabla\times\frac{\boldsymbol{b}}{B}\right)^z} =&
+       {\frac{B{B_{\text{pol}}}}{2{h_\theta}}{\frac{\partial }{\partial x}}\left(\frac{{h_\theta}}{{B_{\text{pol}}}}\right) - \frac{{B_{\text{pol}}}{B_{\text{tor}}}
+       R}{2{h_\theta}B}{\frac{\partial \nu}{\partial x}} - I\frac{B}{2}\left(\nabla\times\frac{\boldsymbol{b}}{B}\right)^x}
+   \end{aligned}
 
 The first and second terms in
-`\left({\boldsymbol{b}}\times{\boldsymbol{\kappa}}\cdot\nabla\right)^z`
+`\frac{B}{2}\left(\nabla\times\frac{\boldsymbol{b}}{B}\right)^z`
 almost cancel, so by expanding out `\nu` a better expression is
 
 .. math::
 
    \begin{aligned}
-   \left({\boldsymbol{b}}\times{\boldsymbol{\kappa}}\right)^z = \frac{{B_{\text{pol}}}^3}{2{h_\theta}
+   \frac{B}{2}\left(\nabla\times\frac{\boldsymbol{b}}{B}\right)^z = \frac{{B_{\text{pol}}}^3}{2{h_\theta}
    B}{\frac{\partial }{\partial x}}\left(\frac{{h_\theta}}{{B_{\text{pol}}}}\right) - \frac{{B_{\text{tor}}}
-   R}{2B}{\frac{\partial }{\partial x}}\left(\frac{{h_\theta}}{{B_{\text{pol}}}}\right)\end{aligned}
+   R}{2B}{\frac{\partial }{\partial x}}\left(\frac{{B_\text{tor}}}{R}\right)\end{aligned}
 
 Curvature of a single line
 --------------------------
@@ -673,7 +900,7 @@ An alternative is to use
 .. math::
 
    \begin{aligned}
-   {\boldsymbol{b}}\times \nabla\phi = \frac{{\sigma_{B\theta}}}{BR^2}\nabla\psi\end{aligned}
+   {\boldsymbol{b}}\times \nabla\phi = \frac{{\sigma_{B\text{pol}}}}{BR^2}\nabla\psi\end{aligned}
 
 and that the tangent vector `{\boldsymbol{T}} =
 {\boldsymbol{b}}`. This gives
@@ -682,7 +909,7 @@ and that the tangent vector `{\boldsymbol{T}} =
    :label: eq:flinenablapsi
 
    \begin{aligned}
-   \nabla\psi = {\sigma_{B\theta}}BR\left[\frac{dR}{ds}{\boldsymbol{Z}} - \frac{dZ}{ds}{\boldsymbol{R}}\right]
+   \nabla\psi = {\sigma_{B\text{pol}}}BR\left[\frac{dR}{ds}{\boldsymbol{Z}} - \frac{dZ}{ds}{\boldsymbol{R}}\right]
    \end{aligned}
 
 and so because
@@ -692,7 +919,7 @@ and so because
    :label: eq:flinekappsi
 
    \begin{aligned}
-   {\boldsymbol{\kappa}}\cdot\nabla\psi = {\sigma_{B\theta}}BR\left[ \left( \frac{{B_{\text{tor}}}^2}{RB^2} -
+   {\boldsymbol{\kappa}}\cdot\nabla\psi = {\sigma_{B\text{pol}}}BR\left[ \left( \frac{{B_{\text{tor}}}^2}{RB^2} -
    {\frac{d^2 R}{d s^2}}\right){\frac{d Z}{d s}} + {\frac{d^2 Z}{d s^2}}\frac{dR}{ds} \right]
    \end{aligned}
 
@@ -714,7 +941,7 @@ The components in field-aligned coordinates can then be calculated:
 .. math::
 
    \begin{aligned}
-   \left({\boldsymbol{b}}\times{\boldsymbol{\kappa}}\right)^x =& {\sigma_{B\theta}}\left({\boldsymbol{b}}\times{\boldsymbol{\kappa}}\right)\cdot\nabla\psi \\ =&
+   \left({\boldsymbol{b}}\times{\boldsymbol{\kappa}}\right)^x =& {\sigma_{B\text{pol}}}\left({\boldsymbol{b}}\times{\boldsymbol{\kappa}}\right)\cdot\nabla\psi \\ =&
        \frac{R{B_{\text{pol}}}^2}{B}\left(2{\frac{d R}{d s}}{\frac{d \phi}{d s}} + R{\frac{d^2 \phi}{d s^2}}\right) -
        R{B_{\text{tor}}}\left({\frac{d R}{d s}}{\frac{d^2 R}{d s^2}} + {\frac{d Z}{d s}}{\frac{d^2 Z}{d s^2}}\right) +
        \frac{{B_{\text{tor}}}^3}{B^2}{\frac{d R}{d s}}\end{aligned}
@@ -750,7 +977,7 @@ Therefore, in terms of unit vectors:
 
    \begin{aligned}
    \nabla\times{\boldsymbol{b}}=
-   \frac{1}{R{h_\theta}}{\frac{\partial }{\partial \theta}}\left(\frac{R{B_{\text{tor}}}}{B}\right){\hat{{\boldsymbol{e}}}}_\psi -
+   \frac{\sigma_{B\text{pol}}}{R{h_\theta}}{\frac{\partial }{\partial \theta}}\left(\frac{R{B_{\text{tor}}}}{B}\right){\hat{{\boldsymbol{e}}}}_\psi -
    {B_{\text{pol}}}{\frac{\partial }{\partial \psi}}\left(\frac{R{B_{\text{tor}}}}{B}\right){\hat{{\boldsymbol{e}}}}_\theta + \frac{{B_{\text{pol}}}
    R}{{h_\theta}}{\frac{\partial }{\partial \psi}}\left(\frac{{h_\theta}{B_{\text{pol}}}}{B}\right){\hat{{\boldsymbol{e}}}}_\phi\end{aligned}
 
@@ -764,7 +991,7 @@ using local derivatives (implemented using DCT).
 .. math::
 
    \begin{aligned}
-   {B_{\text{pol}}}= \frac{\left|\nabla\psi\right|}{R} = \frac{1}{R}\sqrt{\left({\frac{\partial \psi}{\partial R}}\right)^2 +
+   {|B_{\text{pol}}|}= \frac{\left|\nabla\psi\right|}{R} = \frac{1}{R}\sqrt{\left({\frac{\partial \psi}{\partial R}}\right)^2 +
    \left({\frac{\partial \psi}{\partial R}}\right)^2}\end{aligned}
 
 Using
@@ -843,7 +1070,7 @@ From the definition
 Parallel derivative of the B field
 ----------------------------------
 
-To get the parallel nablaients of the `B` field components, start
+To get the parallel gradients of the `B` field components, start
 with
 
 .. math::
@@ -954,36 +1181,36 @@ psi derivative of h
 
 From the expression for curvature (equation :eq:`eq:curvature`),
 and using
-`\nabla x \cdot \nabla \psi = {\sigma_{B\theta}}\left(R{B_{\text{pol}}}\right)^2`
+`\nabla x \cdot \nabla \psi = {\sigma_{B\text{pol}}}\left(R{B_{\text{pol}}}\right)^2`
 and
-`\nabla z\cdot\nabla \psi = -{\sigma_{B\theta}}I \left(R{B_{\text{pol}}}\right)^2`
+`\nabla z\cdot\nabla \psi = -I \left(R{B_{\text{pol}}}\right)^2`
 
 .. math::
 
    \begin{aligned}
-   {\boldsymbol{\kappa}}\cdot\nabla\psi =& -{\sigma_{B\theta}}
+   {\boldsymbol{\kappa}}\cdot\nabla\psi =& -{\sigma_{B\text{pol}}}
        \frac{{B_{\text{pol}}}}{B{h_\theta}}{\left({R{B_{\text{pol}}}}\right)^2}\left[{\frac{\partial }{\partial x}}\left(\frac{B{h_\theta}}{{B_{\text{pol}}}}\right) -
-       {\sigma_{B\theta}}{\frac{\partial }{\partial y}}\left(\frac{{B_{\text{tor}}}IR}{B}\right)\right] \\ &- I{\left({R{B_{\text{pol}}}}\right)^2}
+       {\sigma_y\sigma_{B\text{pol}}}{\frac{\partial }{\partial y}}\left(\frac{{B_{\text{tor}}}IR}{B}\right)\right] \\ &- \sigma_yI{\left({R{B_{\text{pol}}}}\right)^2}
            \frac{{B_{\text{pol}}}}{B{h_\theta}}{\frac{\partial }{\partial y}}\left(\frac{{B_{\text{tor}}}R}{B}\right)\end{aligned}
 
 The second and third terms partly cancel, and using
-`{\frac{\partial I}{\partial y}} = {\sigma_{B\theta}}
+`\sigma_y\sigma_{B\text{pol}}{\frac{\partial I}{\partial y}} = 
 {\frac{\partial \nu}{\partial x}}`
 
 .. math::
 
    \begin{aligned}
      \frac{{\boldsymbol{\kappa}}\cdot\nabla\psi}{{\left({R{B_{\text{pol}}}}\right)^2}} =&
-       -{\sigma_{B\theta}}\frac{{B_{\text{pol}}}}{B{h_\theta}}{\frac{\partial }{\partial x}}\left(\frac{B{h_\theta}}{{B_{\text{pol}}}}\right) +
-       {\sigma_{B\theta}}\frac{{B_{\text{pol}}}}{B{h_\theta}}\frac{{B_{\text{tor}}}R}{B}{\frac{\partial \nu}{\partial x}} \\ =&
-       -{\sigma_{B\theta}}\frac{{B_{\text{pol}}}}{B{h_\theta}}\left[{\frac{\partial }{\partial x}}\left(\frac{B{h_\theta}}{{B_{\text{pol}}}}\right) - \frac{{B_{\text{tor}}}
+       -{\sigma_{B\text{pol}}}\frac{{B_{\text{pol}}}}{B{h_\theta}}{\frac{\partial }{\partial x}}\left(\frac{B{h_\theta}}{{B_{\text{pol}}}}\right) +
+       {\sigma_{B\text{pol}}}\frac{{B_{\text{pol}}}}{B{h_\theta}}\frac{{B_{\text{tor}}}R}{B}{\frac{\partial \nu}{\partial x}} \\ =&
+       -{\sigma_{B\text{pol}}}\frac{{B_{\text{pol}}}}{B{h_\theta}}\left[{\frac{\partial }{\partial x}}\left(\frac{B{h_\theta}}{{B_{\text{pol}}}}\right) - \frac{{B_{\text{tor}}}
        R}{B}{\frac{\partial }{\partial x}}\left(\frac{{B_{\text{tor}}}{h_\theta}}{{B_{\text{pol}}}R}\right)\right] \\ =&
-               -{\sigma_{B\theta}}\frac{{B_{\text{pol}}}}{B{h_\theta}}\left[{h_\theta}{\frac{\partial }{\partial x}}\left(\frac{B}{{B_{\text{pol}}}}\right) -
+               -{\sigma_{B\text{pol}}}\frac{{B_{\text{pol}}}}{B{h_\theta}}\left[{h_\theta}{\frac{\partial }{\partial x}}\left(\frac{B}{{B_{\text{pol}}}}\right) -
                {h_\theta}\frac{{B_{\text{tor}}}R}{B}{\frac{\partial }{\partial x}}\left(\frac{{B_{\text{tor}}}}{{B_{\text{pol}}}R}\right) +
            \frac{B^2}{B{B_{\text{pol}}}}{\frac{\partial {h_\theta}}{\partial x}} -
-       \frac{{B_{\text{tor}}}^2}{B{B_{\text{pol}}}}{\frac{\partial {h_\theta}}{\partial x}}\right] \\ =& -{\sigma_{B\theta}}
-           \frac{{B_{\text{pol}}}}{B^2{h_\theta}}{\frac{\partial {h_\theta}}{\partial x}} -
-           {\sigma_{B\theta}}\frac{{B_{\text{pol}}}}{B^2}\left[B{\frac{\partial }{\partial x}}\left(\frac{B}{{B_{\text{pol}}}}\right) - {B_{\text{tor}}}
+       \frac{{B_{\text{tor}}}^2}{B{B_{\text{pol}}}}{\frac{\partial {h_\theta}}{\partial x}}\right] \\ =& -{\sigma_{B\text{pol}}}
+           \frac{{B_{\text{pol}}^2}}{B^2{h_\theta}}{\frac{\partial {h_\theta}}{\partial x}} -
+           {\sigma_{B\text{pol}}}\frac{{B_{\text{pol}}}}{B^2}\left[B{\frac{\partial }{\partial x}}\left(\frac{B}{{B_{\text{pol}}}}\right) - {B_{\text{tor}}}
            R{\frac{\partial }{\partial x}}\left(\frac{{B_{\text{tor}}}}{{B_{\text{pol}}}R}\right)\right]\end{aligned}
 
 Writing
@@ -1005,7 +1232,7 @@ this simplifies to give
 
    \begin{aligned}
    \frac{{\boldsymbol{\kappa}}\cdot\nabla\psi}{{\left({R{B_{\text{pol}}}}\right)^2}} =
-   -{\sigma_{B\theta}}\frac{{B_{\text{pol}}}^2}{B^2{h_\theta}}{\frac{\partial {h_\theta}}{\partial x}} - {\sigma_{B\theta}}\frac{{B_{\text{tor}}}^2}{B^2
+   -{\sigma_{B\text{pol}}}\frac{{B_{\text{pol}}}^2}{B^2{h_\theta}}{\frac{\partial {h_\theta}}{\partial x}} - {\sigma_{B\text{pol}}}\frac{{B_{\text{tor}}}^2}{B^2
    R}{\frac{\partial R}{\partial x}}
    \end{aligned}
 
@@ -1024,7 +1251,7 @@ Using :eq:`eq:flinenablapsi`,
 .. math::
 
    \begin{aligned}
-   \nabla\psi \cdot \nabla R = -{\sigma_{B\theta}}B R\frac{dZ}{ds}\end{aligned}
+   \nabla\psi \cdot \nabla R = -{\sigma_{B\text{pol}}}B R\frac{dZ}{ds}\end{aligned}
 
 and so
 
@@ -1480,7 +1707,7 @@ We have
 
    \begin{aligned}
        B{\overset{\text{def}}{=}}& \sqrt{{\boldsymbol{B}}\cdot{\boldsymbol{B}}} = \sqrt{J^{-1}{\boldsymbol{e}}_y\cdot
-   J^{-1}{\boldsymbol{e}}_y} = \sqrt{J^{-2}g_{yy}} = J^{-1}\sqrt{g_{yy}}\end{aligned}
+   J^{-1}{\boldsymbol{e}}_y} = \sqrt{J^{-2}g_{yy}}\end{aligned}
 
 Further on
 
@@ -1488,7 +1715,7 @@ Further on
 
    \begin{aligned}
        {\boldsymbol{B}}=&B{\boldsymbol{b}}\\ {\boldsymbol{b}}=&\frac{{\boldsymbol{B}}}{B}
-       =\frac{J^{-1}{\boldsymbol{e}}_y}{J^{-1}\sqrt{g_{yy}}} =\frac{{\boldsymbol{e}}_y}{\sqrt{g_{yy}}}\end{aligned}
+       =\frac{J^{-1}{\boldsymbol{e}}_y}{\sqrt{J^{-2}g_{yy}}} =\frac{\sigma_{B\text{pol}}{\boldsymbol{e}}_y}{\sqrt{g_{yy}}}\end{aligned}
 
 The parallel and perpendicular gradients
 ----------------------------------------
@@ -1535,7 +1762,7 @@ The perpendicular gradients in Laplacian inversion
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
 In the Laplacian inversion BOUT++ currently neglects the parallel
-`y` derivatives if `g_{xy}` and `g_{yz}` are non-zero,
+`y` derivatives even if `g_{xy}` and `g_{yz}` are non-zero,
 thus
 
 .. math::
@@ -1736,8 +1963,8 @@ found in incompressible fluid flow
    \begin{aligned}
        {\boldsymbol{v}}_E =& -\frac{\nabla\phi\times{\boldsymbol{b}}}{B}\\
                 %
-                =&-\frac{\nabla\phi\times{\boldsymbol{e}}_y}{
-   \sqrt{g_{yy}}J^{-1}\sqrt{g_{yy}}}\\
+                =&-\frac{\nabla\phi\times{\sigma_{B\text{pol}}\boldsymbol{e}}_y}{
+   \sqrt{g_{yy}J^{-2}}\sqrt{g_{yy}}}\\
                 %
                 =&-\frac{J}{g_{yy}}\nabla\phi\times{\boldsymbol{e}}_y\\
                 %
@@ -1758,10 +1985,10 @@ found in incompressible fluid flow
    g_{yy}{\boldsymbol{e}}^y\times{\boldsymbol{e}}^z\partial_z + g_{yz}{\boldsymbol{e}}^z\times{\boldsymbol{e}}^z\partial_z
    \right) \phi\\
                 %
-                =&\frac{J}{g_{yy}} \left( - g_{yy}{\boldsymbol{e}}^y\times{\boldsymbol{e}}^x\partial_x +
+                =&\frac{J}{g_{yy}} \left( - g_{yy}{\boldsymbol{e}}^x\times{\boldsymbol{e}}^y\partial_x +
    g_{yz}{\boldsymbol{e}}^z\times{\boldsymbol{e}}^x\partial_x \right.  \\ &\quad +
-   g_{yx}{\boldsymbol{e}}^x\times{\boldsymbol{e}}^y\partial_y - g_{yz}{\boldsymbol{e}}^z\times{\boldsymbol{e}}^y\partial_y
-   \\ &\quad \left.  - g_{yx}{\boldsymbol{e}}^x\times{\boldsymbol{e}}^z\partial_z +
+   g_{yx}{\boldsymbol{e}}^x\times{\boldsymbol{e}}^y\partial_y - g_{yz}{\boldsymbol{e}}^y\times{\boldsymbol{e}}^z\partial_y
+   \\ &\quad \left.  - g_{yx}{\boldsymbol{e}}^z\times{\boldsymbol{e}}^x\partial_z +
    g_{yy}{\boldsymbol{e}}^y\times{\boldsymbol{e}}^z\partial_z \right) \phi\\
                 %
                 =&\frac{1}{g_{yy}} \left( - g_{yy}{\boldsymbol{e}}_z\partial_x +
@@ -1820,7 +2047,7 @@ The pure solenoidal advection is thus
 The brackets operator in BOUT++
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-Notice that the (phi,f)@ operators in BOUT++ returns
+Notice that the `\mathtt{(phi,f)}` operators in BOUT++ returns
 `-\frac{\nabla\phi\times{\boldsymbol{b}}}{B}\cdot\nabla f`
 rather than
 `-\nabla\phi\times{\boldsymbol{b}}\cdot\nabla f`.
@@ -1867,7 +2094,7 @@ velocity can be written as
    \frac{1}{B}{\boldsymbol{b}}\cdot\nabla\times\nabla\phi
    \end{aligned}
 
-The second term on the right is identically zero (curl of a nablaient).
+The second term on the right is identically zero (curl of a gradient).
 The first term on the right can be expanded as
 
 .. math::
@@ -1893,7 +2120,7 @@ this becomes:
      \nabla\cdot\left(\frac{1}{B}{\boldsymbol{b}}\times\nabla\phi\right) =
      &-{\boldsymbol{b}}\times\nabla\left(\frac{1}{B}\right)\cdot\nabla\phi \\ &+
      \frac{1}{B}{\boldsymbol{b}}\times{\boldsymbol{\kappa}}\cdot\nabla\phi \\ &+
-     \left[{\boldsymbol{b}}\cdot\left(\nabla\times{\boldsymbol{b}}\right)\right]{\boldsymbol{b}}\cdot\nabla\phi\end{aligned}
+     \frac{1}{B}\left[{\boldsymbol{b}}\cdot\left(\nabla\times{\boldsymbol{b}}\right)\right]{\boldsymbol{b}}\cdot\nabla\phi\end{aligned}
 
 Alternatively, equation :eq:`eq:exb1` can be expanded as
 
@@ -1904,7 +2131,7 @@ Alternatively, equation :eq:`eq:exb1` can be expanded as
        -B{\boldsymbol{b}}\times\nabla\left(\frac{1}{B^2}\right)\cdot\nabla\phi +
        \frac{1}{B^2}\nabla\times{\boldsymbol{B}}\cdot\nabla\phi \\ =&
        -B{\boldsymbol{b}}\times\nabla\left(\frac{1}{B^2}\right)\cdot\nabla\phi +
-       \frac{1}{B^2}{\boldsymbol{J}}\cdot\nabla\phi\end{aligned}
+       \frac{1}{B^2}{\mu_0\boldsymbol{J}}\cdot\nabla\phi\end{aligned}
 
 .. math::