2 Cap-cyclide geometry

In cylindrical coordinates $(R,Z)$ the vacuum Grad Shafranov equation (Shafranov Reference Shafranov1960) is

(2.1) $$\begin{eqnarray}\unicode[STIX]{x0394}^{\ast }\unicode[STIX]{x1D713}=\frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D713}}{\unicode[STIX]{x2202}R^{2}}-\frac{1}{R}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}{\unicode[STIX]{x2202}R}+\frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D713}}{\unicode[STIX]{x2202}z^{2}}=0.\end{eqnarray}$$

This former equation is formally similar to the Laplace equation, but for a sign (minus instead of plus) in the first derivative part.

(2.2) $$\begin{eqnarray}\unicode[STIX]{x0394}\unicode[STIX]{x1D713}=\frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D713}}{\unicode[STIX]{x2202}R^{2}}+\frac{1}{R}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}{\unicode[STIX]{x2202}R}+\frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D713}}{\unicode[STIX]{x2202}z^{2}}=0.\end{eqnarray}$$

Consequently, it is obvious that any analytical solution found for the Laplace equation will be correlated to the analytical solution of the Grad Shafranov equation. Since the end of the 19th century (Böcher Reference Böcher1894; Wangerin Reference Wangerin1909) analytical solutions to the Laplace equation have been worked out using a large set of system of coordinates, that allows writing of the solution in a ‘separable’ way, i.e.

(2.3) $$\begin{eqnarray}\unicode[STIX]{x1D713}(\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FD},\unicode[STIX]{x1D6FE})=A(\unicode[STIX]{x1D6FC})B(\unicode[STIX]{x1D6FD})G(\unicode[STIX]{x1D6FE}),\end{eqnarray}$$

where ( $\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FD},\unicode[STIX]{x1D6FE}$ ) are independent orthogonal variables. For a subset of the system of coordinates a so-called ‘ $R$ -separable’ solution can be worked out. In particular this is the solution for the coordinate systems with a toroidal symmetry, such as the proper ones for the tokamak experiments

(2.4) $$\begin{eqnarray}\unicode[STIX]{x1D713}(\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FD},\unicode[STIX]{x1D6FE})=A(\unicode[STIX]{x1D6FC})B(\unicode[STIX]{x1D6FD})G(\unicode[STIX]{x1D6FE})/{\mathcal{R}}(\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FD},\unicode[STIX]{x1D6FE})\end{eqnarray}$$

${\mathcal{R}}(\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FD},\unicode[STIX]{x1D706})$ being a function determined by the coordinate system.

For the present toroidal elongated tokamak, the most appropriate coordinate system is the cap-cyclide one (Moon & Spencer Reference Moon and Spencer1971), (figure 1) with

(2.5) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}x=\displaystyle \frac{\unicode[STIX]{x1D6EC}}{a_{s}\unicode[STIX]{x1D6E4}}\,\text{d}n(\unicode[STIX]{x1D708},k_{1})sn(\unicode[STIX]{x1D707},k)\cos \unicode[STIX]{x1D711}\\ y=\displaystyle \frac{\unicode[STIX]{x1D6EC}}{a_{s}\unicode[STIX]{x1D6E4}}\,\text{d}n(\unicode[STIX]{x1D708},k_{1})sn(\unicode[STIX]{x1D707},k)\sin \unicode[STIX]{x1D711}\\ z=\displaystyle \frac{k^{0.5}\unicode[STIX]{x1D6F1}}{2a_{s}\unicode[STIX]{x1D6E4}},\end{array}\right\}\end{eqnarray}$$

where ( $x,y,z$ ) are the standard Cartesian coordinates, $\text{d}n,cn,sn$ are the Jacobi elliptical functions, $a_{s}$ is a dimensional scale parameter and

(2.6) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\unicode[STIX]{x1D6EC}=1-\text{d}n^{2}(\unicode[STIX]{x1D707})sn^{2}(\unicode[STIX]{x1D708})\\ \unicode[STIX]{x1D6E4}=sn^{2}(\unicode[STIX]{x1D707})\,\text{d}n^{2}(\unicode[STIX]{x1D708})+[(\unicode[STIX]{x1D6EC}/k^{0.5})+cn(\unicode[STIX]{x1D707})\,\text{d}n(\unicode[STIX]{x1D707})sn(\unicode[STIX]{x1D708})cn(\unicode[STIX]{x1D708})]^{2}\\ \unicode[STIX]{x1D6F1}=(\unicode[STIX]{x1D6EC}^{2}/k)-[sn^{2}(\unicode[STIX]{x1D707})\,\text{d}n^{2}(\unicode[STIX]{x1D708})+cn^{2}(\unicode[STIX]{x1D707})\,\text{d}n^{2}(\unicode[STIX]{x1D707})sn^{2}(\unicode[STIX]{x1D708})cn^{2}(\unicode[STIX]{x1D708})].\end{array}\right\}\end{eqnarray}$$

The new variables ( $\unicode[STIX]{x1D707},\unicode[STIX]{x1D708},\unicode[STIX]{x1D6F1}$ ) are defined between

(2.7a-c ) $$\begin{eqnarray}0\leqslant \unicode[STIX]{x1D707}\leqslant K;\quad 0\leqslant \unicode[STIX]{x1D708}\leqslant K^{\prime };\quad 0\leqslant \unicode[STIX]{x1D711}\leqslant 2\unicode[STIX]{x03C0}.\end{eqnarray}$$

Here $K$ and $iK^{\prime }$ are respectively the real and the imaginary complete elliptical integrals

(2.8a-c ) $$\begin{eqnarray}\displaystyle K(k)=\int _{0}^{\unicode[STIX]{x03C0}/2}\frac{d\unicode[STIX]{x1D717}}{(1-k^{2}\sin ^{2}\unicode[STIX]{x1D717})^{1/2}};\quad \text{i}K^{\prime }(k_{1})=\text{i}\int _{0}^{\unicode[STIX]{x03C0}/2}\frac{\text{d}\unicode[STIX]{x1D717}}{(1-k_{1}^{2}\sin ^{2}\unicode[STIX]{x1D717})^{1/2}};\quad k^{2}+k_{1}^{2}=1 & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

and $k$ and $k_{1}$ are respectively the parameter and the complementary parameter of the elliptical integral.

Figure 1.

Cap-cyclide coordinate system. Surfaces at constant $\unicode[STIX]{x1D707}$ and $v$ .

The coordinate transformation of (2.5) can be also written in the complex plane as

(2.9) $$\begin{eqnarray}(x^{2}+y^{2})^{1/2}+\text{i}z=\frac{1}{a_{s}sn(\unicode[STIX]{x1D707}+\text{i}\unicode[STIX]{x1D708})+\text{i}a_{s}k^{-0.5}}+\frac{1}{2a_{s}k^{-0.5}}.\end{eqnarray}$$

Writing from here on out $\text{d}n(\unicode[STIX]{x1D707})$ , $cn(\unicode[STIX]{x1D707})$ , $sn(\unicode[STIX]{x1D707})=\text{d}n\unicode[STIX]{x1D707}$ , $cn\unicode[STIX]{x1D707}$ , $sn\unicode[STIX]{x1D707}$ and $\text{d}n(\unicode[STIX]{x1D708})$ , $cn(\unicode[STIX]{x1D708})$ , $sn(\unicode[STIX]{x1D708})=\text{d}n\unicode[STIX]{x1D708}$ , $cn\unicode[STIX]{x1D708}$ , $sn\unicode[STIX]{x1D708}$ , for $\unicode[STIX]{x1D707}=\text{const.}$ the coordinate surfaces are ring-cyclides (see figure 1) with equation

(2.10) $$\begin{eqnarray}(x^{2}+y^{2}+z^{2})^{2}+A_{\unicode[STIX]{x1D707}}(x^{2}+y^{2})+B_{\unicode[STIX]{x1D707}}z^{2}+C_{\unicode[STIX]{x1D707}}=0,\end{eqnarray}$$

with

(2.11a,b ) $$\begin{eqnarray}A_{\unicode[STIX]{x1D707}}=-(R_{\unicode[STIX]{x1D707}}^{1}+R_{\unicode[STIX]{x1D707}}^{2});\quad C_{\unicode[STIX]{x1D707}}=R_{\unicode[STIX]{x1D707}}^{1}R_{\unicode[STIX]{x1D707}}^{2}\end{eqnarray}$$

and

(2.12) $$\begin{eqnarray}B_{\unicode[STIX]{x1D707}}=\frac{4a_{s}^{2}(sn^{2}\unicode[STIX]{x1D707}+1/k)^{2}}{k(1/k-sn^{2}\unicode[STIX]{x1D707})^{2}}\left[-R_{\unicode[STIX]{x1D707}}^{1}R_{\unicode[STIX]{x1D707}}^{2}+(R_{\unicode[STIX]{x1D707}}^{1}+R_{\unicode[STIX]{x1D707}}^{2})\frac{sn^{2}\unicode[STIX]{x1D707}}{a_{s}^{2}(sn^{2}\unicode[STIX]{x1D707}+1/k)^{2}}-\frac{k^{2}}{16a_{s}^{4}}\right],\end{eqnarray}$$

where (note the definitions reported in Moon & Spencer (Reference Moon and Spencer1971) present a small mistake)

(2.13) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle R_{1}(\unicode[STIX]{x1D707})=\frac{ksn^{2}\unicode[STIX]{x1D707}\left(1-{\displaystyle \frac{\text{d}n^{2}\unicode[STIX]{x1D707}}{1+k}}\right)^{2}}{a_{s}^{2}\left\{ksn^{2}\unicode[STIX]{x1D707}+\left[\left(1-{\displaystyle \frac{\text{d}n^{2}\unicode[STIX]{x1D707}}{1+k}}\right){\displaystyle \frac{1}{k^{0.5}}}+{\displaystyle \frac{cn\unicode[STIX]{x1D707}\,\text{d}n\unicode[STIX]{x1D707}}{1+k}}k^{0.5}\right]^{2}\right\}^{2}}\\ \displaystyle R_{2}(\unicode[STIX]{x1D707})=\frac{ksn^{2}\unicode[STIX]{x1D707}\left(1-{\displaystyle \frac{\text{d}n^{2}\unicode[STIX]{x1D707}}{1+k}}\right)^{2}}{a_{s}^{2}\left\{ksn^{2}\unicode[STIX]{x1D707}+\left[\left(1-{\displaystyle \frac{\text{d}n^{2}\unicode[STIX]{x1D707}}{1+k}}\right){\displaystyle \frac{1}{k^{0.5}}}-{\displaystyle \frac{cn\unicode[STIX]{x1D707}\,\text{d}n\unicode[STIX]{x1D707}}{1+k}}k^{0.5}\right]^{2}\right\}^{2}}.\end{array}\right\} & & \displaystyle\end{eqnarray}$$

For $\unicode[STIX]{x1D708}=\text{const.}$ the coordinate surfaces are cap-cyclides (see figure 1) with equation

(2.14) $$\begin{eqnarray}(x^{2}+y^{2}+z^{2})^{2}+A_{\unicode[STIX]{x1D708}}(x^{2}+y^{2})+B_{\unicode[STIX]{x1D708}}z^{2}+C_{\unicode[STIX]{x1D708}}=0,\end{eqnarray}$$

where (note the definitions reported in Moon & Spencer (Reference Moon and Spencer1971) present a small mistake)

(2.15) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle A_{\unicode[STIX]{x1D708}}=\frac{1}{8a_{s}^{2}}\frac{(k+\text{d}n^{2}\unicode[STIX]{x1D708})^{2}}{\text{d}n^{2}\unicode[STIX]{x1D708}}\left[\frac{(cn^{4}\unicode[STIX]{x1D708}+6ksn^{2}\unicode[STIX]{x1D708}cn^{2}\unicode[STIX]{x1D708}+k^{2}sn^{4}\unicode[STIX]{x1D708})(\text{d}n^{2}\unicode[STIX]{x1D708}-k)^{2}}{(cn^{2}\unicode[STIX]{x1D708}-ksn^{2}\unicode[STIX]{x1D708})^{2}(k+\text{d}n^{2}\unicode[STIX]{x1D708})^{2}}-1\right]\\ \displaystyle B_{\unicode[STIX]{x1D708}}=\frac{k}{2a_{s}^{2}}\frac{cn^{4}\unicode[STIX]{x1D708}+6ksn^{2}\unicode[STIX]{x1D708}cn^{2}\unicode[STIX]{x1D708}+k^{2}sn^{4}\unicode[STIX]{x1D708}}{(cn^{2}\unicode[STIX]{x1D708}-ksn^{2}\unicode[STIX]{x1D708})^{2}}\\ \displaystyle C_{\unicode[STIX]{x1D708}}=\frac{k^{2}}{16a_{s}^{4}}.\end{array}\right\} & & \displaystyle\end{eqnarray}$$

By varying the parameter $k$ , the coordinates transformation of equation (2.5) describes a large set of quite different geometries (figure 2). For $\unicode[STIX]{x1D707}\rightarrow 0$ , independently of $k$ , the geometry resembles the standard toroidal geometry; for larger values of $\unicode[STIX]{x1D707}$ the shape of the constant $\unicode[STIX]{x1D707}$ surfaces depend on the value of the $k$ parameter. For $k\rightarrow 0$ (figure 2 a) the surfaces tends to assume a bean shape, like the old PBX-M (England et al. Reference England, Bell, Hirshman, Kaita, Kugel and LeBlanc1997) experiment. For intermediate values of $k$ (figure 2 b) the surfaces can be either $D$ or purely elliptical prolate shaped, like the largest part of the present ongoing tokamak experiments. For $k\rightarrow 1$ (figure 2 c) all the surfaces are similar to the standard toroidal ones, independently of the $\unicode[STIX]{x1D707}$ value.

Figure 2.

(a) For $k\rightarrow 0$ the constant $\unicode[STIX]{x1D707}$ surfaces assume a bean shape; (b) for intermediate $k$ the surfaces can be either D or purely elliptical prolate; (c) for $k\rightarrow 1$ the surfaces tend to the standard toroidal ones.

For $\unicode[STIX]{x1D707}\rightarrow K$ the surface at constant $\unicode[STIX]{x1D707}$ degenerates in a circumference arc of radius

(2.16) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle R_{0}=\sqrt{X_{0}^{2}(K,\unicode[STIX]{x1D708})+Z_{0}^{2}(K,\unicode[STIX]{x1D708})}=\frac{\sqrt{k}}{2a_{s}},\quad \text{where}\\ \displaystyle X_{0}(K,\unicode[STIX]{x1D708})=\frac{k\,\text{d}n\unicode[STIX]{x1D708}}{a_{s}(k+\text{d}n^{2}\unicode[STIX]{x1D708})};\quad Z_{0}(K,\unicode[STIX]{x1D708})=\frac{\sqrt{k}(\text{d}n^{2}\unicode[STIX]{x1D708}-k)}{2a_{s}(\text{d}n^{2}\unicode[STIX]{x1D708}+k)}.\end{array}\right\} & & \displaystyle\end{eqnarray}$$

For $k\rightarrow 1$ , $R_{0}=1/2a_{s}$ , as is clear from the previous discussion, is the centre of the standard toroidal coordinates. For the constant $\unicode[STIX]{x1D707}$ contours it is possible to define the local major radius, the local minus radius and the local aspect ratio.

(2.17) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle R(\unicode[STIX]{x1D707})=\frac{\sqrt{k}(1+ksn^{2}\unicode[STIX]{x1D707})}{2a_{s}(1+k)sn\unicode[STIX]{x1D707}};\quad a(\unicode[STIX]{x1D707})=\frac{\sqrt{k}cn\unicode[STIX]{x1D707}\,\text{d}n\unicode[STIX]{x1D707}}{2a_{s}(1+k)sn\unicode[STIX]{x1D707}}\\ \displaystyle A(\unicode[STIX]{x1D707})=\frac{R(\unicode[STIX]{x1D707})}{a(\unicode[STIX]{x1D707})}=\frac{(1+ksn^{2}\unicode[STIX]{x1D707})}{cn\unicode[STIX]{x1D707}sn\unicode[STIX]{x1D707}};\quad R_{0}=\sqrt{R^{2}(\unicode[STIX]{x1D707})+a^{2}(\unicode[STIX]{x1D707})}.\end{array}\right\} & & \displaystyle\end{eqnarray}$$

By using the expansion of the Jacobi elliptical functions in terms of the hyperbolic coordinates, when $k\rightarrow 1$ , in the limit of large aspect ratio ( $\unicode[STIX]{x1D707}\rightarrow 0$ ) we have that $A(\unicode[STIX]{x1D707})\sim \cosh (\unicode[STIX]{x1D707})$ as in the standard toroidal geometry. We can define the local elongation as $b(\unicode[STIX]{x1D707})/a(\unicode[STIX]{x1D707})$ , where $b(\unicode[STIX]{x1D707})$ is the local maximum height of any contour; after some heavy algebra we have

(2.18a,b ) $$\begin{eqnarray}b(\unicode[STIX]{x1D707})=\frac{\sqrt{k}(1-ksn^{2}\unicode[STIX]{x1D707})}{2a_{s}(1+k)sn\unicode[STIX]{x1D707}};\quad \frac{b(\unicode[STIX]{x1D707})}{a(\unicode[STIX]{x1D707})}=\frac{(1-ksn^{2}\unicode[STIX]{x1D707})}{cn\unicode[STIX]{x1D707}sn\unicode[STIX]{x1D707}}.\end{eqnarray}$$

For $\unicode[STIX]{x1D707}\rightarrow 0$ , $b(\unicode[STIX]{x1D707})/a(\unicode[STIX]{x1D707})\rightarrow 1$ independently of the parameter $k$ , i.e. the contours are circular as for the standard toroidal case. For $\unicode[STIX]{x1D707}\rightarrow K$ , $b(\unicode[STIX]{x1D707})/a(\unicode[STIX]{x1D707})\rightarrow \infty$ independently of the parameter $k$ (but $k=1$ , where $b(\unicode[STIX]{x1D707})/a(\unicode[STIX]{x1D707})\rightarrow 1$ ), this is because, in the centre of the ‘geometry’, the constant surface converges to a circumference arc and not to a point.

3 Geometry metric

For the Laplace equation (2.2) to admit a quasi-separable solution of type (2.4), in a rotational geometry like the cap-cyclides, it is necessary and sufficient that the three following conditions are satisfied (Moon & Spencer Reference Moon and Spencer1971)

(3.1) $$\begin{eqnarray}\frac{1}{f_{\unicode[STIX]{x1D707}}}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}}\left(f_{\unicode[STIX]{x1D707}}\frac{\unicode[STIX]{x2202}{\mathcal{R}}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}}\right)+\frac{1}{f_{\unicode[STIX]{x1D708}}}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D708}}\left(f_{\unicode[STIX]{x1D708}}\frac{\unicode[STIX]{x2202}{\mathcal{R}}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D708}}\right)+\unicode[STIX]{x1D6FC}_{1}(\unicode[STIX]{x1D6F7}_{11}+\unicode[STIX]{x1D6F7}_{21}){\mathcal{R}}=0\end{eqnarray}$$

(3.2a,b ) $$\begin{eqnarray}\frac{g_{\unicode[STIX]{x1D707}}}{g_{\unicode[STIX]{x1D708}}}=-(\unicode[STIX]{x1D6F7}_{13}+\unicode[STIX]{x1D6F7}_{23});\quad \frac{\sqrt{g}}{g_{i}}={\mathcal{R}}^{2}f_{i}(u_{i})F(u_{1},u_{2},u_{3}).\end{eqnarray}$$

Here ${\mathcal{R}}$ is the funcion defined in (2.4) and $g_{i}$ are the coordinate transformation metric factors

(3.3a-c ) $$\begin{eqnarray}\displaystyle g_{\unicode[STIX]{x1D707}}=g_{\unicode[STIX]{x1D708}}=\frac{1}{a_{s}^{2}\unicode[STIX]{x1D6E4}^{2}}\unicode[STIX]{x1D6EC}^{2}\unicode[STIX]{x1D6FA}^{2};\quad g_{\unicode[STIX]{x1D711}}=\frac{1}{a_{s}^{2}\unicode[STIX]{x1D6E4}^{2}}\unicode[STIX]{x1D6EC}^{2}sn^{2}\unicode[STIX]{x1D707}\,\text{d}n^{2}\unicode[STIX]{x1D707};\quad g^{1/2}=\left(\frac{\unicode[STIX]{x1D6EC}}{a_{s}\unicode[STIX]{x1D6E4}}\right)^{3}\unicode[STIX]{x1D6FA}^{2}sn\unicode[STIX]{x1D707}\,\text{d}n\unicode[STIX]{x1D707}, & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D6EC}$ and $\unicode[STIX]{x1D6E4}$ are given in (2.6) and $\unicode[STIX]{x1D6FA}^{2}=(1-sn^{2}\unicode[STIX]{x1D707}\,\text{d}n^{2}\unicode[STIX]{x1D708})(\text{d}n^{2}\unicode[STIX]{x1D708}-k^{2}sn^{2}\unicode[STIX]{x1D707})$ and $\unicode[STIX]{x1D6F7}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}$ are elements of the Stäckel matrix (Morse & Feshbach Reference Morse and Feshbach1953)

(3.4) $$\begin{eqnarray}[S]=\left[\begin{array}{@{}ccc@{}}-k^{2}sn^{2}\unicode[STIX]{x1D707} & -1 & -(k^{2}sn^{2}\unicode[STIX]{x1D707}+sn^{-2}\unicode[STIX]{x1D707})\\ \text{d}n^{2}\unicode[STIX]{x1D708} & 1 & (\text{d}n^{2}\unicode[STIX]{x1D708}+k^{2}\,\text{d}n^{-2}\unicode[STIX]{x1D708})\\ 0 & 0 & 1\end{array}\right].\end{eqnarray}$$

Consequently, from the third condition of (3.2a,b )

(3.5) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}f_{\unicode[STIX]{x1D707}}=sn\unicode[STIX]{x1D707};\quad f_{\unicode[STIX]{x1D708}}=\text{d}n\unicode[STIX]{x1D708};\quad f_{\unicode[STIX]{x1D711}}=1;\quad {\mathcal{R}}=\left({\displaystyle \frac{\unicode[STIX]{x1D6EC}}{a_{s}\unicode[STIX]{x1D6E4}}}\right)^{1/2}\\ F_{\unicode[STIX]{x1D707}}=\text{d}n\unicode[STIX]{x1D708};\quad F_{\unicode[STIX]{x1D708}}=sn\unicode[STIX]{x1D707};\quad F_{\unicode[STIX]{x1D711}}=(\unicode[STIX]{x1D6FA}^{2})/(sn\unicode[STIX]{x1D707}\,\text{d}n\unicode[STIX]{x1D708}).\end{array}\right\}\end{eqnarray}$$

Eventually, by solving the equation (3.1) we find $\unicode[STIX]{x1D6FC}_{1}=1$ .

4 Laplace equation

As a consequence of what was discussed in the previous paragraph, the Laplace equation (2.2) in the cap-cyclide coordinates admits a quasi-separable solution of the type of (2.4).

(4.1) $$\begin{eqnarray}f(\unicode[STIX]{x1D707},\unicode[STIX]{x1D708},\unicode[STIX]{x1D711})=\left(\frac{\unicode[STIX]{x1D6E4}}{\unicode[STIX]{x1D6EC}}\right)^{0.5}M(\unicode[STIX]{x1D707})N(\unicode[STIX]{x1D708})\unicode[STIX]{x1D6F7}(\unicode[STIX]{x1D711}),\end{eqnarray}$$

where the equations for the radial $M(\unicode[STIX]{x1D707})$ and the angular part $N(\unicode[STIX]{x1D708})$ can be written as

(4.2) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}{\displaystyle \frac{\text{d}^{2}M}{\text{d}\unicode[STIX]{x1D707}^{2}}}+{\displaystyle \frac{cn\unicode[STIX]{x1D707}\,\text{d}n\unicode[STIX]{x1D707}}{sn\unicode[STIX]{x1D707}}}{\displaystyle \frac{\text{d}M}{\text{d}\unicode[STIX]{x1D707}}}+\left[k^{2}sn^{2}\unicode[STIX]{x1D707}-\unicode[STIX]{x1D6FC}_{2}-\unicode[STIX]{x1D6FC}_{3}\left(k^{2}sn^{2}\unicode[STIX]{x1D707}+{\displaystyle \frac{1}{sn^{2}\unicode[STIX]{x1D707}}}\right)\right]M=0\\ {\displaystyle \frac{\text{d}^{2}N}{\text{d}\unicode[STIX]{x1D708}^{2}}}-{\displaystyle \frac{k^{2}sn\unicode[STIX]{x1D708}cn\unicode[STIX]{x1D708}}{\text{d}n\unicode[STIX]{x1D708}}}{\displaystyle \frac{\text{d}N}{\text{d}\unicode[STIX]{x1D708}}}+\left[-\text{d}n^{2}\unicode[STIX]{x1D708}+\unicode[STIX]{x1D6FC}_{2}+\unicode[STIX]{x1D6FC}_{3}\left(\text{d}n^{2}\unicode[STIX]{x1D708}+{\displaystyle \frac{k^{2}}{\text{d}n^{2}\unicode[STIX]{x1D708}}}\right)\right]N=0\\ {\displaystyle \frac{\text{d}^{2}\unicode[STIX]{x1D6F7}}{\text{d}\unicode[STIX]{x1D711}^{2}}}+\unicode[STIX]{x1D6FC}_{3}\unicode[STIX]{x1D6F7}=0.\end{array}\right\}\end{eqnarray}$$

By substituting $z_{1}=sn^{2}\unicode[STIX]{x1D707}$ and $z_{2}=\text{d}n^{2}\unicode[STIX]{x1D708}$ the two equations for $M(\unicode[STIX]{x1D707})$ and $N(\unicode[STIX]{x1D708})$ can be written as

(4.3) $$\begin{eqnarray}\frac{\text{d}^{2}Z}{\text{d}z^{2}}+\frac{1}{2}\left[\frac{1}{z-a_{1}}+\frac{1}{z-a_{2}}+\frac{2}{z-a_{3}}\right]\frac{\text{d}Z}{\text{d}z}+\frac{1}{4}\left[\frac{A_{0}+A_{1}z+A_{2}z^{2}}{(z-a_{1})(z-a_{2})(z-a_{3})^{2}}\right]Z=0.\end{eqnarray}$$

Here, the variable $z$ represents both the variables $z_{\text{1}}$ and $z_{2}$ and, consequently, the equations for the functions $M(\unicode[STIX]{x1D707})$ and $N(\unicode[STIX]{x1D708})$ are formally similar. In this equation $a_{1}=(1;1)$ , $a_{2}=(1/k^{2};k^{2})$ , $a_{3}=(0;0)$ , $A_{0}=(-q^{2}/k^{2};-q^{2}k^{2})$ , $A_{1}=(-p^{\prime 2};p^{\prime 2})$ , $A_{2}=(1-q^{2};1-q^{2})$ , where the first term in the brackets refers to the equation for $M(\unicode[STIX]{x1D707})$ , the second to the equation for $N(\unicode[STIX]{x1D708})$ , $p$ and $q$ are integer numbers and $p^{\prime }=p/k$ . Equation (4.3) is a particular case of the general Böcher (Reference Böcher1894), Moon & Spencer (Reference Moon and Spencer1971) equation

(4.4) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \frac{\text{d}^{2}\tilde{Z}}{\text{d}z^{2}}+P(z)\frac{\text{d}\tilde{Z}}{\text{d}z}+Q(z)\tilde{Z}=0\quad \text{with}\\ \displaystyle P(z)=\frac{1}{2}\left[\frac{m_{1}}{z-a_{1}}+\frac{m_{2}}{z-a_{2}}+\cdots +\frac{m_{n-1}}{z-a_{n-1}}\right];\\ \displaystyle Q(z)=\frac{1}{4}\left[\frac{A_{0}+A_{1}z+\cdots A_{l}z^{l}}{(z-a_{1})^{m_{1}}(z-a_{2})^{m_{2}}\cdots (z-a_{n-1})^{m_{n-1}}}\right].\end{array}\right\}\end{eqnarray}$$

In the particular case of (4.3) the general Böcher equation is named a Wangerin equation. This equation has two singularities of the first order for $z=1,a_{2}$ and two singularities of the second order for $z=0$ and $z\rightarrow \infty$ . Consequently the Wangerin equation is characterized by poles $\{1,1,2,2\}$ and the Laplace equation has the solution (Moon & Spencer Reference Moon and Spencer1971).

(4.5) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}M(sn^{2}\unicode[STIX]{x1D707})=AD_{p}^{q}(k,sn^{2}\unicode[STIX]{x1D707})+BF_{p}^{q}(k,sn^{2}\unicode[STIX]{x1D707})\\ N(\text{d}n^{2}\unicode[STIX]{x1D708})=AD_{p^{\prime }}^{q}(k,\text{d}n^{2}\unicode[STIX]{x1D708})+BF_{p^{\prime }}^{q}(k,\text{d}n^{2}\unicode[STIX]{x1D708})\\ \unicode[STIX]{x1D6F7}(\unicode[STIX]{x1D711})=A\sin (q\unicode[STIX]{x1D711})+B\cos (q\unicode[STIX]{x1D711}).\end{array}\right\}\end{eqnarray}$$

5 Grad Shafranov equation

The Grad Shafranov operator of equation (2.1) in axisymmetry ( $\unicode[STIX]{x2202}/\unicode[STIX]{x2202}\unicode[STIX]{x1D711}=0$ ) is similar to the toroidal part of the Laplace operator applied to the potential vector

(5.1) $$\begin{eqnarray}(\unicode[STIX]{x0394}A)_{\unicode[STIX]{x1D711}}=-\frac{1}{R}\unicode[STIX]{x0394}^{\ast }\unicode[STIX]{x1D713}=0.\end{eqnarray}$$

This expression for an orthogonal system of coordinates $(u_{1},u_{2},u_{3})$ assumes the explicit form

(5.2) $$\begin{eqnarray}\displaystyle & & \displaystyle \frac{\unicode[STIX]{x2202}^{2}A_{\unicode[STIX]{x1D711}}}{\unicode[STIX]{x2202}u_{1}^{2}}+\frac{\unicode[STIX]{x2202}^{2}A_{\unicode[STIX]{x1D711}}}{\unicode[STIX]{x2202}u_{2}^{2}}+\frac{1}{2g_{3}}\left(\frac{\unicode[STIX]{x2202}g_{3}}{\unicode[STIX]{x2202}u_{1}}\frac{\unicode[STIX]{x2202}A_{\unicode[STIX]{x1D711}}}{\unicode[STIX]{x2202}u_{1}}+\frac{\unicode[STIX]{x2202}g_{3}}{\unicode[STIX]{x2202}u_{2}}\frac{\unicode[STIX]{x2202}A_{\unicode[STIX]{x1D711}}}{\unicode[STIX]{x2202}u_{2}}\right)\nonumber\\ \displaystyle & & \displaystyle \qquad +\,\frac{A_{\unicode[STIX]{x1D711}}}{2g_{3}}\left\{\left(\frac{\unicode[STIX]{x2202}^{2}g_{3}}{\unicode[STIX]{x2202}u_{1}^{2}}+\frac{\unicode[STIX]{x2202}^{2}g_{3}}{\unicode[STIX]{x2202}u_{2}^{2}}\right)-\frac{1}{g_{3}}\left[\left(\frac{\unicode[STIX]{x2202}g_{3}}{\unicode[STIX]{x2202}u_{1}}\right)^{2}+\left(\frac{\unicode[STIX]{x2202}g_{3}}{\unicode[STIX]{x2202}u_{2}}\right)^{2}\right]\right\}=0.\end{eqnarray}$$

That in our case of the cap-cyclide coordinates becomes

(5.3) $$\begin{eqnarray}\displaystyle (\unicode[STIX]{x0394}A)_{\unicode[STIX]{x1D711}} & = & \displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}}\left[\frac{\unicode[STIX]{x2202}A_{\unicode[STIX]{x1D719}}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}}+\left(\frac{1}{f_{\unicode[STIX]{x1D707}}}\frac{\unicode[STIX]{x2202}f_{\unicode[STIX]{x1D707}}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}}-{\mathcal{R}}^{2}\frac{\unicode[STIX]{x2202}{\mathcal{R}}^{-2}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}}\right)A_{\unicode[STIX]{x1D719}}\right]\nonumber\\ \displaystyle & & \displaystyle +\,\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D708}}\left[\frac{\unicode[STIX]{x2202}A_{\unicode[STIX]{x1D719}}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D708}}+\left(\frac{1}{f_{\unicode[STIX]{x1D708}}}\frac{\unicode[STIX]{x2202}f_{\unicode[STIX]{x1D708}}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D708}}-{\mathcal{R}}^{2}\frac{\unicode[STIX]{x2202}{\mathcal{R}}^{-2}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D708}}\right)A_{\unicode[STIX]{x1D719}}\right]=0.\end{eqnarray}$$

where here ${\mathcal{R}}$ is the separation function defined in (3.5). By remembering that in our coordinates the Laplace equation in axisymmetry can be written as

(5.4) $$\begin{eqnarray}\unicode[STIX]{x0394}f=\frac{1}{{\mathcal{R}}^{2}\unicode[STIX]{x1D6FA}^{2}f_{\unicode[STIX]{x1D707}}f_{\unicode[STIX]{x1D708}}}\left[f_{\unicode[STIX]{x1D708}}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}}\left({\mathcal{R}}^{2}f_{\unicode[STIX]{x1D707}}\frac{\unicode[STIX]{x2202}f}{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}}\right)+f_{\unicode[STIX]{x1D707}}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D708}}\left({\mathcal{R}}^{2}f_{\unicode[STIX]{x1D708}}\frac{\unicode[STIX]{x2202}f}{\unicode[STIX]{x2202}\unicode[STIX]{x1D708}}\right)\right]=0.\end{eqnarray}$$

Equation (5.3) is identical to equation (5.4), plus a non-differential term in $A_{\unicode[STIX]{x1D711}}$

(5.5) $$\begin{eqnarray}\left[\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}}\left(\frac{1}{f_{\unicode[STIX]{x1D707}}}\frac{\unicode[STIX]{x2202}f_{\unicode[STIX]{x1D707}}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}}-{\mathcal{R}}^{2}\frac{\unicode[STIX]{x2202}{\mathcal{R}}^{-2}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D707}}\right)+\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D708}}\left(\frac{1}{f_{\unicode[STIX]{x1D708}}}\frac{\unicode[STIX]{x2202}f_{\unicode[STIX]{x1D708}}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D708}}-{\mathcal{R}}^{2}f_{\unicode[STIX]{x1D708}}\frac{{\mathcal{R}}^{-2}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D708}}\right)\right]=\unicode[STIX]{x1D6F7}_{13}+\unicode[STIX]{x1D6F7}_{23}.\end{eqnarray}$$

Where $\unicode[STIX]{x1D6F7}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}$ are elements of the Stäckel matrix. Consequently the set of separable solutions of (5.4) for the toroidal component of the potential vector is the same as the one of set (4.2), when assuming that the axisymmetry is $(\unicode[STIX]{x2202}/\unicode[STIX]{x2202}\unicode[STIX]{x1D711}=0\rightarrow q^{2}=\unicode[STIX]{x1D6FC}_{3}=0)$ , plus the non-differential term of (5.5).

(5.6) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \frac{\text{d}^{2}M}{\text{d}\unicode[STIX]{x1D707}^{2}}+\frac{cn\unicode[STIX]{x1D707}\,\text{d}n\unicode[STIX]{x1D707}}{sn\unicode[STIX]{x1D707}}\frac{\text{d}M}{\text{d}\unicode[STIX]{x1D707}}+\left(\frac{1}{sn^{2}\unicode[STIX]{x1D707}}-\unicode[STIX]{x1D6FC}_{2}\right)M=0\\ \displaystyle \frac{\text{d}^{2}N}{\text{d}\unicode[STIX]{x1D708}^{2}}-\frac{k^{2}sn\unicode[STIX]{x1D708}cn\unicode[STIX]{x1D708}}{\text{d}n\unicode[STIX]{x1D708}}\frac{\text{d}N}{\text{d}\unicode[STIX]{x1D708}}+\left(\unicode[STIX]{x1D6FC}_{2}+\frac{k^{2}}{\text{d}n^{2}\unicode[STIX]{x1D708}}\right)N=0.\end{array}\right\}.\end{eqnarray}$$

But this set of equations, when assuming $q^{2}=\unicode[STIX]{x1D6FC}_{3}=1$ , is identical to the one of (4.2); consequently the two-dimensional (2-D) solution of the toroidal component of the vectorial Laplace equation for the potential vector is the same as the 3-D scalar Laplace equation when assuming $q^{2}=\unicode[STIX]{x1D6FC}_{3}=1$ . If now we remember that

(5.7) $$\begin{eqnarray}\unicode[STIX]{x1D713}=\oint A_{\unicode[STIX]{x1D711}}\,\text{d}\ell =2\unicode[STIX]{x03C0}A_{\unicode[STIX]{x1D711}}\frac{\unicode[STIX]{x1D6EC}}{a_{s}\unicode[STIX]{x1D6E4}}sn\unicode[STIX]{x1D707}\,\text{d}n\unicode[STIX]{x1D707}.\end{eqnarray}$$

We get that the solution for the Grad Shafranov equation (2.1) can be expressed in terms of the Wangerin functions with $q^{2}=1$

(5.8) $$\begin{eqnarray}\unicode[STIX]{x1D713}=[AD_{p}^{1}(k,sn^{2}\unicode[STIX]{x1D707})+BF_{p}^{1}(k,sn^{2}\unicode[STIX]{x1D707})][AD_{p}^{1}(k,\text{d}n^{2}\unicode[STIX]{x1D708})+BF_{p}^{1}(k,\text{d}n^{2}\unicode[STIX]{x1D708})].\end{eqnarray}$$

It is important to note that in (5.8) all terms, but $A_{\text{p,p'}}$ and $B_{\text{p,p'}}$ , are geometrical terms. The physics information, regarding the topological structure of the magnetic field, depends only on the unknown moments $A_{\text{p,p'}}$ and $B_{\text{p.p'}}$ , constant in the vacuum region. The moments $A_{\text{p,p'}}$ describe the field produced by all the currents flowing internally to the vacuum region, whilst the moments $B_{\text{p,p'}}$ describe the field produced by all the currents flowing externally to the vacuum region. Consequently, the contribution of all the currents flowing externally from the plasma region is intrinsically separated by the internal properties of the confined plasma, which are fully described by the moments $A_{\text{p,p'}}$ . In a region with a current distribution, the solution of the Grad Shafranov equation remains formally the same as (2.1), but now the moments $A_{\text{p,p'}}$ and $B_{\text{p.p'}}$ are not constant anymore, and can be worked out by using a Green function technique, once the current distribution is known.

(5.9) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle A_{p}(r)=\int _{0}^{V}G(r^{\prime })(\unicode[STIX]{x1D707}_{0}r^{\prime }J_{\unicode[STIX]{x1D711}})\,\text{d}^{3}r^{\prime }\\ \displaystyle B_{p}(r)=\int _{V}^{\infty }G(r^{\prime })(\unicode[STIX]{x1D707}_{0}r^{\prime }J_{\unicode[STIX]{x1D711}})\,\text{d}^{3}r^{\prime }.\end{array}\right\}\end{eqnarray}$$

On a given tokamak experiment we have a set of magnetic probes, measuring the local value of the magnetic field $B$ and of the flux function $\unicode[STIX]{x1D713}$

(5.10) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \unicode[STIX]{x1D713}_{l}=\mathop{\sum }_{p.p^{\prime }}[A_{p,p^{\prime }}D_{p}^{1}(k,sn^{2}\unicode[STIX]{x1D707}_{l})+B_{p,p^{\prime }}F_{p}^{1}(k,sn^{2}\unicode[STIX]{x1D707}_{l})]\\ \quad \quad \times ~[A_{p,p^{\prime }}D_{p^{\prime }}^{1}(k,\text{d}n^{2}\unicode[STIX]{x1D708}_{l})+B_{p,p^{\prime }}F_{p^{\prime }}^{1}(k,\text{d}n^{2}\unicode[STIX]{x1D708}_{l})]\\ \displaystyle B_{l}=\frac{1}{2\unicode[STIX]{x03C0}r_{l}}\left(\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}{\unicode[STIX]{x2202}r}\right)_{l}\end{array}\right\}\end{eqnarray}$$

here the subscript $l$ indicates the generic position of a magnetic probe. Clearly this set of magnetic probes allows ua to exactly determine the values of the moments $A_{\text{p,p'}}$ and $B_{\text{p.p' }}$ and, consequently, to have the exact analytical expression of the magnetic topology in the region between the plasma column and the closer poloidal coil. Moreover, this methodology allows us to discriminate between the contribution provided by the external equilibrium currents (given by the moments $B_{\text{p,p'}}$ ) and the one provided by the internal plasma currents distribution (given by the moments $A_{\text{p,p'}}$ ). As described in Alladio & Crisanti (Reference Alladio and Crisanti1986) this approach, alongside the Green technique of (5.9) and the Grad Shafranov equation, allows us to build up a robust semi -analytical reconstructive equilibrium code. The robustness of using this approach to deal with the mathematically ill-posed problem of reconstructing the internal plasma equilibrium only relying on external magnetic measurements has been verified both in the context of the FTU circular tokamak and in the ‘old’ JET, where the divertor coils were missing. The robustness of the solution obtained, against the magnetic measurements error, was widely verified at JET and it allowed us to discover the possibility for JET to realize an X point configuration (Alladio et al. Reference Alladio, Crisanti, Lazzaro and Tanga1984; Alladio & Crisanti Reference Alladio and Crisanti1986). The drawback of the equilibrium code described in Alladio & Crisanti (Reference Alladio and Crisanti1986) was that it was worked out by using the pure toroidal geometry, with a circular poloidal section. When applied to the modern elongated tokamak, this reconstructive equilibrium code was not able to exactly discriminate between the contribution provided by the external coils and the one provided by the plasma. The solution written down in (5.10) will allow us to solve this problem and to build up a very robust semi-analytical equilibrium code, even capable of dealing with the next generation of tokamaks, where the measurements will be far from the plasma border and the signal produced by high internal magnetic moments will be from the very low up to becoming comparable with the intrinsic system noise.