Hostname: page-component-745bb68f8f-hvd4g Total loading time: 0 Render date: 2025-02-11T14:32:29.610Z Has data issue: false hasContentIssue false

Analytical solution of the Grad Shafranov equation in an elliptical prolate geometry

Published online by Cambridge University Press:  05 April 2019

F. Crisanti*
Affiliation:
ENEA for EUROfusion, via E. Fermi 45, 00044 Frascati (Rome), Italy
*
Email address for correspondence: flavio.crisanti@enea.it
Rights & Permissions [Opens in a new window]

Abstract

The analytical solution in toroidal coordinates of the Grad Shafranov equation has been at the origin of the tokamak breakthrough in the fusion development. Unfortunately, the standard toroidal coordinates have a circular poloidal section, which does not fit the elongated cross-section of the present tokamak experiments. In axisymmetry, the vacuum Grad Shafranov equation coincides with the Laplace equation for the toroidal component of the vector potential. In the present paper the solutions for the Laplace equation and that for the vacuum Grad Shafranov equation are tackled in the elliptical prolate toroidal cap-cyclide coordinates framework. The following report of the geometrical properties and of the metric of these coordinates allows us to work out the analytical solution of both equations in terms of the Wangerin functions.

Type
Research Article
Copyright
© Cambridge University Press 2019 

1 Introduction

In plasma physics the magnetic hydro dynamic (MHD) equilibrium equation $\unicode[STIX]{x1D735}P=J\times B$ is one of the most assessed equations. The analytical solution to the toroidal axisymmetric Grad Shafranov equation (Shafranov Reference Shafranov1960; Mukhovatov & Shafranov Reference Mukhovatov and Shafranov1971) was at the basis of the breakthrough realized by the tokamak configuration along the path of fusion energy achievements. In axisymmetry the Grad Shafranov equation

(1.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}}=\unicode[STIX]{x1D707}_{0}RJ_{\unicode[STIX]{x1D711}}\end{eqnarray}$$

where $J_{\unicode[STIX]{x1D711}}$ is the toroidal plasma current, $\unicode[STIX]{x1D713}$ is the flux function and ( $R,z,\unicode[STIX]{x1D711}$ ) are the standard cylindrical coordinates; in the vacuum $(J_{\unicode[STIX]{x1D711}}=0)$ coincides with the Laplace equation for the toroidal component of the vector potential (Bateman Reference Bateman1978). The analytical solution to the related Laplace equation in several different geometries was deeply studied during the second half of the 19th century (Neumann Reference Neumann1864; Böcher Reference Böcher1894; Wangerin Reference Wangerin1909), pointing to the possibility of finding an analytical solution also for oblate and prolate toroidal geometries. A good summary of all these works can be found in Moon and Spencer’s book (Reference Moon and Spencer1971). During the first half of the 20th century the problem of finding an analytical solution to the inductance of a torus was first tackled (Fock Reference Fock1932), and a first analytical solution to the Laplace equation, in elliptical oblate coordinates, was found by Lebedev (Reference Lebedev1937). More recently Chance (Reference Chance1997) has solved the Laplace equation for a general azimuthal geometry by using a Green function technique, even in presence of internal currents. Eventually, a team of Japanese scientists found the analytical solution to the vacuum Grad Shafranov equation for an oblate toroidal system (AikawaI & Takahara Reference AikawaI and Takahara1975; Honma et al. Reference Honma, Kito, Kaji and Seki1979). In the meanwhile, a large amount of scientific papers have been written to solve the tokamak equilibrium problem, adopting several different approaches. Analytical solutions to the nonlinear Grad Shafranov equation have been studied and found for simple cases of the source term (Soloviev Reference Soloviev1968; Greene Reference Greene1988; Zheng, Wootton & Solano Reference Zheng, Wootton and Solano1996; Cerfon & Freidberg Reference Cerfon and Freidberg2010). The employment of a system of toroidal coordinates, which has the flux surfaces as the radial coordinate, has allowed us to develop robust codes with the mapping solution expanded on the momentum (Boozer Reference Boozer1980; Lao, Hirschmann & Wieland Reference Lao, Hirschmann and Wieland1981; Brambilla Reference Brambilla2003). Predictive equilibrium codes have been developed, through the numerical resolution (Feneberg & Lackner Reference Feneberg and Lackner1973; Cenacchi, Galvao & Taroni Reference Cenacchi, Galvao and Taroni1976; Blum, Le Foll & Thooris Reference Blum, Le Foll and Thooris1981; Albanese, Blum & De Barbieri Reference Albanese, Blum and De Barbieri1987; Albanese, Ambrosino & Mattei Reference Albanese, Ambrosino and Mattei2015) and/or the adoption of a semi-analytical approach to deal with (Alladio & Crisanti Reference Alladio and Crisanti1986; Alladio et al. Reference Alladio, Crisanti, Marinucci, Micozzi and Tanga1991) the Grad Shafranov equation, and have been then employed to design the planned plasma configuration of practically all ongoing tokamak experiments. The linearization of the problem (Albanese & Villone Reference Albanese and Villone1998) has allowed us to easily optimize a given equilibrium configuration in terms of different aspects (currents in the poloidal coils, forces on the coils, interaction with the first wall $\ldots$ ). Moreover, the use of the linearized solver has allowed us to optimize the real time control of ongoing experiments (Sartori et al. Reference Sartori, Ambrosino, Ariola, Cenedese, Crisanti and De Tommasi2005, Reference Sartori, Crisanti, Albanese, Ambrosino, Toigo and Hay2008; Albanese et al. Reference Albanese, Ambrosino, Ariola, Artaserse, Bellizio and Coccorese2011) and to design an optimal control system for future devices (Albanese et al. Reference Albanese, Ambrosino, Ariola, Cavinato, De Tommasi and Liu2009; Bachmann et al. Reference Bachmann, Aiello, Albanese, Ambrosino, Arbeiter and Aubert2015). Although it is an ill-posed problem (different internal current distributions can produce the same magnetic signal), the magnetic measurements have been used to constrain the equilibrium codes, in order to reconstruct the experimental plasma boundary and to work out, as much as possible, the most important internal features of the confined plasma; i.e. the internal magnetic and kinetic inductances $\unicode[STIX]{x1D6FD}_{\text{p}}$ and $l_{\text{i}}$ (Shafranov Reference Shafranov1971), various plasma moments (Van Milligen Reference Van Milligen1990; Van Milligen & Lopez Fraguas Reference Van Milligen and Lopez Fraguas1994) (for instance the plasma column centroid and its elongation), the internal mapping of the magnetic surfaces and the identification of the plasma current density. As for the predictive equilibrium codes, both numerical and semi-analytical approaches have been used (Lao et al. Reference Lao, John, Stambaugh, Kellman and Pfeiffer1985; Alladio & Crisanti Reference Alladio and Crisanti1986; Alladio et al. Reference Alladio, Crisanti, Marinucci, Micozzi and Tanga1991). Eventually, the fact that the plasma equilibrium reconstruction is a mathematically ill-posed problem and that in some used codes it depends on the knowledge of external currents, and the fact that experimental noise must be accounted for, has always led to a poor reconstruction of the internal mapping of the magnetic topology. In several codes internal independent measurements have been used to constrain the solution. In particular, within the most used code, EFIT, the polarimetry (Li et al. Reference Li, Lotte, Zwingmann, Gil and Imbeaux2011) and/or the motional Stark effect (Lao et al. Reference Lao, St. John, Peng, Ferron, Strait and Taylor2005) have both been used to further constrain the problem and to improve the reconstructed poloidal flux mapping. This code has also been employed for real time reconstruction, which is to be used on the plasma shape control system (Ferron et al. Reference Ferron, Walker, Lao, St John, Humphreys and Leuer1998). This latter aspect has been strongly improved thanks to the use of powerful modern computer parallelization (Yue et al. Reference Yue, Xiao, Luo and Guo2013). Of course, the use of all these techniques and of modern computers has always led to result improvements. However, the details of the internal parameters, like the current density profile, are still far from what would be necessary for a deeper understanding of the physics and/or for an accurate real time control of the internal profiles. This represents a serious concern when considering that in a DEMO class machine the available diagnostics will be strongly reduced, with the consequent necessity of having very robust and flexible plasma models. Reconstructive equilibrium codes (Alladio & Crisanti Reference Alladio and Crisanti1986; Alladio et al. Reference Alladio, Crisanti, Marinucci, Micozzi and Tanga1991), based on a semi-analytical approach, are in principle more flexible and robust than the completely numerical ones, as they use the multipolar analytical solution of the electromagnetic problem on the most suitable function basis. In particular this approach gives the possibility to fully separate the external world (passive structures, $\text{coils},\ldots$ ) from the internal plasma features. The most important advantage given by this approach is that any ‘unphysical’ solution, which could be worked out by a pure numerical code in the presence of measurement errors, is not intrinsically possible. This aspect was widely discussed in relation to the FTU experiment’s circular toroidal geometry in Alladio & Crisanti (Reference Alladio and Crisanti1986). In that paper the robustness of the approach was also verified on the old JET experiment, without the presence of internal divertor coils. The used reconstructive equilibrium code was the one that showed the possibility in JET of producing X point configurations and the consequent possibility of achieving an improved energy confinement regime (Alladio et al. Reference Alladio, Crisanti, Lazzaro and Tanga1984). Unfortunately, so far, the exact analytical solution of the Grad Shafranov equation for the vacuum region has been found only in relation to the standard circular shaped and for the elliptical oblate toroidal geometries (AikawaI & Takahara Reference AikawaI and Takahara1975; Honma et al. Reference Honma, Kito, Kaji and Seki1979). Both these geometries are not suitable to the present tokamak experiments, all of which are based on an elliptical prolate geometry. An interesting approach to analytically solving the Grad Shafranov equation through the separable variable technique and using not constant source functions has been used in Atanasiu et al. (Reference Atanasiu, Günter, Lackner and Miron2004) and Guazzotto & Freidberg (Reference Guazzotto and Freidberg2007). In both papers the solution has been worked out by using an $(x,y)$ coordinate system. Unluckily this ‘not natural’ system of coordinates results in a solution with a large expansion that does not always converge. Guazzotto & Freidberg partially overcome the problem by linking some geometrical property of the plasma boundary to the solution expansion. Nonetheless, the non-robustness of the solution remains. In the present paper the mathematical problem is summarized and the analytical solution of the vacuum Grad Shafranov equation is worked out for the cap-cyclide coordinate system, which is the natural one for the modern elongated and D-shaped tokamak. Moreover, the paper addresses how the solution could be used to realize a robust reconstructive equilibrium code, capable of better fitting the needs of a future reactor relevant tokamak, like DEMO.

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.

6 Conclusions

Equilibrium reconstructive codes (Lao et al. Reference Lao, John, Stambaugh, Kellman and Pfeiffer1985; Alladio & Crisanti Reference Alladio and Crisanti1986; Alladio et al. Reference Alladio, Crisanti, Marinucci, Micozzi and Tanga1991) are probably the most important diagnostic tool in fusion-based experiments, and their importance will be even stronger in future DEMO class experiments, where the scarce availability of diagnostics will require very robust codes, based on fundamental physics principles. The equilibrium codes based on an analytical solution to the Grad Shafranov equation are very robust and, probably, the most suitable tools to accomplish the necessary quality in terms of robustness and reliability. Unfortunately, so far, an analytical solution to the Grad Shafranov equation was available only for circular and elliptical oblate shaped toroidal geometries, not adaptable to the present elongated tokamak. The analytical solution to the vacuum Grad Shafranov, in the cap-cyclide coordinates, as worked out in this paper, opens the possibility of realizing an equilibrium reconstructive code based on an elongated geometry. The Wangerin functions have been evaluated for the first time in AikawaI & Takahara (Reference AikawaI and Takahara1975). However, in order to use them for an equilibrium code, an independent evaluation and cross-checking will be necessary. Moreover, it will necessary to evaluate their derivatives, to be able to get the poloidal magnetic field and to open the possibility of finding the Green function for this geometry. Eventually, this will allow us to write a robust reconstructive equilibrium code, to be used on the present and future elongated tokamak experiments.

Acknowledgements

The author wants to thank all the ENEA and European colleagues that have helped with fruitful discussions. A particular thank you is due to Dr E. Giovannozzi, for the strong support received in using Matlab and Mathematica tools. This work has been carried out within the framework of the EUROfusion Consortium and has received funding from the Euratom research and training programme 2014-2018 under grant agreement no. 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission.

References

AikawaI, K. & Takahara, M.1975.Wangerin Functions. Yamanashi Univ., Rep. Facul. of Engineering 26.Google Scholar
Albanese, R., Ambrosino, G., Ariola, M., Artaserse, G., Bellizio, T., Coccorese, V. et al. 2011 Overview of modelling activities for plasma control upgrade in JET. Fus. Engin. Des. 86, 1030.Google Scholar
Albanese, R., Ambrosino, G., Ariola, M., Cavinato, M., De Tommasi, G., Liu, Y. et al. 2009 ITER vertical stabilization system. Fus. Engin. Des. 84, 394.Google Scholar
Albanese, R., Ambrosino, R. & Mattei, M. 2015 CREATE-NL $+$ : a robust control-oriented free boundary dynamic plasma equilibrium solver. Fus. Engin. Des. 96–97, 664.Google Scholar
Albanese, R., Blum, J. & De Barbieri, O.1987. Num. Stud. of the Next Europ. Torus via PROTEUS code. 12th Confer. Numeric. Simulat. of Plasmas, San Francisco.Google Scholar
Albanese, R. & Villone, F. 1998 The linearized CREATE-L plasma response model for the control of current, position and shape in tokamaks. Nucl. Fusion 38, 723.Google Scholar
Alladio, F. & Crisanti, F. 1986 Analysis of mhd equilibria by toroidal multipolar expansions. Nucl. Fusion 26, 1143.Google Scholar
Alladio, F., Crisanti, F., Lazzaro, E. & Tanga, A.1984. Observation high $\unicode[STIX]{x1D6FD}p$ effect in JET discharge. In 26th Bull. Am. Phys. Soc., Boston, (p. 3 F11).Google Scholar
Alladio, F., Crisanti, F., Marinucci, M., Micozzi, P. & Tanga, A. 1991 Analysis of tokamak configurations using the toroidal multipole method. Nucl. Fusion 31, 739.Google Scholar
Atanasiu, C. V., Günter, S., Lackner, K. & Miron, I. G. 2004 Analytical solutions to the Grad–Shafranov equation. Phys. Plasmas 11, 3510.Google Scholar
Böcher, M.1894.Üeber Die Reihenentwickelungen Der Potential Theorie. Leipzig. Druck und Verlag von B. G. Teubner.Google Scholar
Bachmann, C., Aiello, G., Albanese, R., Ambrosino, R., Arbeiter, F., Aubert, J. et al. 2015 Initial DEMO tokamak design configuration studies. Fus. Engin. Des. 98–99, 1423.Google Scholar
Bateman, G. 1978 MHD Instabilities. MIT Press.Google Scholar
Blum, J., Le Foll, J. & Thooris, B. 1981 The self-consistent equilibrium and diffusion code SCED. Comput. Phys. Commun. 24, 235.Google Scholar
Boozer, A. 1980 Guiding center drift equations. Phys. Fluids 23, 904.Google Scholar
Brambilla, M. 2003 Iterative solution of the Grad-Shafranov equation in symmetric magnetic coordinates. Phys. Plasmas 10, 3674.Google Scholar
Cenacchi, G., Galvao, R. & Taroni, A. 1976 Numerical computation of axisymmetric mhd-equilibria without conducting shell. Nucl. Fusion 16, 457.Google Scholar
Cerfon, A. J. & Freidberg, J. P. 2010 One size fits all analytic solutions to the Grad–Shafranov equation. Phys. Plasmas 3, 0352502-1.Google Scholar
Chance, M. S. 1997 Vacuum calculations in azimuthally symmetric geometry. Phys. Plasmas 4 (6), 2161.Google Scholar
England, A. C., Bell, R. E., Hirshman, S. P., Kaita, R., Kugel, H. W., LeBlanc, B. L. et al. 1997 Plasma physics and controlled fusion characterization of plasma parameters in shaped PBX-M discharges. Plasma Phys. Control. Fusion 39, 1373.Google Scholar
Feneberg, W. & Lackner, K. 1973 Multipole tokamak equilibria. Nucl. Fusion 13, 549.Google Scholar
Ferron, J. R., Walker, M. L., Lao, L. L., St John, H. E., Humphreys, D. A. & Leuer, J. A. 1998 Real time equilibrium reconstruction for tokamak discharge control. Nucl. Fusion 38 (7), 1055.Google Scholar
Fock, V. 1932 Skin effect in a ring. Fiz. Zh. Sovietunion 1, 215.Google Scholar
Greene, M. J. 1988 An analytic large aspect ratio, high beta equilibrium. Plasma Phys. Control. Fusion 30, 327.Google Scholar
Guazzotto, L. & Freidberg, J. P. 2007 A family of analytic equilibrium solutions for the Grad–Shafranov equation. Phys. Plasmas 14, 112508.Google Scholar
Honma, T., Kito, M., Kaji, I. & Seki, M.1979.The Vacuum Poioidai Flux Functions Satisfying the Grad-Shafranov Equation in the Fiat-Ring Cyclide Coordinate System. Bullet. 94, Hokkaido University, Facul. of Engineering.Google Scholar
Lao, L. L., John, H. S., Stambaugh, R., Kellman, A. & Pfeiffer, W. 1985 Reconstruction of current profile parameters and plasma shapes in tokamaks. Nucl. Fusion 25, 1611.Google Scholar
Lao, L. L., St. John, H. E., Peng, Q., Ferron, J. R., Strait, E. J., Taylor, T. S. et al. 2005 MHD equilibrium reconstruction in the DIII-D Tokamak. Fusion Sci. Technol. 48, 968.Google Scholar
Lao, L., Hirschmann, S. P. & Wieland, R. 1981 Variational moment solutions to the Grad Shafranov equation. Phys. Fluids 24, 1431.Google Scholar
Lebedev, N. 1937 The functions associated with a ring of oval cross-section. Tekh. Fiz. 4, 3.Google Scholar
Li, Y. G., Lotte, P., Zwingmann, W., Gil, C. & Imbeaux, F. 2011 EFIT equilibrium reconstruction including polarimetry measurements on tore supra. Fusion Sci. Technol. 59, 397.Google Scholar
Moon, P. & Spencer, D. 1971 Field Theory Handbook. Springer.Google Scholar
Morse, P. & Feshbach, H. 1953 Methods of Theoretical Physics. MacGraw-Hill.Google Scholar
Mukhovatov, V. S. & Shafranov, V. D. 1971 Plasma equilibrium in a Tokamak. Nucl. Fusion 11, 605.Google Scholar
Neumann, C.1864.Theorie der Elektrizitats-und Warmever-teilung in einem Ringe.Google Scholar
Sartori, F., Ambrosino, G., Ariola, M., Cenedese, A., Crisanti, F., De Tommasi, G. et al. 2005 The system architecture of the new JET shape controller. Fusion Engin. Des. 74, 587.Google Scholar
Sartori, F., Crisanti, F., Albanese, R., Ambrosino, G., Toigo, V., Hay, J. et al. 2008 The JET PCU project: an international plasma control project. Fusion Engin. Des. 83, 202.Google Scholar
Shafranov, V. D. 1960 Equilibrium of a plasma toroid in a magnetic field. Sov. Phys. JETP 37, 775.Google Scholar
Shafranov, V. 1971 Determination of the parameters $\unicode[STIX]{x1D6FD}p$ and li in a Tokamak for arbitrary shape of plasma pinch cross-section. Plasma Phys. 13, 757.Google Scholar
Soloviev, L. S. 1968 The theory of hydromagnetic stability of toroidal plasma configurations. Sov. Phys. 26, 400.Google Scholar
Van Milligen, B. P. 1990 Exact relations between multipole moments of the flux and moments of the toroidal current density in Tokamaks. Nucl. Fusion 30, 157.Google Scholar
Van Milligen, B. P. & Lopez Fraguas, A. 1994 Expansion of vacuum magnetic fields in toroidal harmonics. Comput. Phys. Commun. 81, 74.Google Scholar
Wangerin, A.1909.Theorie des Potentials und der Kugelfunktionen. Leipzig: G. J. Gfischen’sche Verlagshandlung.Google Scholar
Yue, X. N., Xiao, B. J., Luo, Z. P. & Guo, Y. 2013 Plasma physics and controlled fusion paper fast equilibrium reconstruction for Tokamak discharge control based on GPU. Plasma Phys. Control. Fusion 55, 085016.Google Scholar
Zheng, S. B., Wootton, A. J. & Solano, E. R. 1996 Analytical tokamak equilibrium for shaped plasmas. Phys. Plasmas 3, 1176.Google Scholar
Figure 0

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

Figure 1

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.