28. Proton Imaging Unit

Figure 28.1: The Proton Imaging unit directory tree.

The Proton Imaging unit fires proton beams onto the simulation domain and records their deflection due to electric and magnetic fields on detector screens. The function ProtonImaging is the main driver that orchestrates the proton imaging execution. The deflections of the protons are calculated using the Lorentz force. For each cell in the domain the average electric and magnetic fields are used and the electric and magnetic components do not change within each cell. The code allows for setting up several proton beams with possibly different activation times and several detector screens. Each beam is associated with only one detector, but a detector is allowed to record the protons coming from more than one beam. This way one can simulate multienergetic proton beams by splitting them into several computational beams with identical spatial (isospatial) specifications and identical detector target but using different proton energies. Stacked film detectors for resolving different energy components can be simulated by associating separate isospatial detectors to each energy component of multienergetic proton beams. The mapping 'beam $\rightarrow$ many detectors' is currently not allowed.

Tracing of the protons through the domain can be done in two ways: 1) crossing the domain during one time step or 2) allowing for time resolved proton tracing, i.e. tracing the protons through the domain during many time steps, if the velocities of the protons are sufficiently small. Although the proton imaging driver can be called during an entire simulation with proton beams being activated at certain simulation times, it is highly recommended to apply the proton imaging diagnostic as a post processing application on existing checkpoint files. Screen protons can be appended to existing detector files, but the proton imaging code does not write currently any data to the checkpoint files. Restarting a FLASH application from a previous proton imaging run is hence not recommended, especially if the time resolved proton imaging version is used. Old disk protons waiting to be transported though the domain may be out of sync when restarting a previous run.

28.0.1 Proton Deflection by Lorentz Force

In the presence of an electric and magnetic field, the (relativistic and non-relativistic) Lorentz force in CGS units on a proton in motion is:

$${\bf F} = Q({\bf E} + {\bf v}/c\times {\bf B}),$$ (28.1)

where $Q$ is the charge of the proton in esu (g$^{1/2}$ cm$^{3/2}$ s$^{-1}$), ${\bf F}$ the force in dyne (g cm s$^{-2}$), ${\bf v}$ the proton velocity (cm s$^{-1}$), $c$ the speed of light ( $2.998\times 10^{10}$ cm s$^{-1}$) and ${\bf E}$ and ${\bf B}$ are the electric and magnetic fields in Gauss (g$^{1/2}$ cm$^{-1/2}$ s$^{-1}$). The non-relativistic acceleration the proton experiences is:

$${\bf a} = Q_m({\bf E} + {\bf v}/c\times {\bf B}),$$ (28.2)

where $Q_m$ is the mass rescaled charge $Q/m$ of the proton. For uniform fields in each cell (constant $(E_x,E_y,E_z)$ and $(B_x,B_y,B_z)$), this leads to a set of 1st order coupled differential equations in time

$${\dot v}_x(t) = Q_m[E_x + v_y(t)B_z/c - v_z(t)B_y/c]$$ (28.3)

$${\dot v}_y(t) = Q_m[E_y + v_z(t)B_x/c - v_x(t)B_z/c]$$ (28.4)

$${\dot v}_z(t) = Q_m[E_z + v_x(t)B_y/c - v_y(t)B_x/c]$$ (28.5)

for the velocity components $v_x(t)$, $v_y(t)$ and $v_z(t)$. Given the initial proton velocity components $(v_{0x},v_{0y},v_{0z})$ at time $t=0$, the analytical solution to this system of coupled differential equations becomes (assuming $B\neq 0$):

$$v_x(t) = (b_x^2E_x+b_xb_yE_y+b_xb_zE_z)Q_mt$$

$$+ (b_x^2 + b_y^2\cos[BQ_mt/c] + b_z^2\cos[BQ_mt/c])v_{0x}$$

$$+ b_ze_y(1 - \cos[BQ_mt/c])c + b_ye_z(\cos[BQ_mt/c] - 1)c$$

$$+ b_xb_y(1 - \cos[BQ_mt/c])v_{0y} + b_xb_z(1 - \cos[BQ_mt/c])v_{0z}$$

$$+ (b_y^2e_x+b_z^2e_x-b_xb_ye_y-b_xb_ze_z)\sin[BQ_mt/c]c$$

$$+ (b_x^2b_z+b_y^2b_z+b_z^3)\sin[BQ_mt/c]v_{0y}$$

$$- (b_x^2b_y+b_z^2b_y+b_y^3)\sin[BQ_mt/c]v_{0z}$$ (28.6)

$$v_y(t) = [x\rightarrow y, y\rightarrow z, z\rightarrow x]\;v_x(t)$$ (28.7)

$$v_z(t) = [x\rightarrow z, y\rightarrow x, z\rightarrow y]\;v_x(t)$$ (28.8)

where $B$ is the magnitude of the magnetic field vector and $b_x=B_x/B$ and $e_x=E_x/B$ are the magnetic field magnitude rescaled magnetic and electric field components. The operator $[x\rightarrow y, y\rightarrow z, z\rightarrow x]$ stands for 'replace x with y, y with z and z with x' in the formula for $v_x(t)$. In case $B=0$, these equations reduce to: $v_x(t)=Q_mE_xt$, $v_y(t)=Q_mE_yt$ and $v_z(t)=Q_mE_zt$. Note that the units of $E_xQ_mt$ are the same as velocity. The quantities $b_x$, $e_x$ and $BQ_mt/c$ are dimensionless. Integrating $v_x(t)$, $v_y(t)$ and $v_z(t)$ over $t$, we get expressions for the positions $r_x(t)$, $r_y(t)$ and $r_z(t)$ as a function of time. The resulting equations have terms involving $t$, $t^2$, $\cos[BQ_mt/c]$ and $\sin[BQ_mt/c]$, and cannot be solved analytically. We therefore resort to a Runge-Kutta integration approach with the 6-dimensional ODE vector:

$$\frac{d}{dt} \left(\begin{array}{c}{\bf r} {\bf v}\end{array}\right) = \left(\begin{array}{c}{\bf v} {\bf a}\end{array}\right) = \left... ...in{array}{c}{\bf v} Q_m({\bf E} + {\bf v}/c\times {\bf B})\end{array}\right),$$ (28.9)

For certain orientations of ${\bf E}$ and ${\bf B}$, however, analytical solutions for ${\bf r}$ and ${\bf v}$ are possible (see the unit test section).

28.0.1.1 Relativistic Proton Equation of Motion

We state here briefly the relativistic proton equation of motion. Starting from the Lorentz equation 28.1, the force ${\bf F}$ is the time derivative of the relativistic momentum ${\bf p}$:

$${\bf F} = \frac{d{\bf p}}{dt} = \frac{d[\gamma({\bf v}) m{\bf v}]}{dt},$$ (28.10)

where $m$ is the rest mass of the proton and $\gamma({\bf v})$ is the Lorentz factor:

$$\gamma({\bf v}) = 1 / \sqrt{1-({\bf v}\cdot {\bf v}/c)^2} = 1 / \sqrt{1-(v/c)^2}.$$ (28.11)

Performing the differentiation on the r.h.s. of Eq.(28.10) using the chain rule (since both ${\bf v}$ and thus $\gamma$ depend on time), we arrive at:

$${\bf F} = \gamma({\bf v})m{\bf a}_{\perp} + \gamma({\bf v})^3 m{\bf a}_{\parallel},$$ (28.12)

where ${\bf a}_{\parallel}$ and ${\bf a}_{\perp}$ are the parallel and perpendicular components of the acceleration ${\bf a}=d{\bf v}/dt$ with respect to the velocity vector:

$${\bf a} = {\bf a}_{\parallel} + {\bf a}_{\perp} \;\;\;,\;\;\; {\b... ...erp} = 0 \;\;\;,\;\;\; {\bf v}\cdot {\bf a} = {\bf v}\cdot {\bf a}_{\parallel}.$$ (28.13)

Inverting Eq.(28.12) to find the acceleration from the force on a moving proton leads to:

$${\bf a} = \frac{1}{\gamma({\bf v})m}\left({\bf F} - \frac{[{\bf v}\cdot {\bf F}]{\bf v}}{c^2}\right)$$ (28.14)

with magnitude:

$$a = \frac{1}{\gamma({\bf v})m}\sqrt{F^2-\frac{1+\gamma({\bf v})^2}{c^2\gamma({\bf v})^2} [{\bf v}\cdot {\bf F}]^2},$$ (28.15)

and inserting the Lorentz force equation 28.1 we get the relativistic expression for the acceleration:

$${\bf a} = Q_{m\gamma}\left({\bf E} + {\bf v}/c\times {\bf B}-[{\bf v}/c\cdot {\bf E}]\;{\bf v}/c\right),$$ (28.16)

where $Q_{m\gamma}$ is the relativistic mass rescaled charge $Q/m\gamma({\bf v})$ of the proton. Comparing this expression with its non-relativistic version from Eq.(28.2), the main difference is not only the relativistic increase in mass ( $m\gamma({\bf v})$) but also an extra term of order $c^{-2}$ involving the electric field. The relativistic analytical velocity solutions corresponding to Eqs.(28.6-28.8) are much more complicated and involve extra logarithmic velocity component terms due to the extra relativistic velocity term.

28.0.1.2 Approximate Solutions to the Non-Relativistic Proton Equation of Motion

Since Runge-Kutta integration with constraints (due to the cell boundaries) is expensive, we wish to derive conditions under which the proton path through a cell approximates a parabola, for which quadratic solutions to the path positions become available. To this end we expand ${\bf v}(t)$ in terms of $t$ around the cell entry point $t=0$. We get:

$${\bf v}(t) = {\bf v}_0 + \left\vert\frac{d{\bf v}(t)}{dt}\right\vert _{t=0}t + \frac{1}{2}\left\vert\frac{d^2{\bf v}(t)}{dt^2}\right\vert _{t=0}t^2 + O[t^3]$$ (28.17)

$$= {\bf v}_0 + {\bf a}_0t + \frac{{\bf j}_0}{2}t^2 + O[t^3],$$ (28.18)

where

$${\bf a}_0 = Q_m({\bf E} + {\bf v}_0/c\times {\bf B})$$ (28.19)

$${\bf j}_0 = Q_m({\bf a}_0/c\times {\bf B})$$ (28.20)

are the acceleration and jerk vectors at the cell entry point. Integration of Eq.(28.18) yields the path equation:

$${\bf r}(t) = {\bf r}_0 + {\bf v}_0t + \frac{{\bf a}_0}{2}t^2 + \frac{{\bf j}_0}{6}t^3 + O[t^4].$$ (28.21)

The leading parabolic error term is hence the term involving the jerk vector. The largest parabolic error term $\epsilon$ can be estimated by:

$$\epsilon = \frac{\max {\bf j}_0}{6}t^3.$$ (28.22)

If $\epsilon$ is less than the allowed positional error of the computation, the parabolic approximation is used. To get also the velocity components with the same accuracy, the jerk terms are used to compute the velocities. If the jerk terms are omitted from the velocity calculations, the velocity components would only be accurate to first order in $t$ and considerable error in velocity directions can accumulate during a proton imaging application.

In order to get a feeling for the conditions under which the parabolic approximation is applied in FLASH, let us first state the expression for the jerk vector ${\bf j}_0$, which is obtained by noting that ${\bf j}_0=d{\bf a}_0/dt$ and using the acceleration vector expression twice:

$${\bf j}_0 = (Q_m^2/c){\bf E}\times {\bf B} +(Q_m^2/c^2)({\bf v}_0\times {\bf B})\times {\bf B}.$$ (28.23)

The maximum magnitude that the ${\bf j}_0$ vector can achieve is given by the situation in which the pairs ${\bf B},{\bf v}_0$ and ${\bf E},{\bf B}$ are perpendicular and ${\bf E}\times {\bf B}$ is opposite to ${\bf v}_0$. This leads to

$$\vert{\bf j}_0\vert \leq (Q_m^2/c)EB+(Q_m^2/c^2)B^2v_0$$ (28.24)

and the inequality will also hold for the individual components of ${\bf j}_0$. Assuming no electric field, we can estimate the magnitude of $B$ for which the parabolic approximation should be valid. In FLASH, the positional error is defined as a fraction $f$ times the minimum cell side size. For a cubic cell with sides $\Delta$ we can set up the parabolic condition on $B$ as follows:

$$\frac{Q_m^2B^2v_0}{6c^2}t^3 \leq \Delta f.$$ (28.25)

The time it takes to cross the cubic cell is of the order of $t=\Delta/v_0$ and we conservatively extend this time by the factor of $\sqrt[3]{6}$ to get rid of the 6 in the denominator. The parabolic condition on $B$ becomes

$$B \leq \frac{cv_0}{Q_m\Delta}\sqrt{f}.$$ (28.26)

For typical FLASH runs, we have for 20MeV protons $v\approx c/5$ and $\Delta$ is of the order of microns ($10^{-4}$cm) for typical simulations. The accuracy fraction $f$ is typically $10^{-6}$. We also have $Q_m\approx 2.9\times 10^{14}$ and $c=2.99\times 10^{10}$. Hence under these conditions we have:

$$B \leq \approx 10^7 ,$$ (28.27)

i.e., all magnetic fields under $10^7$ Gauss in a cell can be treated parabolically.

28.0.1.3 Approximate Solutions to the Relativistic Proton Equation of Motion

This follows along the lines of the previous section (28.0.1.2). Eq.(28.19) must be replaced by the relativistic expression from Eq.(28.16), replacing ${\bf v}$ by ${\bf v}_0$, and Eq.(28.20) must be replaced by the time derivative of Eq.(28.16), replacing in the final expression ${\bf v}$ by ${\bf v}_0$ and ${\bf a}$ by ${\bf a}_0$:

$${\bf a}_0 = Q_{m\gamma}\left({\bf E} + {\bf v}_0/c\times {\bf B}-[{\bf v}_0/c\cdot {\bf E}]\;{\bf v}_0/c\right),$$ (28.28)

$${\bf j}_0 = Q_{m\gamma}\left({\bf a}_0/c\times {\bf B} -[{\bf a}_0/c\cdot {\bf E}]\;{\bf v}_0/c -2[{\bf v}_0/c\cdot {\bf E}]\;{\bf a}_0/c\right),$$ (28.29)

where $Q_{m\gamma}$ is equal to $Q/m\gamma({\bf v}_0)$. Note the extra terms of order $c^{-2}$ when compared with Eqs.(28.19) and (28.20). The largest parabolic error expression in Eq.(28.22) remains the same. Assuming again no electric field, the same steps leading to the non-relativistic parabolic $B$ condition in Eq.(28.26), leads to:

$$B \leq \frac{cv_0}{Q_{m\gamma}\Delta}\sqrt{f},$$ (28.30)

the only difference being in using the relativistic mass rescaled proton charge. For the 20MeV protons $v\approx c/5$ this would lead to a factor of $\approx 1.02$ larger $B$ than in Eq.(28.27), due to the relativistically reduced magnitude of ${\bf j}_0$.

28.0.2 Setting up the Proton Beam

The setup of each proton beam follows closely the setup of the laser beams in 18.4.6. The major difference is that the proton imaging beams originate from 3D capsules instead of 2D planar crossection areas. This allows for simulating proton capsule aberrations in the proton beams. The radius of the capsule as well as its center location is specified for each beam at runtime, as well as each beams full conical aperture angle. Target coordinates are only needed for directional purposes and the target area is circular and perpendicular to the beam's direction. Additional runtime parameters needed for each beam are: 1) the launching time (the beam fires once the simulation time exceeds the launching time), 2) its target detector screen and 3) the number of protons to be launched and their energies. Once each beam has been set up, it is checked, if its capsule volume is completely outside of the computational domain. A beam capsule (partially) inside the domain is not allowed. A local set of 3D rectangular grid unit vectors ${\bf u}_1$,${\bf u}_2$,${\bf u}_3$ is constructed for each beam, such that ${\bf u}_3$ points from the capsule center to the target center and ${\bf u}_1\times {\bf u}_2={\bf u}_3$. The 3D grid unit vectors serve for generating statistical 3D points inside the spherical capsule as well as statistical 2D points inside the circular target area.

28.0.3 Creating the Protons

The protons inside each beam are created by connecting one-to-one the statistical 3D grid points inside the capsule with the statistical 2D grid points inside the circular target area. The proton directions inside each beam are therefore not conically collinear and are allowed to cross. The proton initial positions and velocities on the domain surface are calculated using the same strategy as presented in 18.4.7.6. Protons missing the domain can either be simply ignored or directly recorded on the detector screen.

28.0.3.1 Moving Protons through the Domain

As the total number of protons generated can be very large, special techniques have to be used to overcome the computational memory constraints. The Proton Imaging unit has been designed to work with a proton batch and screen proton bucket combination. Both the batch and the bucket are of much smaller size than the total protons and screen protons generated and can be seen as temporary storage devices for both entities. During a time step, the active proton beams start filling the generated protons into the beam proton batch. Once filled, the beam protons in the batch are sent through the domain, where they are either collected as disk protons into the disk proton batch (if time resolved proton imaging was requested) or as screen protons into the screen proton bucket. If the screen proton bucket is full, its contents (screen protons) are written to the corresponding detector files and emptied. The refilling of both the beam/disk proton batch and the screen bucket are independent of each other and both batch and bucket can be of different sizes. The main function ProtonImaging processes beam/disk proton batch after batch until no more active beam/disk protons are present. As an example, if at a time step we have 3 active beams with $1,000,000$ protons each and the dimension of the proton batch is set to $100,000$, a total of 30 beam proton batches will be generated and sent through the domain.

28.0.3.2 Proton Specifications

Each beam/disk proton needs to know its position $x,y,z$ and velocity $v_x,v_y,v_z$ as it traverses the domain. In addition to these 6 components it needs to know: 1) the time spent travelling in domain, 2) the current block and processor numbers, 3) a global identification tag, 4) the originating beam and target detector numbers. Screen protons on the other hand are only defined on the detector screen and hence need no specific information for domain travelling. Their information is reduced to a screen position coordinate pair $x,y$ and a detector number to which detector they are associated with.

Additionally, upon request, the proton imaging code is able to calculate extra magnetic diagnostic quantities ${\bf k}$ and $J$ for each beam/disk and screen proton:

$${\bf k} = \int {\bf B}\;d\ell,$$ (28.31)

$$J = \int {\bf v}_n \cdot (\nabla\times {\bf B})\;d\ell,$$ (28.32)

where $\int d\ell$ denotes path integration, ${\bf B}$ is the magnetic field vector and ${\bf v}_n$ is the unit normal velocity vector along the proton path. Transforming the path integrals into time integrals using the relation $d\ell=\vert{\bf v}(t)\vert dt$, we get the differentials:

$$\frac{d{\bf k}}{dt} = {\bf B}\vert{\bf v}(t)\vert,$$ (28.33)

$$\frac{dJ}{dt} = {\bf v}(t)\cdot (\nabla\times {\bf B}).$$ (28.34)

leading to extra 4 entries in the Runge-Kutta vector in Eq.(28.9):

$$\frac{d}{dt} \left(\begin{array}{c}{\bf r} {\bf v} {\bf k} J\end{array}\right) = \left(\begin{array}{c}{\bf v} Q_m({\bf E} + {\bf v}/c\times {\bf B}) {\bf B}v {\bf v}\cdot (\nabla\times {\bf B})\end{array}\right),$$ (28.35)

If the parabolic path approximation is used for a cell, then ${\bf v}(t)={\bf v}_0+{\bf a}_0t$, where ${\bf v}_0$ and ${\bf a}_0$ are the velocity and acceleration at cell entry. For a $J$ cell contribution we have (note this is a general result not depending on the parabolic approximation):

$$J = \int^T_0 {\bf v}(t)\cdot (\nabla\times {\bf B}) dt = \int^T_0 d{\... ...\times {\bf B}) = \left[{\bf r}(T)-{\bf r}_0\right]\cdot (\nabla\times {\bf B})$$ (28.36)

where $T$ is the cell crossing time and ${\bf r}_0$, ${\bf r}(T)$ are the cell entry and exit locations of the proton. For the parabolic ${\bf k}$ cell contribution we have;

$${\bf k} = {\bf B}\int^T_0 \vert{\bf v}_0+{\bf a}_0t\vert dt.$$ (28.37)

The parabolic path length integral can be solved exactly:

$$\int^T_0 \vert{\bf v}_0+{\bf a}_0t\vert dt = \frac{1}{2}\left[ \vert{\bf a}_0T+{\bf v}_0\vert T + \frac{{\bf a... ...0\vert^2}\left(\vert{\bf a}_0T+{\bf v}_0\vert-\vert{\bf v}_0\vert\right)\right.$$

$$+ \left.\frac{\vert{\bf a}_0\times {\bf v}_0\vert^2}{\vert{\bf a}... ...\bf a}_0\cdot {\bf v}_0+\vert{\bf a}_0\vert\vert{\bf v}_0\vert}\right)}\right].$$ (28.38)

The second term is the offperpendicular and the third term the offcolinear contribution to the path length. Taylor series expansion of this expression around ${\bf v}_0$ in terms of the components of ${\bf a}_0$ gives the zeroth order term as $\vert{\bf v}_0\vert T$ corresponding to the path length for ${\bf a}_0={\bf0}$. Likewise the zeroth order term for the expansion around ${\bf a}_0$ in terms of the components of ${\bf v}_0$ is $\vert{\bf a}_0\vert T^2/2$ which is the path length for ${\bf v}_0={\bf0}$. Application of the above general formula must be done with some care to avoid computational pitfalls. The cases where ${\bf v}_0$ and ${\bf a}_0$ are colinear must be handled separately. It is also easy to show, for example, that for the offcolinear cases the natural logarithm argument is always $\geq 1$.

28.0.4 Setting up the Detector Screens

Each detector screen is defined as a square planar area with a specific side length $s$ and an optional circular pinhole. The orientation of the screen in 3D space is specified by three components: 1) the center location of the screen ${\bf C}$, 2) the normal vector of the screen plane ${\bf n}$ and 3) the orthogonal unit vector pair ${\bf u}_x$,${\bf u}_y$ with origin at the center and oriented in such a way that each unit vector is parallel to two opposite sides of the detector screen. Location of the pinhole ${\bf H}$ is fixed by giving the distance between the two pinhole and detector centers and is always located on the detector side opposite to ${\bf n}$. The circular pinhole area is coplanar to the detector screen area. While ${\bf C}$ and ${\bf n}$ are sufficient to characterize the screen plane, the exact location of the detector square area is only specified once the orientation of ${\bf u}_y$ (or ${\bf u}_x$) is fixed. For that purpose we define a tilting angle $\alpha$ of ${\bf u}_y$ wrt one of the global 3D axes along the normal vector ${\bf n}$. A tilting angle of $\alpha=0$ would thus mean that the chosen global axis and ${\bf u}_y$ are coplanar.

28.0.4.1 Recording Protons on the Detector Screen

After each proton leaves the domain, it is located on the domain surface at a position ${\bf P}$ and with (outward from domain) velocity ${\bf v}$. The goal is to see, if the proton actually hits the screen plane and to record the local coordinates $(x,y)$ on the screen plane with coordinate basis $[{\bf u}_x,{\bf u}_y]$. After some algebra we arrive at the following conditions

$${\bf v}\cdot {\bf n} =$$ 0   proton moves parallel to screen plane (28.39)

$$\frac{({\bf C}-{\bf P})\cdot {\bf n}}{{\bf v}\cdot {\bf n}} >$$ 0   proton hits screen plane (28.40)

and the following expressions for the local coordinates

$$x = \frac{{\bf u}_y\cdot [{\bf v}\times ({\bf C}-{\bf P})]}{{\bf v}\cdot {\bf n}}$$

$$y = \frac{{\bf u}_x\cdot [({\bf C}-{\bf P})\times {\bf v}]}{{\bf v}\cdot {\bf n}}.$$ (28.41)

For better handling of the screen output data, the local coordinates are shifted and rescaled such that $(x,y)\in [0,1]$ when the proton is on the detector screen. Using the detectors side length $s$, we achieve this by the following local coordinate transformation:

$$x \longrightarrow (x + s/2)/s$$ (28.42)

$$y \longrightarrow (y + s/2)/s.$$ (28.43)

If a pinhole is present, the code has to check, if the proton makes it through the hole. Since the detector screen area and the circular pinhole area are coplaner, we can use the local coordinate basis $[{\bf u}_x,{\bf u}_y]$ located at the pinhole center ${\bf H}$ to calculate a local pinhole coordinate for each proton, in analogy to Eq.(28.41):

$$h_x = \frac{{\bf u}_y\cdot [{\bf v}\times ({\bf H}-{\bf P})]}{{\bf v}\cdot {\bf n}}$$ (28.44)

$$h_y = \frac{{\bf u}_x\cdot [({\bf H}-{\bf P})\times {\bf v}]}{{\bf v}\cdot {\bf n}}.$$ (28.45)

The proton makes it through the pinhole, if $h_x^2+h_y^2\leq r_h^2$, where $r_h$ is the radius of the pinhole.

The user is able to specify via a runtime parameter, if he wants only the detector screen protons $(x,y)\in [0,1]$ to be written to the detector file or if he also wants to record the offscreen protons $(x,y)\notin [0,1]$. Each detector file is currently an ascii file (formatted) and named:

$$\mbox{$<$basename$>$ProtonDetectorFile$<$detectorID$>$\_$<$timestamp$>$}$$ (28.46)

where $<$basename$>$ is the name of the simulation, $<$detectorID$>$ is the detector number (currently limited to the range [01-99]) and (optionally) $<$timestamp$>$ is the simulation time when the file was created. Only the $(x,y)$ pairs are written (first $x$ followed by $y$ on each line) without any headers. Thus the file can be directly incorporated and read by other software for graphical output (see for example the proton imaging ASCII $->$ PGM greyscale converter). Additionally, if requested, the detector file will contain in columns 3 to 6 the proton magnetic diagnostic quantities $J,k_x,k_y,k_z$ from Eqs.(28.31) and (28.32), in that order.

28.0.5 Time Resolved Proton Imaging

For some applications (slow moving protons, fast domain changing) it is desireable to adjust the domain environment as the protons move through it. The proton imaging unit has been given this capability by monitoring the time each proton spends in the domain and compare it to the computational time step. When the proton time exceeds the time step value, it is stored as a disk proton which will be written to disk. All disk protons from a previous time step must be first processed during a new time step before any new beam protons are generated. During a time step, the main driver ProtonImaging calls the following two tasks, each of which does the following operations:

Transport disk protons:

• Read all disk protons from old disk proton file
• Trace all disk protons through domain:
  • To new disk proton file (if proton time $>$ time step)
  • To screen detector
• Rename new disk proton file $\longrightarrow$ old disk proton file

Transport beam protons:

• Create all beam protons
• Trace all beam protons through domain:
  • To new disk proton file (if proton time $>$ time step)
  • To screen detector
• Append new disk proton file $\longrightarrow$ old disk proton file

The time resolved application requires extra memory and disk storage due to the presence of disk protons as well as extra reading/writing time from/to disk. If only a non-time resolved proton imaging run is requested, the code skips all array allocations and operations associated with the disk protons, and only the underlined actions of the transport beam protons task get activated.

28.0.6 Usage

To include the use of the Proton Imaging unit, the following should be included into the setup line command:  

+protonImaging [pi_maxBeams=<number> pi_maxDetectors=<number> threadProtonTrace=True]

The +protonImaging shortcut handles all the logistics for properly including the ProtonImaging unit. The default settings are: maximum number of beams = 1, maximum number of detectors = 1 and no use of threading during tracing of protons through the domain, which can be changed by the following three setup variables:

• pi_maxBeams: The maximum number proton beams for the simulation.
• pi_maxDetectors: The maximum number of proton imaging detectors for the simulation.
• threadProtonTrace: Enables threading during proton domain tracing.

Using proton imaging runtime parameters in the flash.par file, the user can set up the proper proton imaging configuration for his specific needs. The following list of runtime parameters is currently available for the user.

28.0.6.1 Proton Imaging Beams Runtime Parameters

The following are the runtime parameters for the proton imaging beams. The _n at the end of each runtime parameter characterizes the beam number, hence replace _n by _1 for the first beam, _n by _2 for the second beam, etc.

• pi_numberOfBeams: The number of proton imaging beams that are going to be used.
• pi_beamCapsule[X,Y,Z]_n: The global 3D components of the capsule center position vector ${\bf C}$.
• pi_beamCapsuleRadius_n: The radius of the spherical capsule.
• pi_beamTarget[X,Y,Z]_n: The global 3D components of the target position vector ${\bf T}$ for directional purpose.
• pi_beamApertureAngle_n: The conical full aperture angle.
• pi_beamProtonEnergy_n: The energy (in MeV) of the protons in the beam.
• pi_beamTime2Launch_n: The time during the simulation, when the beam should be fired.
• pi_beamDetector_n: The detector number, where the protons of the beam should be recorded.
• pi_beamNumberOfProtons_n: Number of protons to be launched from the beam.
• pi_beamNoBoundaryCondition_n: Ignore domain boundary conditions (reflective) when protons enter the domain?

28.0.6.2 Proton Imaging Detectors Runtime Parameters

The following are the runtime parameters for the proton imaging detectors. As for the beams, the _n at the end of each runtime parameter characterizes the detector number.

• pi_numberOfDetectors: The number of proton imaging detectors that are going to be used.
• pi_detectorCenter[X,Y,Z]_n: The global 3D components of the detector center position vector ${\bf C}$.
• pi_detectorNormal[X,Y,Z]_n: The local 3D components of the detector normal vector ${\bf n}$.
• pi_detectorSideLength_n: The side length (cm) of the detector square screen.
• pi_detectorSideTiltingAngle_n: The tilting angle (degrees) of two sides parallel to the detector screen unit axis ${\bf u}_y]$ with respect to the tilting axis.
• pi_detectorSideTiltingAxis_n: The global tilting axis ('x','y' or 'z') for the detector screen.
• pi_detectorAlignWRTbeamNr_n: A useful shortcut to place the detector screen right along a beam path with the screen surface orthogonal to the beam path. This, together with the specified detector distance from the beam capsule center, automatically calculates the center vector ${\bf C}$ and the normal vector ${\bf n}$. If a beam number $\leq 0$ is specified, no alignment is performed and the user has to supply center vector and normal vector coordinates.
• pi_detectorDist2BeamCapsule_n: The detector distance from the beam capsule center, used only if the detector is alligned with respect to a beam.
• pi_detectorPinholeRadius_n: The radius of the pinhole. If $\leq 0$ no pinhole is used.
• pi_detectorPinholeDist2Det_n: The pinhole center distance from the detector center, useful only if a pinhole radius $>0$ was specified.

28.0.6.3 Proton Imaging General Runtime Parameters

• pi_cellStepTolerance: This factor times the smallest dimension of each cell is taken as the positional error tolerance for the protons during their path (parabolic or Runge-Kutta) tracing through each cell.
• pi_cellWallThicknessFactor: Controls the (imaginary) thickness of the cell walls to ensure computational stability of the proton imaging code. The cell thickness is defined as this factor times the smallest cell dimension along all geometrical axes. The factor is currently set to $10^{-6}$ and should only very rarely be changed.
• pi_detectorFileNameTimeStamp: If set to true, a time stamp will be added to each detector file name. This allows for time splitting of detectors.
• pi_detectorXYwriteFormat: Controls formatted ascii output of the proton $(x,y)$ pairs on the detector screens (default 'es20.10').
• pi_detectorDGwriteFormat: Controls formatted ascii output of the magnetic diagnostic quantities $k_x,k_y,k_z,J$ on the detector screens (default 'es15.5').
• pi_flagDomainMissingProtons: If true, any proton missing the domain will abort the program. If false, protons missing the domain will be recorded directly on the detector screen.
• pi_ignoreElectricalField: If true, the electrical field in each cell is ignored and the proton movement is only governed by the magnetic field (deflection only).
• pi_IOaddDetectorScreens: If true, the square frame of each detector will be added to the plotfile.
• pi_IOaddProtonsCapsule2Domain: If true, the proton paths from the capsule to the domain entry point will be added to the plotfile.
• pi_IOaddProtonsDomain2Screen: If true, the proton paths from the domain exit point to the detector screen will be added to the plotfile.
• pi_IOmaxBlockCrossingNumber: This is an estimate of the maximum proton path length while travelling through the block. This estimate is given as an integer number in units of block sides (default 5). Only in extremely rare cases, when the magnetic and electric fields are strong and sufficiently warped inside a block and the protons start to circulate inside the block, does one have to increase this number.
• pi_IOnumberOfProtons2Plot: The number of protons that are to be plotted and written to the plotfile.
• pi_maxProtonCount: The maximum number of beam/disk protons that can be created/read-in on one processor per batch in the domain.
• pi_opaqueBoundaries: If true, the protons do not go through cells marked as opaque boundaries.
• pi_printBeams: If true, it prints detailed information about the proton beams to a file with name $<$basename$>$ProtonBeamsPrint.txt, where $<$basename$>$ is the base name of the simulation.
• pi_printDetectors: If true, it prints detailed information about the proton detectors to a file with name $<$basename$>$ProtonImagingDetectors.txt, where $<$basename$>$ is the base name of the simulation.
• pi_printMain: If true, it prints general information regarding the proton imaging setup to a file with name $<$basename$>$ProtonImagingMainPrint.txt, where $<$basename$>$ is the base name of the simulation.
• pi_printProtons: If true, it prints detailed information about the all protons generated during the simulation. Protons are 'generated' on the domain surface and the info is written to several files labeled by batch number BID, processor rank number PID and a time stamp. Each processor writes its own file(s) with name(s): $<$basename$>$printProtonsBatch$<$BID$>$Proc$<$PID$>$.txt, where $<$basename$>$ is the base name of the simulation. Use of this feature is reserved ONLY for debugging purposes and is currently limited to 10 batches and 100 processors per time stamp. Usage of a larger number of batches/processors during a simulation does not abort the run, but protons on batches with BID $> 10$ and processors with PID $> 99$ are simply ignored and not printed. Users other than code developers should not activate this feature.
• pi_protonDeterminism: If true, the Grid Unit will be forced to use the sieve algorithm to move the proton particle data. Forcing this algorithm will result in a slower movement of data, but will fix the order the processors pass data and eliminate round off differences in consecutive runs.
• pi_recordOffScreenProtons: If true, protons not hitting their target detector screen are also recorded on the detector files. These protons will have screen coordinate pairs $(x,y)\notin [0,1]$.
• pi_relativisticTracing: If true, the motion of the protons through the domain will be treated relativistically using Eq.(28.16) for the relativistic acceleration and Eq.(28.29) for the relativistic jerk terms.
• pi_RungeKuttaMethod: Specifies which Runge Kutta method to use for proton tracing. Current options are: 'CashKarp45' (order 4, default), 'EulerHeu12' (order 1), 'BogShamp23' (order 2), 'Fehlberg34' (order 3) and 'Fehlberg45' (order 4).
• pi_screenProtonBucketSize: Sets the bucket size for flushing screen protons out to disk.
• pi_screenProtonDiagnostics: If true, the magnetic diagnostic quantities $k_x,k_y,k_z,J$ are evaluated for each proton and recorded on the detector screens.
• pi_timeResolvedProtonImaging: If true, it activates the time resolved proton imaging part of the code. Protons might need several time steps to cross the domain.
• pi_useIOprotonPlot: If true, protons are plotted to the plotfile for visualization pusposes.
• pi_useParabolicApproximation: If true, the code traces protons parabolically through cells for low $B$ / high $v$ combinations (section 28.0.1.2).
• useProtonImaging: If false, no proton imaging will be performed, even if the code was compiled to do so. Bypasses the need to rebuild the code.
• threadProtonTrace: If true, proton tracing through a block is threaded. This runtime parameter can only be set during setup of the code.

28.0.7 Unit Test

The unit test for the proton imaging unit consists in sending a ring of protons perpendicular onto either a uniform circular $B$ field or a radially outward uniform $E$ field and measuring the deflection (radial increase) of the ring as well as the magnitudes of the magnetic diagnostic quantities $k$ and $J$ from Eqs.(28.31) and (28.32). Due to the circular $B$ field or the radial outward $E$ field, each proton experiences a radial outward force and the radial increase of the protons can be calculated analytically. Exact analytical expressions for $k$ and $J$ are also possible.

28.0.7.1 Deflection of a Proton Ring by a Uniform Circular Magnetic Field

Consider the case of a uniform circular magnetic field with constant magnetic flux density $B$ pointing tangentially in clock- or anticlock-wise orientation around an axis, which we chose to be the $y$ axis. If a proton has only a $y$ component velocity $v_{0y}\neq 0$ initially, then it will be deflected radially outward or inward, depending on the clock- or anticlockwise orientation of ${\bf B}$. To simplify the calculations, we assume that the proton is located on the positive $x$ axis, such that at this point we have only a $z$ magnetic component, that is $B_z=\pm B$ and $B_x,B_y=0$. Further we deal with the non-relativistic situation first, as the extension to the relativistic case is trivial. We will use Eqs.(28.6-28.8) for that purpose. Since no electric field is present, we have that:

$$E_x,E_y,E_z,e_x,e_y,e_z,b_x,b_y =$$ 0 (28.47)

$$b_z = B_z/B = \pm 1$$ (28.48)

$$v_{0x},v_{0z} = 0.$$ (28.49)

Putting these values into Eqs.(28.6-28.8), the proton velocities become:

$$v_x(t) = \pm v_{0y} \sin[BQ_mt/c]$$ (28.50)

$$v_y(t) = v_{0y} \cos[BQ_mt/c]$$ (28.51)

$$v_z(t) =$$ 0 (28.52)

showing that the proton moves within the same $z$-plane. The $+$ sign in $v_x(t)$ applies, if $B_z=+B$ (anticlockwise orientation) and the $-$ sign, if $B_z=-B$ (clockwise orientation). Integration over time gives the position equations:

$$r_x(t) = r_{0x} \pm \frac{cv_{0y}}{BQ_m}(1-\cos[BQ_mt/c])$$ (28.53)

$$r_y(t) = r_{0y} + \frac{cv_{0y}}{BQ_m}\sin[BQ_mt/c]$$ (28.54)

$$r_z(t) = r_{0z}$$ (28.55)

Taking the anticlockwise orientation of $B$, looking from above onto the fized $z$-plane, the proton will perform a half-circle of radius $cv_{0y}/BQ_m$. We are interested in the (radial) deflection along the $x$ axis. The proton will travel through the domain in the $y$ direction, where on exit the detector screen is placed. If we place the detector right at the exit of the $y$ direction, then we can calculate from the $r_y(t)$ equation the time it takes for the proton to travel through (exit) the domain:

$$t_{exit} = \frac{c}{BQ_m}\arcsin\left(\frac{D_yBQ_m}{cv_{0y}}\right),$$ (28.56)

where $D_y$ is the length of the domain in $y$ direction. If $D_yBQ_m/cv_{0y}>1$, then the proton will not exit the domain but rather curve back inside the domain. This happens if either $D_y$ is too large (the domain extends too far into the $y$ direction), $B$ is too large or the initial velocity $v_{0y}$ is too low. For a given $D_y$ and $v_{0y}$, the limiting $B$ would be:

$$B_{limit} = \frac{cv_{0y}}{D_yQ_m}.$$ (28.57)

Inserting the $t_{exit}$ expression into the $r_x(t)$ equation, we get:

$$r_x(t_{exit}) = r_{0x} + \frac{cv_{0y}}{BQ_m} \left(1-\sqrt{1-\left(\frac{D_yBQ_m}{cv_{0y}}\right)^2}\right).$$ (28.58)

The radial deflection $\Delta r$ of the proton on the detector screen will be equal to $r_x(t_{exit})-r_{0x}$ and therefore:

$$\Delta r = \frac{cv_{0y}}{BQ_m}\left(1-\sqrt{1-\left(\frac{D_yBQ_m}{cv_{0y}}\right)^2}\right).$$ (28.59)

The limiting radial deflection (i.e., the maximal radial deflection) that can be detected by the screen at the domains exit is obtained by inserting the expression of $B_{limit}$ into the last equation, which gets:

$$\Delta r_{limit} = D_y,$$ (28.60)

and is hence equal to the domains extension in $y$ direction.

To evaluate the magnitude of the magnetic diagnostic quantities $k$ and $J$, we start from the differential equations derived from Eqs.(28.33) and (28.34):

$$\frac{dk}{dt} = Bv_{0y},$$ (28.61)

$$\frac{dJ}{dt} = v_y(t)\frac{\partial B_x}{\partial z} = - B\frac{v_y(t)}{r_x(t)},$$ (28.62)

where the first equation follows from the fact that the velocity magnitude does not change and the second equation utilized the fact that for a proton positioned on the $x$ axis only the term $\partial Bx/\partial z$ from $\nabla\times {\bf B}$ is nonzero (for anticlockwise rotation of $B$ we have $B_x=-zB/\sqrt{x^2+z^2}$, $B_y=0$ and $B_z=xB/\sqrt{x^2+z^2}$, which gives at $z=0$ locations $\partial B_x/\partial z=-B/x$). Performing the time integration from $t=0$ to $t=t_{exit}$ and using the expressions for $v_y(t)$, $r_x(t)$ and $t_{exit}$ from Eqs.(28.51), (28.53) and (28.56), we get:

$$k(t_{exit}) = BD_y\frac{\arcsin(x)}{x},$$ (28.63)

$$J(t_{exit}) = B\left(\arcsin(x)-2\frac{1+1/y}{\sqrt{1+2/y}} \arctan\left(\frac{x\sqrt{1+2/y}}{1+\sqrt{1-x^2}}\right)\right),$$ (28.64)

where we have used the dimensionless quantities:

$$x = \frac{D_yBQ_m}{cv_{0y}} , y = \frac{r_{0x}BQ_m}{cv_{0y}}.$$ (28.65)

Since the electrical field is zero, inclusion of relativity into the deflection equations is trivial. There is no increase in velocity as the proton travels inside the domain, which means that $\gamma({\bf v})$ stays constant. To account for relativistic effects, one simply has to replace $Q_m$ in the radial deflection and magnetic diagnostic equations by $Q_{m\gamma}$, where $\gamma({\bf v})=1/\sqrt{1-v_{0y}^2/c^2}$. For fixed parameters, the relativistic radial deflection is always less than the non-relativistic one. This can be shown by evaluating the ratio between the two radial deflections. Using the non-relativistic $x$ as defined in Eq.(28.65), we have:

$$\frac{{\Delta r}_{non-rel}}{{\Delta r}_{rel}} = \frac{1-\sqrt{1-x^2}}{\gamma - \sqrt{\gamma^2-x^2}} \;=\; \gamma\frac{1+\sqrt{1-(x/\gamma)^2}}{1+\sqrt{1-x^2}},$$ (28.66)

and, since $x^2\leq 1$ and $\gamma\geq 1$, it follows that:

$$\frac{{\Delta r}_{non-rel}}{{\Delta r}_{rel}} \geq \gamma.$$ (28.67)

28.0.7.2 Deflection of a Proton Ring by a Uniform Radial Outward Electric Field

In this section we state the radial deflection equations in the presence of a radial outward directed electric field. To simplify again the calculations, we assume that the proton is located such that at this point we have only a $x$ electric component, that is $E_x=\pm E$ and $E_y,E_z=0$. To include the relativistic case, we start from the relativistic acceleration in Eq.(28.16) and solve the set of differential equations corresponding to the unit test setup:

$${\dot v}_x(t) = Q_{m\gamma}[E - v_x^2(t)E/c^2]$$ (28.68)

$${\dot v}_y(t) = Q_{m\gamma}[- v_x(t)v_y(t)E/c^2]$$ (28.69)

$${\dot v}_z(t) = Q_{m\gamma}[- v_x(t)v_z(t)E/c^2]$$ (28.70)

with initial conditions $v_y(0)=v_{0y}$ and $v_x(0),v_z(0)=0$. The results are:

$$v_x(t) = EQ_{m\gamma}t/\sqrt{1+\delta (EQ_{m\gamma}t/c)^2}$$ (28.71)

$$v_y(t) = v_{0y}/\sqrt{1+\delta (EQ_{m\gamma}t/c)^2}$$ (28.72)

$$v_z(t) =$$ 0 (28.73)

where $Q_{m\gamma}$ is evaluated using $\gamma({\bf v})=1/\sqrt{1-v_{0y}^2/c^2}$ with $Q_{m\gamma}=Q_m$ for the non-relativistic case, and $\delta$ is 1 for the relativistic case and 0 otherwise. Integration over time gives (note the relativistic corrections to the classical expressions):

$$r_x(t) = r_{0x} + \frac{EQ_{m\gamma}t^2}{2}\left[2\left(\frac{EQ_{m\gamma}... ... \left(\sqrt{1+\left(\frac{EQ_{m\gamma}t}{c}\right)^2}-1\right)\right]^{\delta}$$ (28.74)

$$r_y(t) = r_{0y} + v_{0y}t\left[\left(\frac{EQ_{m\gamma}t}{c}\right)^{-1}\operatorname{arsinh} \left(\frac{EQ_{m\gamma}t}{c}\right)\right]^{\delta}$$ (28.75)

$$r_z(t) = r_{0z}$$ (28.76)

From the equation for $r_y(t)$ we get the domain crossing time in the $y$ direction as:

$$t_{exit} = \frac{D_y}{v_{0y}}\left[\left(\frac{D_yEQ_{m\gamma}}{cv_{0y}}\right)^{-1} \sinh\left(\frac{D_yEQ_{m\gamma}}{cv_{0y}}\right)\right]^{\delta}$$ (28.77)

and the resulting radial reflection $r_x(t_{exit})-r_{0x}$ is:

$$\Delta r = \frac{D_y^2EQ_{m\gamma}}{2v_{0y}^2} \left[\left(\frac{D_yEQ_{m\ga... ...ght)^{-2} \sinh^2\left(\frac{D_yEQ_{m\gamma}}{2cv_{0y}}\right)\right]^{\delta}.$$ (28.78)

In contrast to the magnetic field radial deflection, there is, for a fixed $E$ a certain $v_{0y}$ and for a fixed $v_{0y}$ a certain $E$, beyond which the relativistic radial deflection is greater than the non-relativistic one. Using the simplification:

$$x = \frac{D_yEQ_{m\gamma}}{2cv_{0y}},$$ (28.79)

we have for the ratio:

$$\frac{{\Delta r}_{non-rel}}{{\Delta r}_{rel}} = \frac{\gamma x^2}{\sinh^2(x)}.$$ (28.80)

But the function $f(x)=x^2/\sinh^2(x)$ is a monotonic decaying function between $f(x)=1$ for $x = 0$ and $f(x)=0$ for $x\rightarrow \infty$. Hence at some $x$ we must have $f(x)=1/\gamma$, beyond which we have ${\Delta r}_{rel}>{\Delta r}_{non-rel}$. The crossover radial deflecion ${\Delta r}_{cross}$, at which both the non-relativistic and the relativistic radial reflections coincide, can be calculated by finding the solutions for $x_{cross}$ and ${\Delta r}_{cross}$ in the following two equations:

$$\gamma x_{cross}^2 = \sinh^2(x_{cross})$$ (28.81)

$${\Delta r}_{cross} = \gamma x_{cross} D_y c / v_{0y}.$$ (28.82)

Note that the first equation, defining $x_{cross}$ as a function of $\gamma$, cannot be solved analytically and one must resort to numerical approaches.

28.0.7.3 FLASH Setup of the Unit Test

We set up a 3D cartesian domain consisting of a rectangular box with dimensions (in cm): $x$ axis $[0,3]$, $y$ axis $[0,1]$ and $z$ axis $[0,3]$. A circular beam of 10000 protons with radius 0.5cm is shot along the $y$ axis perpendicular to the $xz$ plane and centered on the $xz$ plane at (1.5,1.5). The square detector screen was placed extremely close to the domain exit in $y$ direction and was of the same size as the domain's $xz$ plane, i.e., of side length 3cm and centered at (1.5,1.5). The variable parameters used were: 1) the proton energies (currently only 20MeV), 2) the uniform refinement levels (3,4,5) and the magnetic / electric fields (in Gauss). The following radial deflections $\Delta r$ were recorded: 1) the theoretical $\Delta r_{thr}$, 2) the maximum $\Delta r_{max}$, 3) the minimum $\Delta r_{min}$ and 4) the average $\Delta r_{avg}$.

28.0.7.4 FLASH Unit Test Results for Uniform Circular Magnetic Fields

From the FLASH setup parameters, we can calculate the following values for 20 MeV protons:

$$v_{0y} = 0.2032392740948\;c$$

$$\gamma = 1.0213157779052$$

$$B_{limit} = 636085.837576\;$$

$${\Delta r}_{limit} = 1.0\; .$$

Tables 28.1, 28.2 and 28.3 show the obtained FLASH results for the radial deviation, the magnetic diagnostic $J$ value and the magnetic diagnostic ${\bf k}$ magnitude, including a run with $B_{limit}$.

Table: Magnetic outward radial deflections: maximum $\Delta r_{max}$, minimum $\Delta r_{min}$, average $\Delta r_{avg}$ and theoretical $\Delta r_{thr}$ from Eq.(28.59) on a ring of 20 MeV protons for uniform domain refinement levels 3,4,5. Subscript numbers indicate multiplication by negative exponents: $n_e$ means $n\times 10^{-e}$. The radial deflection values are in cm. Italic numbers represent the relativistic results.

Proton energy 20 MeV (0.203 c), Refinement Levels = 3,4,5
$B$	$\Delta r_{max}$	$\Delta r_{min}$	$\Delta r_{avg}$	$\Delta r_{thr}$
$10,000$	$7.86106_3$	$7.85565_3$	$7.85991_3$
	$7.86106_3$	$7.85995_3$	$7.86075_3$	$7.86106_3$
	$7.86106_3$	$7.86081_3$	$7.86098_3$
	$\it 7.69697_3$	$\it 7.69166_3$	$\it 7.69585_3$
	$\it 7.69697_3$	$\it 7.69588_3$	$\it 7.69667_3$	$\it 7.69697_3$
	$\it 7.69697_3$	$\it 7.69672_3$	$\it 7.69690_3$
$100,000$	$7.90975_2$	$7.90635_2$	$7.90862_2$
	$7.90975_2$	$7.90891_2$	$7.90946_2$	$7.90975_2$
	$7.90975_2$	$7.90956_2$	$7.90968_2$
	$\it 7.74266_2$	$\it 7.73932_2$	$\it 7.74154_2$
	$\it 7.74266_2$	$\it 7.74183_2$	$\it 7.74237_2$	$\it 7.74266_2$
	$\it 7.74266_2$	$\it 7.74246_2$	$\it 7.74258_2$
$636,085$	$9.98376_1$	$9.80750_1$	$9.86908_1$
	$9.98376_1$	$9.90287_1$	$9.93061_1$	$9.98375_1$
	$9.98376_1$	$9.94942_1$	$9.96196_1$
	$\it 8.13739_1$	$\it 8.12900_1$	$\it 8.13357_1$
	$\it 8.13739_1$	$\it 8.13531_1$	$\it 8.13638_1$	$\it 8.13739_1$
	$\it 8.13739_1$	$\it 8.13687_1$	$\it 8.13713_1$

Table: Magnetic diagnostic values: maximum $J_{max}$, minimum $J_{min}$, average $J_{avg}$ and theoretical $J_{thr}$ from Eq.(28.64) on a ring of 20 MeV protons for uniform domain refinement levels 3,4,5. Subscript numbers indicate multiplication by positive exponents: $n_e$ means $n\times 10^{+e}$. All $J$ values are in Gauss. Italic numbers represent the relativistic results.

Proton energy 20 MeV (0.203 c), Refinement Levels = 3,4,5
$B$	$J_{max}$	$J_{min}$	$J_{avg}$	$J_{thr}$
$10,000$	$-1.93499_4$	$-2.06117_4$	$-1.99244_4$
	$-1.96351_4$	$-2.01282_4$	$-1.98985_4$	$-1.98962_4$
	$-1.98027_4$	$-1.99475_4$	$-1.98977_4$
	$\it -1.93499_4$	$\it -2.06117_4$	$\it -1.99272_4$
	$\it -1.96351_4$	$\it -2.01295_4$	$\it -1.99008_4$	$\it -1.98983_4$
	$\it -1.98027_4$	$\it -1.99504_4$	$\it -1.98997_4$
$100,000$	$-1.88339_5$	$-1.92096_5$	$-1.90503_5$
	$-1.89600_5$	$-1.90952_5$	$-1.90412_5$	$-1.90377_5$
	$-1.90084_5$	$-1.90575_5$	$-1.90385_5$
	$\it -1.88506_5$	$\it -1.92296_5$	$\it -1.90692_5$
	$\it -1.89761_5$	$\it -1.91166_5$	$\it -1.90596_5$	$\it -1.90563_5$
	$\it -1.90277_5$	$\it -1.90757_5$	$\it -1.90571_5$
$636,085$	$-9.58761_5$	$-9.69203_5$	$-9.64623_5$
	$-9.62174_5$	$-9.65740_5$	$-9.64229_5$	$-9.64102_5$
	$-9.63384_5$	$-9.64623_5$	$-9.64133_5$
	$\it -9.65418_5$	$\it -9.75959_5$	$\it -9.71335_5$
	$\it -9.68864_5$	$\it -9.72467_5$	$\it -9.70939_5$	$\it -9.70812_5$
	$\it -9.70087_5$	$\it -9.71338_5$	$\it -9.70844_5$

Table: Magnetic diagnostic values: maximum $k_{max}$, minimum $k_{min}$, average $k_{avg}$ and theoretical $k_{thr}$ from Eq.(28.63) on a ring of 20 MeV protons for uniform domain refinement levels 3,4,5. Subscript numbers indicate multiplication by positive exponents: $n_e$ means $n\times 10^{+e}$. All $k$ values are in units of Gauss cm. Italic numbers represent the relativistic results.

Proton energy 20 MeV (0.203 c), Refinement Levels = 3,4,5
$B$	$k_{max}$	$k_{min}$	$k_{avg}$	$k_{thr}$
$10,000$	$1.00004_4$	$9.99561_3$	$1.00001_4$
	$1.00004_4$	$9.99924_3$	$1.00003_4$	$1.00004_4$
	$1.00004_4$	$1.00001_4$	$1.00004_4$
	$\it 1.00004_4$	$\it 9.99560_3$	$\it 1.00001_4$
	$\it 1.00004_4$	$\it 9.99922_3$	$\it 1.00003_4$	$\it 1.00004_4$
	$\it 1.00004_4$	$\it 1.00001_4$	$\it 1.00003_4$
$100,000$	$1.00417_5$	$1.00378_5$	$1.00406_5$
	$1.00417_5$	$1.00407_5$	$1.00413_5$	$1.00417_5$
	$1.00417_5$	$1.00414_5$	$1.00416_5$
	$\it 1.00399_5$	$\it 1.00360_5$	$\it 1.00388_5$
	$\it 1.00399_5$	$\it 1.00389_5$	$\it 1.00396_5$	$\it 1.00399_5$
	$\it 1.00400_5$	$\it 1.00397_5$	$\it 1.00398_5$
$636,085$	$9.98129_5$	$9.86877_5$	$9.90802_5$
	$9.98129_5$	$9.92972_5$	$9.94739_5$	$9.98129_5$
	$9.98129_5$	$9.95942_5$	$9.96739_5$
	$\it 8.87495_5$	$\it 8.86913_5$	$\it 8.87219_5$
	$\it 8.87495_5$	$\it 8.87350_5$	$\it 8.87422_5$	$\it 8.87496_5$
	$\it 8.87495_5$	$\it 8.87459_5$	$\it 8.87477_5$

28.0.7.5 FLASH Unit Test Results for Uniform Radial Outward Electric Fields

From the FLASH setup parameters, we can calculate the following values for 20 MeV protons:

$$v_{0y} = 0.2032392740948\;c$$

$$\gamma = 1.0213157779052$$

$$x_{cross} = 0.2518110572090\;$$ (28.83)

$$E_{cross} = 327175.337727\;$$

$${\Delta r}_{cross} = 1.2653981713131\; .$$

Table 28.4 shows the obtained FLASH results, including a run with $E_{cross}$.

Table: Electric outward radial deflections: maximum $\Delta r_{max}$, minimum $\Delta r_{min}$, average $\Delta r_{avg}$ and theoretical $\Delta r_{thr}$ from Eq.(28.78) on a ring of 20 MeV protons for uniform domain refinement levels 3,4,5. Subscript numbers indicate multiplication by negative exponents: $n_e$ means $n\times 10^{-e}$. The radial deflection values are in cm. Italic numbers represent the relativistic results.

Proton energy 20 MeV (0.203 c), Refinement Levels = 3,4,5
$E$	$\Delta r_{max}$	$\Delta r_{min}$	$\Delta r_{avg}$	$\Delta r_{thr}$
$10,000$	$3.86765_2$	$3.86592_2$	$3.86708_2$
	$3.86765_2$	$3.86722_2$	$3.86750_2$	$3.86765_2$
	$3.86765_2$	$3.86754_2$	$3.86761_2$
	$\it 3.78700_2$	$\it 3.78530_2$	$\it 3.78644_2$
	$\it 3.78700_2$	$\it 3.78658_2$	$\it 3.78685_2$	$\it 3.78700_2$
	$\it 3.78700_2$	$\it 3.78689_2$	$\it 3.78696_2$
$100,000$	$3.86765_1$	$3.86647_1$	$3.86718_1$
	$3.86765_1$	$3.86736_1$	$3.86752_1$	$3.86765_1$
	$3.86765_1$	$3.86757_1$	$3.86761_1$
	$\it 3.79441_1$	$\it 3.79325_1$	$\it 3.79394_1$
	$\it 3.79441_1$	$\it 3.79412_1$	$\it 3.79429_1$	$\it 3.79441_1$
	$\it 3.79441_1$	$\it 3.79434_1$	$\it 3.79438_1$
$327,176$	$1.26540_0$	$1.24971_0$	$1.25927_0$
	$1.26540_0$	$1.24978_0$	$1.25936_0$	$1.26540_0$
	$1.26540_0$	$1.24980_0$	$1.25938_0$
	$\it 1.26540_0$	$\it 1.24893_0$	$\it 1.25896_0$
	$\it 1.26540_0$	$\it 1.24899_0$	$\it 1.25904_0$	$\it 1.26540_0$
	$\it 1.26540_0$	$\it 1.24901_0$	$\it 1.25907_0$