Computer Physics Communications 294 (2024) 108930

Contents lists available at ScienceDirect

Computer Physics Communications
journal homepage: www.elsevier.com/locate/cpc

Computational Physics

A framework for synthetic diagnostics using energetic-particle orbits intokamaks ✩H. Järleblad a,∗ , L. Stagner b , M. Salewski a , J. Eriksson c , M. Nocente d , B.S. Schmidt a ,M. Rud Larsen aaDepartment of Physics, Technical University of Denmark, DK-2800 Kgs. Lyngby, Denmark
General Atomics, P.O. Box 85608, San Diego, CA 92186-5608, USAc
Department of Physics and Astronomy, Uppsala University, 751 20 Uppsala, SwedendDepartment of Physics, University of Milano-Bicocca, 20126 Milano, Italyb

a r t i c l e

i n f oArticle history:
Received 4 November 2022
Received in revised form 3 August 2023
Accepted 5 September 2023
Available online 9 September 2023
Keywords:
Nuclear fusion
Fast ionsOrbits
Weight functionsa b s t r a c t
In fusion plasma physics, the large-scale trajectories of energetic particles in magnetic conﬁnement
devices are known as orbits. To effectively and eﬃciently be able to work with orbits, the Orbit Weight
Computational Framework (OWCF) was developed. The OWCF constitutes a set of scripts, functions andapplications capable of computing, visualizing and working with quantities related to fast-ion (FI) orbits
in toroidally symmetric fusion devices. The current version is highly integrated with the DRESS code,
which enables the OWCF to compute and analyze the orbit sensitivity for arbitrary neutron- and gammadiagnostics. However, the framework is modular in the sense that any future codes (e.g. FIDASIM) canbe easily integrated. The OWCF can also compute projected velocity spectra for FI orbits, which play a
key role in many FI diagnostics. Via interactive applications, the OWCF can function both as a tool forinvestigative research but also for teaching. The OWCF will be used to analyze and simulate the diagnostic
results of current and future fusion experiments such as ITER. The orbit weight functions computedwith the OWCF can be used to reconstruct the FI distribution in terms of FI orbits from experimental
measurements using tomographic inversion.
© 2023 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY license(http://creativecommons.org/licenses/by/4.0/).

1. Introduction
In magnetic conﬁnement fusion devices, such as tokamaks, thepopulation of particles can be categorized into two sets; thermaland energetic [1]. In tokamaks of size similar to the Joint European Torus (JET, major radius ∼ 3 meters) [2], the following distinction can be made. If we take three times the thermal velocity
as a threshold, 3v th , the energetic fast ions are ∼ 10 times moreenergetic than the thermal (bulk) ions. For a typical plasma temperature of 1 − 10 keV (≈ 1 − 11 × 107 K), the lower limit forenergetic particles should arguably be placed around 10 keV. In
reality, the transition from thermal to fast is naturally continuous.
Despite their relatively small population [3], fast ions can cause
damage to ﬁrst-wall components [4] and drive plasma instabilities [5,6]. Fast ions also need to be conﬁned for a suﬃcient lengthof time (several slowing-down times, i.e. slowing down from fast✩

The review of this paper was arranged by Prof. David W. Walker.

* Corresponding author.

E-mail address: henrikj@dtu.dk (H. Järleblad).

to slow speeds) so as to heat the thermal ions and enable a selfsustained plasma burn [4].
Because of their relatively low energy, the thermal ions staymostly conﬁned to the magnetic ﬂux surfaces, while the fast ionshave more exotic large-scale motion [7,3]. This motion is commonly referred to as a fast-ion orbit (FIO). Depending on the energy,
starting position and angle of the particle velocity with respect
to the magnetic ﬁeld, the path of the FIOs can vary substantially.
Therefore, it is of great importance to know how sensitive a diagnostic is to different FIOs, since various diagnosable particles
or quanta (e.g. fusion-born neutrons and γ -rays) might originate
more strongly from certain FIOs [3]. Depending on the heating
scheme of the tokamak plasma, certain FIOs will be more populated than others [3] and the exact shape of the diagnostic signal
is thus dependent on a multitude of factors, including the magnetic equilibrium (which determines the shape of the topological
FIO regions in phase space) [8] and the heating scheme, as well as
the density and temperature proﬁles of the fusion plasma [9].
By knowing the sensitivity of a plasma diagnostic in phasespace, the diagnostic signal can be directly related to the fast-ion(FI) distribution, given the assumption of a linear relationship be-https://doi.org/10.1016/j.cpc.2023.108930
0010-4655/© 2023 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).

H. Järleblad, L. Stagner, M. Salewski et al.

Computer Physics Communications 294 (2024) 108930

Fig. 1. a) A basic building block of the OWCF. The script has some input data visualized as an ‘IN’ circle in the top-left corner. The script has some output data visualized as an‘OUT’ circle in the top-right corner. b) An illustration of script inter-dependency. Dotted circle means ‘optional’. Black dot means ‘depends on’, and the black line shows which
ﬁle that contains required functions. These schematics are similar to the schematics of e.g. the Uniﬁed Modeling Language [29], but with the inter-block arrow schematics
replaced by number ‘bubbles’ for clarity.

tween the signal and the FI distribution [10]. This can be expressedas the matrix equation

S=WF

FI distribution exist. These include expressing the FI distribution
in terms of a basis of slowing-down distribution functions [24,25].
For orbit-space grid resolutions of interest, 3D OWs usually need
to be computed using a computational cluster to be able to complete the calculation within a reasonable timeframe, due to both
RAM and CPU requirements. Naturally, to be able to calculate, analyze and work with these 3D OWs, a speed-optimized framework
written in a high-performance programming language is useful. To
this end, the orbit weight computational framework (OWCF) was
developed.
In short, the OWCF is a collection of scripts, utilities, structuresand functions written in the Julia programming language [26]. The
eﬃcient computation of guiding-center FIOs (where the ﬁnitenessof the Larmor radius is not taken into account), made possible by
the Julia language, serves as the foundation for the framework. The
OWCF provides answers to several questions, such as:(1)

where S is the diagnostic signal, W is the weight matrix modeling the phase-space sensitivity for each measurement bin and F isthe fast-ion distribution in phase space. As will be shown, equation(1) is one of the problems that our framework will help solve. The
sensitivity of the FI diagnostic can be expressed with varying dimensionality depending on the phase space of interest. For example, in velocity-space tomography [11–15], the sensitivity is twodimensional and can be expressed as a so-called weight functionof the parallel and perpendicular velocity components ( v || , v ⊥ ) of
the ion with respect to the magnetic ﬁeld (or, equivalently, via theenergy E and pitch p = v || / v [16]). In terms of computational resources, 2D weight functions are usually possible to compute ona regular laptop, for any phase-space grid resolutions of interest;in terms of both RAM and CPU power. For example, on an Intel(R)Core(TM) i7-8665U CPU @ 1.90 GHz laptop, single-thread computing 2D weight functions using the DRESS code [9] for the neutron emission spectroscopy diagnostic TOFOR [17] with assumed
80 measurement bins and a 400x200 (E,p) grid takes about 10
minutes and requires roughly 1 GB of RAM (including all required
input data). Assume we add another dimension that we discretize
into N grid points. Then the computation time would be roughly
N times as long. The RAM requirement would increase due to several factors, e.g. the need to hold more synthetic diagnostic signals
in memory as well as modeling the sightline of the fast-ion diagnostic in position space. This is because 2D weight functions canrelate a diagnostic measurement to the FI distribution only at asingle point ( R , z) in position space [18,19,15,20–22], where R is
the major radius cylindrical coordinate and z the vertical coordinate. To be able to relate a measurement to the FI distribution atall ( R , z) points via FIOs, three-dimensional so-called orbit weightfunctions (OWs) can be used [8,3,19,23]. The ( E , pm , R m ) so-called
orbit space is the 3D coordinate space of choice for the current
version of the OWCF, due to its semi-bounded space and favorability for tomographic reconstructions [3]. E is the energy of thefast ion and pm is the pitch at the maximum major radius position
R m of the fast ion as it traverses its orbit. All quantities computedwith the OWCF can be mapped to the standard ( E , μ, P φ ; σ ) coordinates (where μ is the magnetic moment, P φ is the toroidal
canonical angular momentum and σ is a binary coordinate). This
provides the user with the ﬂexibility to analyze and export computed results in either ( E , pm , R m ) or ( E , μ, P φ ; σ ) coordinates.
This is useful for e.g. stability analysis, which is often performed inthe ( E , μ, P φ ; σ ) coordinates. It should be mentioned that, in addition to OWs, other ways of relating a diagnostic signal to the full

• How to compute orbit weight functions eﬃciently,• How to analyze orbit weight functions effectively and intuitively,• How to streamline transformations of arbitrary fast-ion distributions into fast-ion orbits,• How to investigate populatable fast-ion orbits for every ( R , z)point,• How to compute synthetic diagnostic signals using orbit
weight functions with an exceptionally large amount of orbitspace grid points (> 100 × 100 × 100),• How to decompose diagnostic signals into the contributions
from fast-ion orbit types,• How to streamline the possibility of reconstructing fast-ion orbit distributions from measurements.
Before being structured into a shareable framework, the toolsof the OWCF have been used to develop and analyze orbit weightfunctions for neutron emission spectroscopy [8] and one-step
gamma-ray spectroscopy diagnostics [3]. The tools have also beenused to facilitate discussion of an optimal fast-ion diagnostics setup [27], and to analyze the most likely orbit-type constituents ofa decrease in neutron measurements during a discharge in the JET
deuterium-tritium campaign [28].
In addition, even though the OWCF is currently heavily integrated with the DRESS [9] code for computing synthetic diagnostic
neutron and γ -ray spectra, the OWCF maintains a modular approach towards synthetic diagnostics codes in general. This enables
future codes to be easily integrated with the OWCF, for computing
OWs via synthetic diagnostic FIO spectra.
This paper is structured as follows. In section 2, an overview of
the OWCF is given. In sections 3 and 4, discussions on topological maps and orbit weight functions are provided, respectively. In
2

H. Järleblad, L. Stagner, M. Salewski et al.

Computer Physics Communications 294 (2024) 108930

Fig. 2. A graphical overview of the OWCF. (For interpretation of the colors in the ﬁgure(s), the reader is referred to the web version of this article.)

5, the capability of the OWCF to transform between ( E , pm , R m )and ( E , p , R , z) is outlined and in section 6 the tools of the OWCF
for decomposing diagnostic signals in terms of orbit types are presented. Finally, a conclusion follows in section 7.
It should be mentioned that the OWCF contains more utilities
and tools than can be included in this paper. For a complete anddetailed breakdown of the whole OWCF, the reader is referred tothe manual included with the OWCF.

yet maintains high-performance via eﬃcient code compilation andexecution [26].
To help the reader understand the structure of the OWCF, a basic building block of the framework has been visualized in Fig. 1a,
and a basic relationship between different scripts has been visualized in Fig. 1b. We can observe how the inputs and outputs form
an I/O chain between the scripts, including optional inputs andintegrated usage of dependencies. Together, this simple example
serves as a basic illustration of the much larger system that is thefull OWCF.
A graphical overview of the full OWCF is given in Fig. 2. There
are four main groups of code (inside dotted lines): main scripts,
apps, data and misc. The main scripts (colored purple) constitutethe bulk of the computational tools of the OWCF. The apps pro-

2. Overview of the OWCFThe OWCF comprises many scripts and ﬁle structures, and it is
written almost entirely in the Julia programming language. It waschosen as a suitable language for the OWCF since it is high-level,3

H. Järleblad, L. Stagner, M. Salewski et al.

Computer Physics Communications 294 (2024) 108930

vide interactive analysis of computed quantities using web-basedplotting. We shall explore the user interface and capabilities of
the OWCF main scripts and apps in later sections. Furthermore,
as we can observe in Fig. 2, there are ﬁve forms of basic input
data for the OWCF: a magnetic equilibrium, a fast-ion (FI) distribution in ( E , p , R , z) coordinates, a diagnostic line-of-sight (LOS), a.cdf-ﬁle with TRANSP [30] discharge (shot) data and a .cdf-ﬁle with
TRANSP NUBEAM [31] FI data. TRANSP is a 1.5D equilibrium andtransport solver for interpretation and prediction of tokamak discharges. Its simulations and results are widely used in the plasma
physics community, thus making its output a natural basic data input for the OWCF. NUBEAM is the FI module of TRANSP.
Because of its versatile set of tools, the OWCF can be used in a
number of ways depending on the needs of the user. One integral
part of the framework that is likely to be used by almost all users
is the easy and eﬃcient computation of topological maps.

poloidal ﬂux is decreasing towards the edge, depending on deﬁnition) for a given R = R m . That iszm = maxz (ψ p ( R m , z))Topological maps enable the user to easily interpret and analyze computed quantities in 3D ( E , pm , R m ) orbit space [7,32,
23,19,8] and 4D ( E , p , R , z) guiding-center phase space. For each
( E , pm , R m ) or ( E , p , R , z) point, the corresponding guiding-center
drift orbit is computed by solving the equations-of-motion [33]p  B∗

γ m B ∗ṗ  = qE∗ ·+  E∗ ×

B∗

c b̂

(2)

B ∗

(3)

B ∗

where Ẋ is the time derivative of the guiding-center particle position, p  is the particle momentum parallel to the magnetic ﬁeld,
γ is the Lorentz factor, m is the particle mass,  is an ordering parameter, c is the speed of light, b̂ is the unit vector in the directionof the magnetic ﬁeld and q is the particle charge. B∗ and E∗ are
the effective magnetic and electric ﬁelds, respectively, deﬁned asB∗ ≡ B + E∗ ≡ E −


q

cp 
q

∇ × b̂p∂ b̂+ mc 2 ∇ γ∂t(6)where ψ p ( R , z) is the poloidal ﬂux as a function of R and z, zaxis isthe vertical coordinate value of the magnetic axis and dz is an arbitrary length within which the maximum ψ p value is to be foundwith certainty. dz naturally varies depending on machine.
In addition to the six basic orbit types mentioned in the previous section, the OWCF also deﬁnes ‘lost’, ‘incomplete’ and ‘invalid’
orbits. Lost orbits are particle trajectories that, when obtained via
integration of equations (2), (3), intersect the wall of the machine.
Incomplete orbits are orbits that the integration algorithm cannot
complete (for any reason e.g. numerical errors). Invalid orbits are( E , pm , R m ) coordinates that do not correspond to physically possible orbits. When examining the orbit-space topology, it makes
most sense to look at slices of constant energy E, and not pm or
R m . This is because only slices of constant E contain all orbit types[8,3]. An overview of how the OWCF computes topological maps in( E , pm , R m ) orbit space is shown in Fig. 3.
In addition to topological maps, the OWCF also allows theuser to compute maps of the poloidal and toroidal transit timein ( E , pm , R m ) orbit space. This can be useful for e.g. examining
magnetic equilibria in detail or resonances between particles andAlfvén eigenmodes [6,36].
The orbitsWebApp.jl (A4 in Fig. 2) can be used to interactively
visualize the ( E , pm , R m ) topological maps, poloidal and toroidaltransit times. A screenshot of the web application can be seen in
Fig. 4. We can see how the web application also visualizes the current pm , as well as a top-view and a poloidal projection of theorbit trajectory for the current ( E , pm , R m ) coordinate. In the plotof the toroidal transit time τt (middle-left plot), it is interesting
to note the narrow curved region of relatively large values of τt .
This corresponds to trapped orbits that, due to the speciﬁc magnetic equilibrium, precess toroidally approximately as much in the
co-current as in the counter-current direction. With the orbitsWebApp.jl web application, the user can change the values of theE, pm and R m sliders interactively and analyze ( E , pm , R m ) orbit
space in real-time. ( E , μ, P φ ; σ ) constants-of-motion space (where
μ is the magnetic moment, P φ is the toroidal canonical angularmomentum and σ is a binary coordinate) can be examined via a
toggle button. This switch has been illustrated in Fig. 5.

3. Topological mapsẊ =∀ z ∈ [zaxis − dz, zaxis + dz]

(4)


(5)

where B is the magnetic ﬁeld, E is the electric ﬁeld and ∂∂b̂t is
the time derivative of the unit vector in the direction fo themagnetic ﬁeld. The equations-of-motion are relativistic to enable
users of the OWCF to work with particles with speeds close toc, e.g. runaway electrons [34,35]. Once the equations-of-motion
have been integrated once poloidally, we can identify the FIO typefor our phase-space point of interest. By identifying which partsof phase space correspond to the six basic types of FIOs (copassing, counter-passing, trapped, stagnation, counter-stagnationand potato) [8,3,36], we can identify topological regions and use
them to gain further insight about the phase-space quantities.
Topological maps also serve as a ‘discrete prism’ when splitting
synthetic diagnostic signals into their orbit-type constituents, as
we shall discuss in later sections.

3.2. ( E , p , R , z) phase space
The OWCF can also compute maps for the orbit topology,
poloidal and toroidal transit time in ( E , p , R , z) phase space, via
the calcEpRzTopoMap.jl main script (M1 in Fig. 2). Topologicalmaps for ( E , p , R , z) phase space are useful when the user would
like to know which orbit types pass through a certain ( R , z) point,
what the ( E , p ) topological map looks like and how that changes
from one ( R , z) point to another.
Data computed in ( E , p , R , z) phase space can be visualizedusing the EpRzWebApp.jl (A3) web application, as illustrated inFig. 6. We can observe how, at the ( R , z) point of interest, the valid
orbits consist mostly (> 80%) of co- and counter-passing orbits. At
pitch values close to zero, the valid orbits are exclusively trapped
orbits for all energies.
The EpRzWebApp.jl web application is envisioned to assist fastion diagnostics that observe small (compared with minor radius)
measurement volumes, e.g. fast-ion D-α diagnostics, in determining which orbit types are able to produce signals at given ( R , z)points. Several other usage areas can also be imagined, such as forteaching and a deeper insight into the populated orbits for a fastion distribution given in ( E , p , R , z) coordinates. For example, it

3.1. ( E , pm , R m ) orbit space
The computation of topological maps in orbit space is done byintegrating equations (2), (3) for all of ( E , pm , R m ) orbit space until the particle has completed its orbit poloidally and then inferring
the orbit type. The vertical starting point of the integration, zm isfound by ﬁnding the minimum poloidal ﬂux (or maximum, if the4

H. Järleblad, L. Stagner, M. Salewski et al.

Computer Physics Communications 294 (2024) 108930

Fig. 3. A graphical overview of the process by which the OWCF computes topological maps for the ( E , pm , R m ) coordinate space. All the energy slices (E 1 , E 2 , ..., E n ) arecomputed sequentially and then put together into a 3D topological map of the whole ( E , pm , R m ) coordinate space.

can be imagined how students can use the OWCF apps to perform
numerical experiments by modifying the orbit parameters interactively. This could provide a quick and effective way of teaching FIOs
and orbit-space topologies, and help the students faster develop an
intuition for FIOs and FIO-related quantities.

where w ( E 1,d , E 2,d ) is the orbit weight function for the diagnostic measurement bin with lower and upper limits E 1,d and E 2,drespectively. τ p is the poloidal transit time of the orbit, S is the
expected measurement signal of the orbit, t is time, γ is the gyroangle and φ0 is the initial toroidal angle of the orbit. The orbit
weight function is the expected signal of an orbit averaged over all
reducible coordinates [23,19].
The calcOrbWeights.jl main script computes orbit weight functions using equation (7) while also taking the geometry of thediagnostic line-of-sight (LOS) into account. When computing OWs
for e.g. neutron emission spectroscopy [37] or gamma-ray spectroscopy [38,39] by coupling the OWCF to the DRESS code, there
are several fusion reactions for the user to choose between. The
following fusion reactions are currently supported by the OWCF
(via DRESS) [40,41]:

4. Orbit weight functions
Orbit weight functions (OWs) quantify how sensitive a diagnostic is to different fast-ion drift orbits. OWs are arguably one of themain features of the OWCF. Hence, the discussion concerning OWs
will be given extra care, and it has therefore been split into two
subsections: computation and visualization.
4.1. Computation
The computation of OWs is done via the calcOrbWeights.jl
main script (M2 in Fig. 2). As has been discussed in [8,3,23], theprocess of computing orbit weight functions in ( E , pm , R m ) orbit
space can be mathematically expressed as• D + D → n (2.45 MeV) + 3 He (0.82 MeV)
• D + T → n (14.1 MeV) + 4 He (3.5 MeV)
• D + 3 He → p (14.7 MeV) + 4 He (3.6 MeV)
• T + p → γ (19.8 MeV) + 4 He

w ( E 1,d , E 2,d )=1
4π 2 τ p

Please note however, that any other fusion reaction can be
added to the list, by simply implementing the required crosssectional data. Also, please note that the fusion-product proton
from the D(3 He,p)4 He reaction has a non-zero charge. The pro-2π2πτ p
S ( E 1,d , E 2,d , E , pm , R m , t , γ , φ0 )dtdγ dφ0 (7)
0

0

0
5

H. Järleblad, L. Stagner, M. Salewski et al.

Computer Physics Communications 294 (2024) 108930

Fig. 4. A screenshot of the orbitsWebApp.jl web application for interactively visualizing FIOs, topological maps, poloidal and toroidal transit time maps in ( E , pm , R m ) orbit
space and ( E , μ, P φ ; σ ) constants-of-motion space. ‘τ p (7.6, 0.26) = 13.5 microseconds’ means that τ p at ( R m , pm ) = (7.6 m, 0.26) is 13.5 microseconds.

ton trajectory would therefore have to be followed post-fusion to
enable diagnostic investigation. As of the current version of theOWCF, this is not implemented. Instead, the resulting proton energy distribution for the plasma as a whole is returned (assuming
4π emission, i.e. ignoring any diagnostic LOS input). In future versions of the OWCF, implementation could utilize e.g. the methods
developed in [42].
When computing OWs numerically, the temperature and density proﬁles of the thermal (slow) particle species in the fusion
reaction must be speciﬁed. With the OWCF, they are speciﬁed asfunctions of ρ pol ,ρ pol =

ψ p − ψ p ,axis
ψ p ,sep − ψ p ,axis

e.g. neutron or gamma-ray energies. The projected velocity of afast-ion can be written as [11]u = v || cos φ + v ⊥ sin φ cos ,where v || and v ⊥ are the fast-ion velocity components parallel andperpendicular to the magnetic ﬁeld, respectively. φ is the observation angle between the line-of-sight and the magnetic ﬁeld vector
and is the gyrophase [43]. The projected velocity is a good proxyfor spectral fast-ion diagnostics since it reﬂects essential featuresof the spectrum formation, due to the conservation of energy andmomentum in a fusion reaction [11]. u is an analytic measure of
the sensitivity of a diagnostic in the sense that u does not require
any thermal density and temperature data. When u has been computed, a value of w is added to its corresponding velocity bin. The
binning weights w are

1/2
,

(9)

(8)

where ψ p is the poloidal magnetic ﬂux function, ψ p ,axis is thevalue of ψ p at the magnetic axis and ψ p ,sep is the value of ψ p
at the plasma separatrix. If no proﬁles are speciﬁed, the OWCF willuse default (re-scalable) proﬁles. These are illustrated in Fig. 7.
Continuing, calcOrbWeights.jl also gives the user the option ofcomputing weights binned into projected velocities u instead ofw=τp  φ,τ p 2π(10)where τ p /τ p is the fraction of the total poloidal transit time for
the FIO spent at the point of interest,  is the solid angle of the6

H. Järleblad, L. Stagner, M. Salewski et al.

Computer Physics Communications 294 (2024) 108930

Fig. 5. An example of how the OWCF can be used to map a) topological regions and b) poloidal transit times (τ p ) from ( E , pm , R m ) orbit coordinates to ( E , μ, P φ ; σ )constants-of-motion coordinates. The slices in b) correspond to the same energy slices as in a). We can observe how the apparent size of the topological regions change
between the phase spaces. σ = −1 corresponds to counter-going (counter-current) orbits and σ = +1 corresponds to co-going (co-current) orbits. The lost region (brown)
was not included for the poloidal transit times.

diagnostic viewing cone voxel and φ/2π is the toroidal angle
fraction that the diagnostic viewing cone voxel of interest occupies. As mentioned earlier, the ﬁniteness of the Larmor radius is
not taken into account by the OWCF. The process of creating theresulting projected velocity spectrum has been illustrated in Fig. 8.
Let u i , j denote the projected velocity of sample j from gyro-centeri. In Fig. 8: A Two samples (yellow), with projected velocities u 1,1
and u 1,2 , along the gyro-orbit of the same guiding-center (red) end
up in the same u-bin, even though u 1,1 = u 1,2 . This is because of
the ﬁnite discretization of u, and u 1,1 ≈ u 1,2 so the samples end up
in the same bin. B A sample with u = v ,1 cos φ . Because of cos
in eq. (9), not many samples end up here. C Samples from bothgyro-centers can end up in this bin. From the ﬁrst point becauseof up-shift from v ,1 cos φ and from the second point because ofdown-shift from v ,2 cos φ , where v ,1 < v ,2 . D Again, because ofcos a lot of samples will end up here, from both gyro-centers.
E A sample with v ,2 cos φ . F Even though this sample is outside of the LOS, the OWCF still includes it in the binning, since its
guiding-center is inside of the LOS.

4.2. Visualization
The orbit weight functions can be visualized with the weightsWebApp.jl web application (A7 in Fig. 2). With the OWCF, the orbit
weight functions can be visualized both as a ( E d , E , pm , R m ) 4D
matrix or as an equivalent ( E d , E , μ, P φ ; σ ) 5D matrix. In Fig. 9,
we can see a screenshot of the weightsWebApp.jl (A7) web application. We can observe how the S and WF signals match, verifyingthat the orbit weight functions are an accurate model of the diagnostic [11]. We can also observe how the product of the weightfunction (top-right plot) and the fast-ion distribution (middle-left)
result in a detailed view of where the diagnostic signal is likely tohave originated in orbit space (middle-right).
Another way of thinking about OWs is: given an orbit with thecoordinate ( E , pm , R m ), what is the expected signal? What signal isthat particular orbit likely to produce? To answer these questions,the user can employ the weightWebApp.jl web application (A8 in
Fig. 2, also please note the singular form weight in the name). The
user can scan the valid orbits of orbit space and examine what7

H. Järleblad, L. Stagner, M. Salewski et al.

Computer Physics Communications 294 (2024) 108930

Fig. 6. A screenshot of the EpRzWebApp.jl web application for interactively visualizing FIOs and topological maps in ( E , p , R , z) orbit space. The app controls have been
omitted for brevity. The top right bar plot shows the populatable orbits at ( R , z) as fractions. The bottom right plot shows the evolution of pitch as the fast ion traverses the
FIO. JET discharge 96100 at 13 s.

Fig. 7. The default bulk temperature and density proﬁles used by the OWCF, in casethe user provided none. ρ pol is the normalized poloidal ﬂux coordinate as deﬁned
in (8). The proﬁles are based on the temperature and density proﬁles of JET shot No
96100 at 13 seconds.

their expected signals look like. In Fig. 10, we can observe how
the user can examine the expected signals for various orbits andorbit types by changing the ( E , pm , R m ) coordinate with the sliders. In Fig. 10, the trapped orbit in focus has both its ‘banana tips’
inside the vertical sightline of the diagnostic, resulting in a doublehumped expected signal, plotted in the top-left plot.
5. Transforming to orbit space
A guiding-center fast-ion distribution f can be equivalently parameterized in several different sets of coordinates. The fast-ion
data from TRANSP (D5 in Fig. 2) is given on an irregular grid spiraling outwards from the magnetic axis, as shown in Fig. 11. To
ultimately transform the fast-ion distribution to ( E , pm , R m ) orbit
space, the OWCF ﬁrst transforms the data onto a rectangular gridin ( E , p , R , z) coordinates by creating a Delaunay tessellation forthe spiral ( R , z) grid [44,19] and then utilize linear barycentric interpolation [45]. The barycentric weights used by the OWCF has
been chosen as suggested by [46]. That is,Fig. 8. An illustration of the computation of projected velocity u spectrum by the
OWCF. The orbit trajectory is depicted as a pink line.

μ( R , z, R j , z j ) =λ( R , z, R j , z j )
3

j =1

(11)

λ( R , z, R j , z j )withλ( R , z, R j , z j ) =81( R − R j )2 + ( z − z j )2

,

(12)

H. Järleblad, L. Stagner, M. Salewski et al.

Computer Physics Communications 294 (2024) 108930

Fig. 9. A screenshot of the weightsWebApp.jl web application for interactively visualizing orbit weight functions. The app controls have been omitted for brevity.

Fig. 10. A screenshot of the weightWebApp.jl web application for interactively visualizing orbit weight functions. The top left-hand plot shows the expected signal for theorbit visualized in the bottom plots. The app controls have been omitted for brevity.
9H. Järleblad, L. Stagner, M. Salewski et al.

Computer Physics Communications 294 (2024) 108930

=
h

pm

h

Rm

=

where ( R j , z j ) are the individual coordinates of the three verticesof the Delaunay triangle containing the query point ( R , z). Query
points outside of the Delaunay tessellation will be treated withnearest neighbor extrapolation, followed by forcing all f ( E , p ) outside of the separatrix to be identically zero.
With f ( E , p , R , z) obtained from the TRANSP data, the fast-ion
distribution in ( E , pm , R m ) orbit space can be obtained via MonteCarlo sampling. This is done as follows. The 4D ( E , p , R , z) array is‘spaghettiﬁed’ into a 1D vector. The cumulative sum vector is thencomputed and an evenly distributed random number between 0
and the last cumulative sum vector element is drawn. The corresponding element in the vector is identiﬁed, as well as the pertaining ( E , p , R , z) coordinate. The FIO is then computed and the( E , pm , R m ) coordinate of the FIO is identiﬁed. This is done for asuﬃcient number of Monte-Carlo samples and the f ( E , pm , R m )distribution is obtained. The f ( E , p , R , z) → f ( E , pm , R m ) process
has been visualized in Fig. 12. Other methods of transformingf ( E , p , R , z) → f ( E , pm , R m ) might be integrated in the OWCF in
updated versions of the framework. This is discussed in the future
work section of this paper.
6. Orbit constituents of diagnostic signals(13)

where ‘h’ stands for the different orbit types (co-passing, trappedetc). To investigate the dependency on ( E , pm , R m ) in terms of orbit types, we can further expand (13) asWh Fh =

W h, E F h, E
h(16)

The OWCF is a framework of scripts, functions, routines andapps that combine to enable computation, visualization andanalysis of fast-ion orbits and quantities in ( E , pm , R m ) (orbit space) ( E , μ, P φ ; σ ) (constants-of-motion) and ( E , p , R , z)
(guiding-center) coordinates, such as orbit weight functions andfast-ion distribution functions. It is written in the Julia programming language to enable eﬃcient computation of fast-ion orbits,
pushing the computation time down to ∼ 1 ms (deeply passing∼ 1 ms, marginally trapped ∼ 3 ms). With the OWCF, the user cancompute topological maps of the different orbit types and view
them in any of the supported coordinate spaces. The topologicalboundaries can then be superimposed onto e.g. orbit weight functions to analyze the sensitivity of fast-ion diagnostics to different
orbit types. The orbits and orbit-space quantities can be visualizedhh

W h, R m F h, R m

7. Conclusion

With the OWCF, the user can also decompose synthetic diagnostic signals into orbit type constituents, to investigate whatfractions of orbit types are most likely to have produced a certain
diagnostic measurement. Mathematically, this can be expressed as[3]

Wh Fh,(15)

where e.g. F h, E is the fast-ion distribution for a particular orbit
type h at a certain fast-ion energy E.
To compute these kinds of signal splits, the user can employ
the ps2WF.jl (M6 in Fig. 2) main script. In addition to being able
to compute orbit split signals, the ps2WF.jl tool can also compute
W F signals for ultra high-resolution grids in orbit space. That is,
grids that consist of tens of millions of valid orbits and that require
a computational cluster to compute within a reasonable timeframe(minutes instead of hours, assuming one orbit takes ∼ 1 ms tocompute). These kinds of grids can be useful when the S = W F
identity needs to be validated. S is the synthetic signal from e.g. a
regular Monte-Carlo code and the validity of the OWs depends on
the perfect match between the W F vector and the S vector. For
too coarse grids, the identity S = W F will likely not hold. This isbecause accurate representation of the orbit-space sensitivity (i.e.
W ) of the diagnostic requires a discretization of ( E , pm , R m ) orbit
space with a suﬃciently high resolution.
The resulting output of the ps2WF.jl script can be visualizedusing the signalWebApp.jl (A6 in Fig. 2) web application. This has
been illustrated in Fig. 13. We can observe how the OWCF predicts
that about 60% of the measurement signal in the 2500 keV measurement bin (indicated by a red dot in the top row left-hand plot)
is likely to have originated from co-passing orbits (top row righthand plot)and about 20% from stagnation orbits. The OWCF also
shows that stagnation orbits are more likely to have contributed
to the measurement signal than trapped orbits, for fast-ion energies above approximately 250 keV (second row plot). This isbecause the fast-ion distribution is likely to consist of more stagnation orbits than trapped orbits above that energy (third row plot),and because this particular diagnostic sightline (oblique sightline
in JET) is more sensitive to stagnation orbits than trapped orbits
above the same energy level, approximately (bottom row plot).
Furthermore, if a simple overview of the orbit-type constituentsof a diagnostic signal is sought, the OWCF provides quick and easyto-use tools in gui.jl (E3 in Fig. 2) that can decompose a synthetic
diagnostic signal into orbit-type origin, as shown in Fig. 14. We
can observe how the diagnostic signal is likely to have originatedmostly from co-passing orbits and about a third from trapped orbits. To know the likely orbit-type origin of diagnostic signals canbe useful when e.g. verifying certain heating schemes, such as thethree-ion heating scheme [47]. It can also be useful when investigating the likely orbit types of the birth distribution of alphaparticles in future burning plasmas.

Fig. 11. An example of a Delaunay tesselation performed on TRANSP fast-ion distribution data known only at discrete, irregular ( R , z) points. The energy E andpitch p dependency of the distribution was integrated out to enable illustrative
plotting, i.e. f ( R i , zi ) = f ( E , p , R i , zi )dEdp. f ( R i , zi ) was then normalized so thatmax( f ( R i , zi ))= 1.0 to enable clean plotting. The true maximum value of f ( R i , zi )is shown in the plot title as max( f ( R i , zi ))=1.0.

WF =W h, pm F h, pm

(14)

E
10

H. Järleblad, L. Stagner, M. Salewski et al.

Computer Physics Communications 294 (2024) 108930

Fig. 12. An example of f ( E , p , R , z) → f ( E , pm , R m ). In a) the f ( R , z) distribution has been illustrated, where the E- and p-dependences have been integrated out. In b) the
opposite is shown, i.e. f ( E , p ). In c) the resulting fast-ion distribution in orbit-space f ( E , pm , R m ) has been split into slices of constant E.

interactively via the OWCF web applications, where fast-ion orbits are computed in real-time based on the user’s controls, madepossible by highly optimized guiding-center codes [44,48].
The OWCF has been designed to be modular in terms of synthetic diagnostics codes. In addition to the already integrated
DRESS code, future codes such as FIDASIM can be added to the
OWCF framework by simply altering a few lines of code in the
OWCF.
Furthermore, the OWCF also contains tools for splitting synthetic diagnostic signals into their likely orbit-type constituents.
This can be analyzed at varying levels of detail, depending on theneeds of the user.
Tools that allow for transforming quantities between ( E , pm ,R m ), ( E , μ, P φ ; σ ) and ( E , p , R , z) are incorporated into the OWCF.
This enables the user to compute orbit-related quantities in one
coordinate space and analyze them in any of the other supported
spaces.
Future work for the OWCF includes continuing to make thescripts and apps more user-friendly by improving their eﬃciency
and speed, making the OWCF compatible with additional forward
models and computational cluster architectures, as well as adding
further transformation capabilities to be able to transform more
eﬃciently back and forth between ( E , pm , R m ), ( E , μ, P φ ; σ ) and( E , p , R , z). Such capabilities include the methods developed by[49] and [50], that transforms via grid-optimized automatic differentiation and marker distributions, respectively.
Finally, the versity and standard of the OWCF will provide a
useful and easy-to-use toolbox for understanding the complexityof fast-ion orbits in tokamaks. It will help pave the way for future
tomographic reconstructions of the fast-ion distribution using orbit

weight functions, and contribute to the realization of fusion as a
virtually limitless source of sustainable energy.
Declaration of competing interest
The authors declare the following ﬁnancial interests/personal
relationships which may be considered as potential competing interests:
The authors report ﬁnancial support was provided by Royal
Danish Academy of Sciences and Letters. The authors report ﬁnancial support was provided by European Consortium for the Development of Fusion Energy.
Data availability
Raw data were generated at the JET large scale facility. Derived
data supporting the discussions of this study are available from thecorresponding author upon reasonable request.
Acknowledgements
This work has received support from the Niels Bohr Foundation which is a merger of The Niels Bohr Grant, The Emil Herborg
Grantpart, The Grant of MA Marcus Lorenzen, The Ole Rømer Foundation and The Julie Marie Vinter Hansen Travel Grant. The views
and opinions expressed herein do not necessarily reﬂect those of
The Royal Danish Academy of Sciences and Letters.
This work has been carried out within the framework of theEUROfusion Consortium, funded by the European Union via the
Euratom Research and Training Programme (Grant Agreement No.
11

H. Järleblad, L. Stagner, M. Salewski et al.

Computer Physics Communications 294 (2024) 108930

Fig. 13. An example screenshot of the signalWebApp.jl (A6) web application. The user controls have been omitted for brevity. JET shot 94701 at 10.8 seconds, TRANSP ID V01.

Fig. 14. An example of a) a diagnostic signal W F split into its orbit-type constitutents. In b) the signal amplitude has been set to 1.0 for all diagnostic measurement bins, toallow the orbit-type fractions to be examined in detail. The absolute and normalized diagnostic signal has been superimposed for reference in a) and b), respectively.

101052200 — EUROfusion). Views and opinions expressed are

References

however those of the author(s) only and do not necessarily reﬂect[1] B. Madsen, Reconstructing the fast-ion velocity distribution in fusion plasmas
from sparse datasets, Ph.D. thesis, Technical University of Denmark, 2020.
[2] J. Mailloux, et al., Nucl. Fusion 62 (2022) 042026, https://doi.org/10.1088/17414326/ac47b4.

those of the European Union or the European Commission. Neither
the European Union nor the European Commission can be held responsible for them.
12

H. Järleblad, L. Stagner, M. Salewski et al.

Computer Physics Communications 294 (2024) 108930[3] H. Järleblad, et al., Nucl. Fusion 62 (2022) 112005, https://doi.org/10.1088/
1741-4326/ac63d3.
[4] S.D. Pinches, et al., Plasma Phys. Control. Fusion 46 (2004) B187, https://doi.
org/10.1088/0741-3335/46/12B/017.
[5] S.E. Sharapov, et al., Nucl. Fusion 53 (2013) 104022, https://doi.org/10.1088/0029-5515/53/10/104022.
[6] W.W. Heidbrink, Phys. Plasmas 15 (2008) 055501, https://doi.org/10.1063/1.
2838239.
[7] J.A. Rome, Y-K.M. Peng, Nucl. Fusion 19 (1979) 1193, https://doi.org/10.1088/0029-5515/19/9/003.
[8] H. Järleblad, et al., Rev. Sci. Instrum. 92 (2021) 043526, https://doi.org/10.1063/
5.0040696.
[9] J. Eriksson, et al., Comput. Phys. Commun. 1996 (2016) 40–46, https://doi.org/
10.1016/j.cpc.2015.10.010.
[10] M. Salewski, et al., Nucl. Fusion 56 (2016) 106024, https://doi.org/10.1088/0029-5515/56/10/106024.
[11] M. Salewski, Fast-ion diagnostic in fusion plasmas by velocity-space tomography, Dr. Techn. Thesis, January 2020.
[12] M. Salewski, et al., Nucl. Fusion 52 (2012) 103008, https://doi.org/10.1088/0029-5515/52/10/103008.
[13] M. Salewski, et al., Nucl. Fusion 57 (2017) 056001, https://doi.org/10.1088/
1741-4326/aa60e9.
[14] M. Salewski, et al., Fusion Sci. Technol. 74 (1-2) (2018) 23–36, https://doi.org/
10.1080/15361055.2017.1380482.
[15] M. Salewski, et al., Nucl. Fusion 58 (2018) 096019, https://doi.org/10.1088/
1741-4326/aace05.
[16] B. Madsen, et al., Nucl. Fusion 60 (2020) 066024, https://doi.org/10.1088/17414326/ab82b5.
[17] M. Gatu Johnson, et al., Nucl. Instrum. Methods Phys. Res., Sect. A 591 (2008)
417–430, https://doi.org/10.1016/j.nima.2008.03.010.
[18] W.W. Heidbrink, et al., Plasma Phys. Control. Fusion 49 (2007) 1457–1475,
https://doi.org/10.1088/0741-3335/49/9/008.
[19] L. Stagner, et al., Nucl. Fusion 62 (2022) 026033, https://doi.org/10.1088/17414326/ac3ed2.
[20] M. Salewski, et al., Nucl. Fusion 54 (2014) 023006, https://doi.org/10.1088/0029-5515/54/2/023005.
[21] M. Salewski, et al., Plasma Phys. Control. Fusion 56 (10) (2014) 105005, https://doi.org/10.1088/0741-3335/56/10/105005.
[22] M. Salewski, et al., Nucl. Fusion 53 (6) (2013) 063019, https://doi.org/10.1088/0029-5515/53/6/063019.
[23] L. Stagner, W.W. Heidbrink, Phys. Plasmas 24 (2017) 092505, https://doi.org/10.
1063/1.4990391.
[24] B. Madsen, et al., Plasma Phys. Control. Fusion 62 (2020) 115019, https://doi.
org/10.1088/1361-6587/abb79b.
[25] B.S. Schmidt, et al., Nucl. Fusion 63 (7) (2023) 076016, https://doi.org/10.1088/
1741-4326/acd6a6.
[26] J. Bezanson, et al., SIAM Rev. 59 (1) (2017) 65–98, https://doi.org/10.1137/
141000671.

[27] H. Järleblad, Orbit-space sensitivity of fast-ion diagnostics, Ph.D. thesis, Technical University of Denmark, Kgs. Lyngby, 2022.
[28] H. Järleblad, et al., Fast-ion orbit origin of neutron measurements in the JET DT
campaign, to be submitted to Nucl. Fusion.
[29] I.J. Grady Booch, James Rumbaugh, Uniﬁed Modeling Language User Guide, 2ndedition, Addison-Wesley, Boston, 2005.
[30] J. Breslau, M. Gorelenkova, F. Poli, J. Sachdev, A. Pankin, G. Perumpilly, Transp,[Computer Software], https://doi.org/10.11578/dc.20180627.4, jun 2018.
[31] A. Pankin, D. McCune, R. Andre, et al., Comput. Phys. Commun. 159 (3) (2004)
157–184, https://doi.org/10.1016/j.cpc.2003.11.002.
[32] L. Stagner, Inference of the fast-ion distribution function, Ph.D. thesis, University of California, Irvine, 2018.
[33] X. Tao, A.A. Chan, A.J. Brizard, Phys. Plasmas 14 (2007) 092107, https://doi.org/
10.1063/1.2773702.
[34] J. Connor, R. Hastie, Nucl. Fusion 15 (3) (1975) 415, https://doi.org/10.1088/0029-5515/15/3/007.
[35] M. Rosenbluth, S. Putvinski, Nucl. Fusion 37 (10) (1997) 1355, https://doi.org/
10.1088/0029-5515/37/10/I03.
[36] W.W. Heidbrink, R.B. White, Phys. Plasmas 27 (2020) 030901, https://doi.org/
10.1063/1.5136237.
[37] J. Eriksson, et al., Plasma Phys. Control. Fusion 61 (1) (2018) 014027, https://doi.org/10.1088/1361-6587/aad8a6.
[38] M. Nocente, et al., Rev. Sci. Instrum. 81 (2010) 10D321, https://doi.org/10.1063/
1.3501386.
[39] M. Nocente, et al., Nucl. Fusion 57 (7) (2017) 076016, https://doi.org/10.1088/
1741-4326/aa6f7d.
[40] M. Kikuchi, K. Lackner, M.Q. Tran, Fusion Physics, International Atomic Energy
Agency, Vienna, Austria, 2012.
[41] J.E. Perry, S.J. Bame Jr., Phys. Rev. 99 (5) (1955) 1386.
[42] W.W. Heidbrink, et al., Plasma Phys. Control. Fusion 63 (2021) 055008, https://doi.org/10.1088/1361-6587/abeda0.
[43] M. Salewski, et al., Nucl. Fusion 55 (9) (2015) 093029, https://doi.org/10.1088/0029-5515/55/9/093029.
[44] L. Stagner, Inference of the fast-ion distribution function, Ph.D. thesis, University of California, Irvine, 2018.
[45] J.P. Berrut, G. Klein, J. Comput. Appl. Math. 259 (2014) 95–107, https://doi.org/
10.1016/j.cam.2013.03.044.
[46] F.F. Little, et al., Surf. Comput. Aided Geom. Des. (1983) 99–107.
[47] Ye.O. Kazakov, et al., Phys. Plasmas 28 (2021) 020501, https://doi.org/10.1063/
5.0021818.
[48] L. Stagner, H. Järleblad, https://github.com/JuliaFusion/GuidingCenterOrbits.jl,
2022.
[49] S. Benjamin, Distribution transforms for orbit tomography, Master’s thesis, The
Australian National University, July 2021.
[50] A. Bierwage, et al., Comput. Phys. Commun. 275 (2022) 108305, https://doi.org/
10.1016/j.cpc.2022.108305.

13