Simsopt provides several representations of magnetic fields, including the Biot-Savart field associated with coils, as well as others. Magnetic fields are represented as subclasses of the base class :obj:`simsopt.field.magneticfield.MagneticField`. The various field types can be scaled and summed together; for instance you can add a purely toroidal field to the field from coils. Each field object has an associated set of evaluation points. At these points, you can query the field B, the vector potential A, and their derivatives with respect to position or the field parameters.
In simsopt, a filamentary coil is represented by the class
:obj:`simsopt.field.coil.Coil`. A Coil is a pairing of a
:obj:`~simsopt.geo.curve.Curve` with a current magnitude. The latter
is represented by a :obj:`simsopt.field.coil.Current` object. For
information about Curve objects see :ref:`the page on geometric
objects <curves>`. If you wish for several filamentary curves to carry
the same current, you can reuse a :obj:`~simsopt.field.coil.Current`
object in multiple :obj:`~simsopt.field.coil.Coil` objects.
Once :obj:`~simsopt.field.coil.Coil` objects are defined, the magnetic field they produce can be evaluated using the class :obj:`simsopt.field.biotsavart.BiotSavart`. This class provides an implementation of the Biot-Savart law
B(\mathbf{x}) = \frac{\mu_0}{4\pi} \sum_{k=1}^{n_\mathrm{coils}} I_k \int_0^1 \frac{(\Gamma_k(\phi)-\mathbf{x})\times \Gamma_k'(\phi)}{\|\Gamma_k(\phi)-\mathbf{x}\|^3} d\phi
where \mu_0=4\pi \times 10^{-7} is the vacuum permitivity, \Gamma_k is the position vector of coil k, and I_k indicates the electric currents.
Example:
import numpy as np from simsopt.geo.curvexyzfourier import CurveXYZFourier from simsopt.field.coil import Current, Coil from simsopt.field.biotsavart import BiotSavart curve = CurveXYZFourier(100, 1) # 100 = Number of quadrature points, 1 = max Fourier mode number curve.x = [0, 0, 1., 0., 1., 0., 0., 0., 0.] # Set Fourier amplitudes coil = Coil(curve, Current(1.0e4)) # 10 kAmpere-turns field = BiotSavart([coil]) # Multiple coils can be included in the list field.set_points(np.array([[0.5, 0.5, 0.1], [0.1, 0.1, -0.3]])) print(field.B())
For a more complex example of a
:obj:`~simsopt.field.biotsavart.BiotSavart` object used in coil
optimization, see
examples/2_Intermediate/stage_two_optimization.py.
The :obj:`simsopt.field.magneticfieldclasses.ToroidalField` class represents a purely toroidal magnetic field. The field is given by \mathbf B = B_0 \frac{R_0}{R} \mathbf e_\phi, where R_0 and B_0 are input scalar quantities, with R_0 representing a reference major radius, B_0 the magnetic field at R_0, and R is the radial coordinate of the cylindrical coordinate system (R,Z,\phi). The vector \mathbf e_\phi is a unit vector pointing in the direction of increasing \phi, with \phi the standard azimuthal angle. Given Cartesian coordinates (x,y,z), \mathbf e_\phi is calculated as \mathbf e_\phi=-\sin \phi \mathbf e_x+\cos \phi \mathbf e_y with \phi=\arctan(y/x).
The :obj:`simsopt.field.magneticfieldclasses.PoloidalField` class represents a poloidal magnetic field acording to the formula \mathbf B = B_0 \frac{r}{q R_0} \mathbf e_\theta, where R_0, q and B_0 are input scalar quantities. R_0 represents the major radius of the magnetic axis, B_0 the magnetic field at r=R_0 q and q the safety factor associated with the sum of a poloidal magnetic field and a toroidal magnetic field with major radius R_0 and magnetic field on-axis B_0. r is the radial coordinate of the simple toroidal coordinate system (r,\phi,\theta). Given a set of points (x,y,z), r is calculated as r=\sqrt{(\sqrt{x^2+y^2}-R_0)^2+z^2}. The vector \mathbf e_\theta is a unit vector pointing in the direction of increasing \theta, with \theta the poloidal angle in the simple toroidal coordinate system (r,\phi,\theta). Given a set of points (x,y,z), \mathbf e_\theta is calculated as \mathbf e_\theta=-\sin \theta \cos \phi \mathbf e_x+\sin \theta \sin \phi \mathbf e_y+\cos \theta \mathbf e_z with \phi=\arctan(y/x) and \theta=\arctan(z/(\sqrt{x^2+y^2}-R_0)).
The
:obj:`simsopt.field.magneticfieldclasses.ScalarPotentialRZMagneticField`
class initializes a vacuum magnetic field \mathbf B = \nabla
\Phi defined via a scalar potential \Phi in cylindrical
coordinates (R,Z,\phi). The field \Phi is specified as
an analytical expression via a string argument. Simsopt performs the
necessary partial derivatives in order find \mathbf B and its
derivatives. For example, the function
ScalarPotentialRZMagneticField("2*phi") represents a toroidal
magnetic field \mathbf B = \nabla (2\phi)=2/R \mathbf e_\phi.
Note: this functions needs the library sympy for the analytical
derivatives.
The :obj:`simsopt.field.magneticfieldclasses.CircularCoil` class represents a magnetic field created by a single circular coil. It takes four input quantities: a, the radius of the coil, \mathbf c=[c_x,c_y,c_z], the center of the coil, I, the current flowing through the coil and \mathbf n, the normal vector to the plane of the coil centered at the coil radius, which could be specified either with its three Cartesian components \mathbf n=[n_x,n_y,n_z] or as \mathbf n=[\theta,\phi] with the spherical angles \theta and \phi.
The magnetic field is calculated analitically using the following expressions (reference)
- B_x=\frac{\mu_0 I}{2\pi}\frac{x z}{\alpha^2 \beta \rho^2}\left[(a^2+r^2)E(k^2)-\alpha^2 K(k^2)\right]
- B_y=\frac{y}{x}B_x
- B_z=\frac{\mu_0 I}{2\pi \alpha^2 \beta}\left[(a^2-r^2)E(k^2)+\alpha^2 K(k^2)\right]
where \rho^2=x^2+y^2, r^2=x^2+y^2+z^2, \alpha^2=a^2+r^2-2a\rho, \beta^2=a^2+r^2+2 a \rho, k^2=1-\alpha^2/\beta^2.
The :obj:`simsopt.field.magneticfieldclasses.Dommaschk` class represents a vacuum magnetic field \mathbf B = \nabla \Phi with basis functions for the scalar potential \Phi described in W. Dommaschk (1986), Computer Physics Communications 40, 203-218. This representation provides explicit analytic formulae for vacuum fields with a mixture of flux surfaces, islands, and chaos. Following the original reference, a toroidal field with B_0=R_0=1 is already included in the definition. As input parameters, it takes two arrays:
- The first array is an N\times2 array [(m_1,n_1),(m_2,n_2),...] specifying which harmonic coefficients m and n are non-zero.
- The second array is an N\times2 array [(b_1,c_1),(b_2,c_2),...] with b_i=b_{m_i,n_i} and c_i=c_{m_i,n_i} the coefficients used in the Dommaschk representation.
The :obj:`simsopt.field.magneticfieldclasses.Reiman` provides the magnetic field model in section 5 of Reiman and Greenside, Computer Physics Communications 43 (1986) 157—167. It is an analytical magnetic field representation that allows the explicit calculation of the width of the magnetic field islands.
The :obj:`simsopt.field.magneticfieldclasses.InterpolatedField` function takes an existing field and interpolates it on a regular grid in r,\phi,z. This resulting interpolant can then be evaluated very quickly. This is useful for efficiently tracing field lines and particle trajectories.
Magnetic field objects can be added together, either by using the
+ operator, or by creating an instance of the class
:obj:`simsopt.field.magneticfield.MagneticFieldSum`. (The +
operator creates the latter.)
Magnetic fields can also be scaled by a constant. This can be accomplished either using the * operator,
or by creating an instance of the class
:obj:`simsopt.field.magneticfield.MagneticFieldMultiply`. (The *
operator creates the latter.)
Example:
from simsopt.field.magneticfieldclasses import ToroidalField, CircularCoil field1 = CircularCoil(I=1.e7, r0=1.) field2 = ToroidalField(R0=1., B0=1.) total_field = field1 + 2.5 * field2
Magnetic field objects have a large number of functions available. Before evaluating the field, you must
set the evaluation points. This can be done using either Cartesian or cylindrical coordinates.
Let m be a :obj:`~simsopt.field.magneticfield.MagneticField` object, and suppose there are n points
at which you wish to evaluate the field.
m.set_points_cart()takes a numpy array of size(n, 3)with the Cartesian coordinates(x, y, z)of the points.m.set_points_cyl()takes a numpy array of size(n, 3)with the cylindrical coordinates(r, phi, z)of the points.m.set_points()is shorthand form.set_points_cart().m.get_points_cart()returns a numpy array of size(n, 3)with the Cartesian coordinates(x, y, z)of the points.m.get_points_cyl()returns a numpy array of size(n, 3)with the cylindrical coordinates(r, phi, z)of the points.
A variety of functions are available to return the magnetic field B, vector potential A, and their gradients. The most commonly used ones are the following:
m.B()returns an array of size(n, 3)with the Cartesian coordinates of B.m.B_cyl()returns an array of size(n, 3)with the cylindrical(r, phi, z)coordinates of B.m.A()returns an array of size(n, 3)with the Cartesian coordinates of A.m.AbsB()returns an array of size(n, 1)with the field magnitude |B|.m.dB_by_dX()returns an array of size(n, 3, 3)with the Cartesian coordinates of \nabla B. Denoting the indices by (i,j,l), the result contains \partial_j B_l(x_i).m.d2B_by_dXdX()returns an array of size(n, 3, 3, 3)with the Cartesian coordinates of \nabla\nabla B. Denoting the indices by (i,j,k,l), the result contains \partial_k \partial_j B_l(x_i).m.dA_by_dX()returns an array of size(n, 3, 3)with the Cartesian coordinates of \nabla A. Denoting the indices by (i,j,l), the result contains \partial_j A_l(x_i).m.d2A_by_dXdX()returns an array of size(n, 3, 3, 3)with the Cartesian coordinates of \nabla\nabla A. Denoting the indices by (i,j,k,l), the result contains \partial_k \partial_j A_l(x_i).m.GradAbsB()returns an array of size(n, 3)with the Cartesian components of \nabla |B|.
Example:
import numpy as np from simsopt.field.magneticfieldclasses import CircularCoil field = CircularCoil(I=1.e7, r0=1.) points = np.array([[0.5, 0.5, 0.1], [0.1, 0.1, -0.3]]) field.set_points(points) print(field.B()) print(field.dB_by_dX())