Concepts

Ways to use simsopt

Simsopt is a collection of classes and functions that can be used in several ways. One application is to solve optimization problems involving stellarators, similar to STELLOPT. You could also define an objective function using simsopt, but use an optimization library outside simsopt to solve it. Or, you could use the simsopt optimization infrastructure to optimize your own objects, which may or may not have any connection to stellarators. Alternatively, you can use the stellarator-related objects in a script for some purpose other than optimization, such as plotting how some code output varies as an input parameter changes, or evaluating the finite-difference gradient of some code outputs. Or, you can manipulate the objects interactively, at the python command line or in a Jupyter notebook.

Input files

Simsopt does not use input data files to define optimization problems, in contrast to STELLOPT. Rather, problems are specified using a python driver script, in which objects are defined and configured. However, objects related to specific physics codes may use their own input files. In particular, a simsopt.mhd.vmec.Vmec object can be initialized using a standard VMEC input.* input file, and a simsopt.mhd.spec.Spec object can be initialized using a standard SPEC *.sp input file.

Optimization

To do optimization using simsopt, there are four basic steps:

  1. Create some objects that will participate in the optimization.

  2. Define the independent variables for the optimization, by choosing which degrees of freedom of the objects are free vs fixed.

  3. Define an objective function.

  4. Solve the optimization problem that has been defined.

This pattern is evident in the examples in this documentation and in the examples directory of the repository.

Some typical objects are a MHD equilibrium represented by the VMEC or SPEC code, or some electromagnetic coils. To define objective functions, a variety of additional objects can be defined that depend on the MHD equilibrium or coils, such as a simsopt.mhd.boozer.Boozer object for Boozer-coordinate transformation, a simsopt.mhd.spec.Residue object to represent Greene’s residue of a magnetic island, or a simsopt.geo.objectives.LpCurveCurvature penalty on coil curvature.

More details about setting degrees of freedom and defining objective functions can be found on the Defining optimization problems page.

For the solution step, two functions are provided presently, simsopt.solve.serial.least_squares_serial_solve() and simsopt.solve.mpi.least_squares_mpi_solve(). The first is simpler, while the second allows MPI-parallelized finite differences to be used in the optimization.

MPI Partitions and worker groups

To use MPI parallelization with simsopt, you must install the mpi4py package.

MPI simsopt programs can generally be launched using mpirun or mpiexec, e.g.

mpiexec -n 32 ./name_of_your_python_driver_script

where 32 can be replaced by the number of processors you want to use. There may be specific instructions to run MPI programs on your HPC system, so consult the documentation for your system.

Given a set of \(M\) MPI processes, an optimization problem can be parallelized in different ways. At one extreme, all \(M\) processes could work together on each evaluation of the objective function, with the optimization algorithm itself only initiating a single evaluation of the objective function \(f\) at a time. At the other extreme, the optimization algorithm could request \(M\) evaluations of \(f\) all at once, with each evaluation performed by a single processor. This type of parallelization can be useful for instance when finite differences are used to evaluate gradients. Or, both types of parallelization could be used at the same time. To allow all of these types of parallelization, simsopt uses the concepts of worker groups and MPI partitions.

A worker group is a set of MPI processes that works together to evaluate the objective function. An MPI partition is a division of the total set of \(M\) MPI processors into one or more worker groups. If there are \(W\) worker groups, each group has approximately \(M/W\) processors (approximate because one may need to round up or down.) Simsopt has a class simsopt.util.mpi.MpiPartition to manage this division of processors into worker groups. Each MPI-aware simsopt object, such as simsopt.mhd.vmec.Vmec, keeps an instance of MpiPartition which can be set from its constructor. The number of worker groups can be set by an argument to MpiPartition. For instance, to tell a Vmec object to use four worker groups, one could write

from simsopt.mhd import Vmec
from simsopt.util.mpi import MpiPartition

mpi = MpiPartition(4)
equil = Vmec("input.li383_low_res", mpi=mpi)

The same MpiPartition instance should be passed to the solver:

# ... code to define an optimization problem "prob" ...

from simsopt.solve.mpi import least_squares_mpi_solve

least_squares_mpi_solve(prob, mpi, grad=True)

Many optimization algorithms that do not use derivatives do not support concurrent evaluations of the objective. In this case, the number of worker groups should be \(W=1\). Any algorithm that uses derivatives, such as Levenberg-Marquardt, can take advantage of multiple worker groups to evaluate derivatives by finite differences. If the number of parameters (i.e. independent variables) is \(N\), you ideally want to set \(W=N+1\) when using 1-sided finite differences, and set \(W=2N+1\) when using centered differences. These ideal values are not required however - simsopt will evaluate finite difference derivatives for any value of \(W\), and results should be exactly independent of \(W\). Other derivative-free algorithms intrinsically support parallelization, such as HOPSPACK, though no such algorithm is available in simsopt yet.

An MPI partition is associated with three MPI communicators, “world”, “groups”, and “leaders”. The “world” communicator represents all \(M\) MPI processors available to the program. (Normally this communicator is the same as MPI_COMM_WORLD, but it could be a subset thereof if you wish.) The “groups” communicator also contains the same \(M\) processors, but grouped into colors, with a different color representing each worker group. Therefore operations such as MPI_Send and MPI_Bcast on this communicator exchange data only within one worker group. This “groups” communicator is therefore the one that must be used for evaluation of the objective function. Finally, the “leaders” communicator only includes the \(W\) processors that have rank 0 in the “groups” communicator. This communicator is used for communicating data within a parallel optimization algorithm (as opposed to within a parallelized objective function).

Given an instance of simsopt.util.mpi.MpiPartition named mpi, the number of worker groups is available as mpi.ngroups, and the index of a given processor’s group is mpi.group. The communicators are available as mpi.comm_world, mpi.comm_groups, and mpi.comm_leaders. The number of processors within the communicators can be determined from mpi.nprocs_world, mpi.nprocs_groups, and mpi.nprocs_leaders. The rank of the present processor within the communicators is available as mpi.rank_world, mpi.rank_groups, and mpi.rank_leaders. To determine if the present processor has rank 0 in a communicator, you can use the variables mpi.proc0_world or mpi.proc0_groups, which have type bool.

Geometric Objects

Simsopt contains implementations of two commonly used geometric objects: curves and surfaces.

Curves

A simsopt.geo.curve.Curve is modelled as a function \(\Gamma:[0, 1] \to \mathbb{R}^3\). A curve object stores a list of \(n_\phi\) quadrature points \(\{\phi_1, \ldots, \phi_{n_\phi}\} \subset [0, 1]\) and returns all information about the curve, at these quadrature points.

  • Curve.gamma(): returns a (n_phi, 3) array containing \(\Gamma(\phi_i)\) for \(i\in\{1, \ldots, n_\phi\}\), i.e. returns a list of XYZ coordinates along the curve.

  • Curve.gammadash(): returns a (n_phi, 3) array containing \(\Gamma'(\phi_i)\) for \(i\in\{1, \ldots, n_\phi\}\), i.e. returns the tangent along the curve.

  • Curve.kappa(): returns a (n_phi, 1) array containing the curvature \(\kappa\) of the curve at the quadrature points.

  • Curve.torsion(): returns a (n_phi, 1) array containing the torsion \(\tau\) of the curve at the quadrature points.

The different curve classes, such as simsopt.geo.curverzfourier.CurveRZFourier and simsopt.geo.curvexyzfourier.CurveXYZFourier differ in the way curves are discretized. Each of these take an array of dofs (e.g. Fourier coefficients) and turn these into a function \(\Gamma:[0, 1] \to \mathbb{R}^3\). These dofs can be set via Curve.set_dofs(dofs) and dofs = Curve.get_dofs(), where n_dofs = dofs.size. Changing dofs will change the shape of the curve. Simsopt is able to compute derivatives of all relevant quantities with respect to the discretisation parameters. For example, to compute the derivative of the coordinates of the curve at quadrature points, one calls Curve.dgamma_by_dcoeff(). One obtains a numpy array of shape (n_phi, 3, n_dofs), containing the derivative of the position at every quadrature point with respect to every degree of freedom of the curve. In the same way one can compute the derivative of quantities such as curvature (via Curve.dkappa_by_dcoeff()) or torsion (via Curve.dtorsion_by_dcoeff()).

A number of quantities are implemented in simsopt.geo.curveobjectives and are computed on a simsopt.geo.curve.Curve:

  • CurveLength: computes the length of the Curve.

  • LpCurveCurvature: computes a penalty based on the \(L_p\) norm of the curvature on a curve.

  • LpCurveTorsion: computes a penalty based on the \(L_p\) norm of the torsion on a curve.

  • MinimumDistance: computes a penalty term on the minimum distance between a set of curves.

The value of the quantity and its derivative with respect to the curve dofs can be obtained by calling e.g., CurveLength.J() and CurveLength.dJ().

Surfaces

The second large class of geometric objects are surfaces. A simsopt.geo.surface.Surface is modelled as a function \(\Gamma:[0, 1] \times [0, 1] \to \mathbb{R}^3\) and is evaluated at quadrature points \(\{\phi_1, \ldots, \phi_{n_\phi}\}\times\{\theta_1, \ldots, \theta_{n_\theta}\}\). Here, \(\phi\) is the toroidal angle and \(\theta\) is the poloidal angle. Note that \(\phi\) and \(\theta\) go up to 1, not up to \(2 \pi\)!

In practice, you almost never use the base Surface class. Rather, you typically use one of the subclasses corresponding to a specific parameterization. Presently, the available subclasses are SurfaceRZFourier, SurfaceGarabedian, SurfaceHenneberg, SurfaceXYZFourier, and SurfaceXYZTensorFourier. In many cases you can convert a surface from one type to another by going through SurfaceRZFourier, as most surface types have to_RZFourier() and from_RZFourier() methods. Note that SurfaceRZFourier corresponds to the surface parameterization used internally in the VMEC and SPEC codes. However when using these codes in simsopt, any of the available surface subclasses can be used to represent the surfaces, and simsopt will automatically handle the conversion to SurfaceRZFourier when running the code.

The points \(\phi_j\) and \(\theta_j\) are used for evaluating the position vector and its derivatives, for computing integrals, and for plotting, and there are several available methods to specify these points. For \(\theta_j\), you typically specify a keyword argument ntheta to the constructor when instantiating a surface class. This results in a grid of ntheta uniformly spaced points between 0 and 1, with no endpoint at 1. Alternatively, you can specify a list or array of points to the quadpoints_theta keyword argument when instantiating a surface class, specifying the \(\theta_j\) directly. If both ntheta and quadpoints_theta are specified, an exception will be raised. For the \(\phi\) coordinate, you sometimes want points up to 1 (the full torus), sometimes up to \(1/n_{fp}\) (one field period), and sometimes up to \(1/(2 n_{fp})\) (half a field period). These three cases can be selected by setting the range keyword argument of the surface subclasses to "full torus", "field period", or "half period". Equivalently, you can set range to the constants S.RANGE_FULL_TORUS, S.RANGE_FIELD_PERIOD, or S.RANGE_HALF_PERIOD, where S can be simsopt.geo.surface.Surface or any of its subclasses. For all of these cases, the nphi keyword argument can be set to the desired number of \(\phi\) grid points. Alternatively, you can pass a list or array to the quadpoints_phi keyword argument of the constructor for any Surface subclass to specify the \(\phi_j\) points directly. An exception will be raised if both nphi and quadpoints_phi are specified. For more information about these arguments, see the SurfaceRZFourier API documentation.

The methods available to each surface class are similar to those of the Curve class:

  • Surface.gamma(): returns a (n_phi, n_theta, 3) array containing \(\Gamma(\phi_i, \theta_j)\) for \(i\in\{1, \ldots, n_\phi\}, j\in\{1, \ldots, n_\theta\}\), i.e. returns a list of XYZ coordinates on the surface.

  • Surface.gammadash1(): returns a (n_phi, n_theta, 3) array containing \(\partial_\phi \Gamma(\phi_i, \theta_j)\) for \(i\in\{1, \ldots, n_\phi\}, j\in\{1, \ldots, n_\theta\}\).

  • Surface.gammadash2(): returns a (n_phi, n_theta, 3) array containing \(\partial_\theta \Gamma(\phi_i, \theta_j)\) for \(i\in\{1, \ldots, n_\phi\}, j\in\{1, \ldots, n_\theta\}\).

  • Surface.normal(): returns a (n_phi, n_theta, 3) array containing \(\partial_\phi \Gamma(\phi_i, \theta_j)\times \partial_\theta \Gamma(\phi_i, \theta_j)\) for \(i\in\{1, \ldots, n_\phi\}, j\in\{1, \ldots, n_\theta\}\).

  • Surface.area(): returns the surface area.

  • Surface.volume(): returns the volume enclosed by the surface.

A number of quantities are implemented in simsopt.geo.surfaceobjectives and are computed on a simsopt.geo.surface.Surface:

  • ToroidalFlux: computes the flux through a toroidal cross section of a Surface.

The value of the quantity and its derivative with respect to the surface dofs can be obtained by calling e.g., ToroidalFlux.J() and ToroidalFlux.dJ_dsurfacecoefficients().

Caching

The quantities that Simsopt can compute for curves and surfaces often depend on each other. For example, the curvature or torsion of a curve both rely on Curve.gammadash(); to avoid repeated calculation geometric objects contain a cache that is automatically managed. If a quantity for the curve is requested, the cache is checked to see whether it was already computed. This cache can be cleared manually by calling Curve.invalidate_cache(). This function is called everytime Curve.set_dofs() is called (and the shape of the curve changes).

Magnetic Field Classes

Simsopt contains several magnetic field classes available to be called directly. Any field can be summed with any other field and/or multiplied by a constant parameter. To get the magnetic field (or its derivatives) at a set of points, first, an instance of that particular magnetic field is created, then all its properties are evaluated internally at those points and, finally, those properties can be outputed. Below is an example that prints the components of a magnetic field and its derivatives of a sum of a circular coil in the xy-plane with current I=1.e7 and a radius r0=1 and a toroidal field with a magnetic field B0=1 at major radius R0=1. This field is evaluated at the set of points=[[0.5, 0.5, 0.1],[0.1, 0.1, -0.3]].

from simsopt.field.magneticfieldclasses import ToroidalField, CircularCoil

Bfield1 = CircularCoil(I=1.e7, r0=1.)
Bfield2 = ToroidalField(R0=1., B0=1.)
Bfield = Bfield1 + Bfield2
points=[[0.5, 0.5, 0.1], [0.1, 0.1, -0.3]]
Bfield.set_points(points)
print(Bfield.B())
print(Bfield.dB_by_dX())

Below is a similar example where, instead of calculating the magnetic field using analytical functions from the circular coil class, it is calculated using the BiotSavart class

from simsopt.field.magneticfieldclasses import ToroidalField
from simsopt.field.biotsavart import BiotSavart
from simsopt.geo.curvexyzfourier import CurveXYZFourier

coil = CurveXYZFourier(300, 1)
coil.set_dofs([0, 0, 1., 0., 1., 0., 0., 0., 0.])
Bfield1 = BiotSavart([coil], [1.e7])
Bfield2 = ToroidalField(R0=1., B0=1.)
Bfield = Bfield1 + Bfield2
points=[[0.5, 0.5, 0.1], [0.1, 0.1, -0.3]]
Bfield.set_points(points)
print(Bfield.B())
print(Bfield.dB_by_dX())

BiotSavart

The simsopt.field.biotsavart.BiotSavart class initializes a magnetic field vector induced by a list of closed curves \(\Gamma_k\) with electric currents \(I_k\). The field is given by

\[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 10^{-7}\) is the vacuum permitivity. As input, it takes a list of closed curves and the corresponding currents.

ToroidalField

The simsopt.field.magneticfieldclasses.ToroidalField class initializes a toroidal magnetic field vector acording to the formula \(\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 the major radius of the magnetic axis and \(B_0\) the magnetic field at \(R_0\). \(R\) is the radial coordinate of the cylindrical coordinate system \((R,Z,\phi)\), so that \(B\) has the expected \(1/R\) dependence. Given a set of points \((x,y,z)\), \(R\) is calculated as \(R=\sqrt{x^2+y^2}\). The vector \(\mathbf e_\phi\) is a unit vector pointing in the direction of increasing \(\phi\), with \(\phi\) the standard azimuthal angle of the cylindrical coordinate system \((R,Z,\phi)\). Given a set of points \((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)\).

PoloidalField

The simsopt.field.magneticfieldclasses.PoloidalField class initializes a poloidal magnetic field vector 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))\).

ScalarPotentialRZMagneticField

The 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 given as an analytical expression and simsopt performed the necessary partial derivatives in order find \(\mathbf B\) and its derivatives. Example: the function

ScalarPotentialRZMagneticField("2*phi")

initializes 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.

CircularCoil

The simsopt.field.magneticfieldclasses.CircularCoil class initializes 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\).

Dommaschk

The simsopt.field.magneticfieldclasses.Dommaschk class initializes a vacuum magnetic field \(\mathbf B = \nabla \Phi\) with a representation for the scalar potential \(\Phi\) as proposed in W. Dommaschk (1986), Computer Physics Communications 40, 203-218. It allows to quickly generate magnetic fields with islands with only a small set of scalar quantities. 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.

Reiman

The simsopt.field.magneticfieldclasses.Reiman initializes 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. It takes as input arguments: \(\iota_0\), the unperturbed rotational transform, \(\iota_1\), the unperturbed global magnetic shear, \(k\), an array of integers with that specifies the Fourier modes used, \(\epsilon_k\), an array that specifies the coefficient in front of the Fourier modes, \(m_0\), the toroidal symmetry parameter (usually 1).

InterpolatedField

The 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. As input arguments, it takes field: the underlying simsopt.field.magneticfield.MagneticField to be interpolated, degree: the degree of the piecewise polynomial interpolant, rrange: a 3-tuple of the form (rmin, rmax, nr), phirange: a 3-tuple of the form (phimin, phimax, nphi), zrange: a 3-tuple of the form (zmin, zmax, nz), extrapolate: whether to extrapolate the field when evaluate outside the integration domain or to throw an error, nfp: Whether to exploit rotational symmetry, stellsym: Whether to exploit stellarator symmetry.

Particle Tracing

Simsopt is able to follow particles in a magnetic field. The main function to use in this case is simsopt.field.tracing.trace_particles (click the link for more information on the input and output parameters) and it is able to use two different sets of equations depending on the input parameter mode:

  • In the case of mode='full' it solves

\[[\ddot x, \ddot y, \ddot z] = \frac{q}{m} [\dot x, \dot y, \dot z] \times \mathbf B\]

  • In the case of mode='gc_vac' it solves the guiding center equations under the assumption \(\nabla p=0\), that is

\[\begin{split}[\dot x, \dot y, \dot z] &= v_{||}\frac{\mathbf B}{B} + \frac{m}{q|B|^3} \left(\frac{v_\perp^2}{2} + v_{||}^2\right) \mathbf B\times \nabla B\\ \dot v_{||} &= -\mu \mathbf B \cdot \nabla B\end{split}\]

where \(v_\perp^2 = 2\mu B\). See equations (12) and (13) of Guiding Center Motion, H.J. de Blank.

Below is an example of the vertical drift experienced by two particles in a simple toroidal magnetic field.

from simsopt.field.magneticfieldclasses import ToroidalField
from simsopt.geo.curvexyzfourier import CurveXYZFourier
from simsopt.util.constants import PROTON_MASS, ELEMENTARY_CHARGE, ONE_EV
from simsopt.field.tracing import trace_particles_starting_on_curve

bfield = ToroidalField(B0=1., R0=1.)
start_curve = CurveXYZFourier(300, 1)
start_curve.set_dofs([0, 0, 1.01, 0, 1.01, 0., 0, 0., 0.])
nparticles = 2
m = PROTON_MASS
q = ELEMENTARY_CHARGE
tmax = 1e-6
Ekin = 100*ONE_EV
gc_tys, gc_phi_hits = trace_particles_starting_on_curve(
    start_curve, bfield, nparticles, tmax=tmax, seed=1, mass=m, charge=q,
    Ekin=Ekin, umin=0.2, umax=0.5, phis=[], mode='gc_vac', tol=1e-11)
z_particle_1 = gc_tys[0][:][2]
z_particle_2 = gc_tys[1][:][2]
print(z_particle_1)
print(z_particle_2)

We note that SIMSOPT draws initial data for particles consisting of the guiding center position, the parallel, and the perpendicular speed.

  • To compute the speed, the user specifies the total kinetic energy and an interval of pitch-angles \([u_\min, u_\max]\). Given a pitch angle \(u\) and total speed \(v\) (computed from the kinetic energy and the particle mass), the parallel speed is given by \(v_{||} = u v\) and, and \(v_\perp = \sqrt{v^2-v_{||}^2}\).

  • In the case of full orbit simulations, we need velocity initial data and hence this only defines the initial data up to the phase. To specify the phase, one can pass the phase_angle variable to the tracing functions. A value in \([0, 2\pi]\) is expected.