I would like to plot the magnetic field lines of a tokamak.
Introduction
Since the magnetic field is solenoidal: $\nabla \cdot \vec{B} = 0$, it is possible to decompose it as:
$\vec{B} = \vec{B}_{T} + \vec{B}_{P} \tag{1}$
where $\vec{B}_{T}$ is the Toroidal magnetic field and $\vec{B}_{P}$ is the Poloidal magnetic field. This vector decomposition is called Toroidal-Poloidal decomposition $[1]$.
Now, in toroidal fusion reactors, or Tokamaks, the intuitive picture of the magnetic field lines of $\vec{B}$, $\vec{B}_{T}$ and $\vec{B}_{P}$ is given by Figure $1$:
Figure $1$: The magnetic fields of a Tokamak system.
My question
An illustrative example:
Well, if you want to plot the magnetic field lines of a magnetic dipole, you simply take the field $[2]$:
$\vec{B}_{\mathrm{dip}} = \frac{\mu_{0}}{4\pi r^{3}}2cos\theta \hat{r} + \frac{\mu_{0}}{4\pi r^{3}}sin\theta \hat{\theta} \tag{2}$
and perform a coordinate transformation $(r,\theta,\phi) \to (x,y,z)$ to get $[2]$, $[3]$:
$\vec{B}_{\mathrm{dip}} = \frac{3 \mu m x z}{4 \pi \left(x^2+y^2+z^2\right)^{5/2}}\hat{x}+\frac{3 \mu m y z}{4 \pi \left(x^2+y^2+z^2\right)^{5/2}}\hat{y}-\frac{\mu m \left(x^2+y^2-2 z^2\right)}{4 \pi \left(x^2+y^2+z^2\right)^{5/2}} \hat{z} \tag{3}$
with $(3)$ you can use your favorite software to plot the magnetic configuration. Using mathematica you should get (Figure $2$):
Figure $2$: The magnetic field lines of theoretical magnetic dipole.
For the present case,
The question:
The magnetic field of a tokamak is far more complicated than that of a magnetic dipole. Assuming an axisymmetric ideal MHD equilibrium leads to the Grad-Shafranov equation. This can be solved numerically for a specified toroidal current density, $\mathbf{j}_{\phi}$. Assuming the current density is linearly related to the poloidal flux yields the Solov’ev solution. For large aspect ratio, Wesson's Tokamaks provides the approximate large aspect ratio solution:
$B_{\phi} = \frac{B_{\phi 0} R_0}{R} \tag{4}$
$B_{\theta} = B_{\theta a} \frac{1 - \left (1 - \frac{r^2}{a^2} \right )^{\nu+1}}{r / a} \tag{5}$
where its assumed current $j \propto \left (1 - \frac{r^2}{a^2} \right )^{\nu}$ and pressure $p \propto 1 - \frac{r^2}{a^2} $. Here the simple toroidal coordinates are used. Also, $B_{\phi0}$ is the toroidal magnetic field at $R=R_{0}$ (center of the cross section) and $B_{\theta a}$ is the poloidal magnetic field at $r=a$ (edge of plasma). $R$ and $R_{0}$ are respectively distance from axis of symmetry and torus major radius ( simple toroidal coordinates ).
The code that I've written is:
Manipulate[
StreamPlot3D[{0, (c*f)/Sqrt[x^2 + y^2], (
a (1 - (1 - Sqrt[x^2 + y^2 + z^2]/a^2)^(1 + ν)) *d)/
Sqrt[x^2 + y^2 + z^2]}, {x, -2, 2}, {y, -2, 2}, {z, -2, 2}], {a, 1,
5}, {ν, 1, 5}, {d, 1, 5}, {c, 1, 5}, {f, 1, 5}]
But I simply do not get the right results.
So: how can I plot the magnetic field lines of $(4)$ and $(5)$ using the StreamPlot3D?
$$--------------------$$
$[1]$ https://en.wikipedia.org/wiki/Poloidal%E2%80%93toroidal_decomposition
$[2]$ GRIFFITHS.D. Introduction to Electrodynamics. page 255.