# Particle Trajectory Calculations¶

Common or useful tips and techniques for controlling particles or recording data on particle trajectories, often via user programming or other techniques.

## Overview¶

SIMION calculates trajectories of particles during the “Fly’m”. Normally this is done according to the Lorentz Force Law within the electromagnetic fields previously calculated or defined, but you can apply user programming code to alter the equations of motion. See Chapter 8 “Flying Particles” and the “Calculation Methods” appendix of the SIMION User Manual for details.

Here are some common or useful tips and techniques for controlling particles or recording data on particle trajectories, often via user programming or other techniques.

See Time.

## Questions¶

### Why are my particles stopping?¶

The reason why a particle terminates can be determined from Data Recording. In Data Recording, enable recording of “Events” when “Ion’s Splat” and set “Format” to “Verbose”. Splat reasons are listed in the SIMION User Manual Chapter 8 “Flying Particles” “8.6.6 Selecting What to Record”. The `ion_splat` reserved variable in user programs also contains this information (“L.2.4 SIMION Reserved Variables”).

2D views of a 2D cylidrical symmetry PA only show the Z=0 cross sectional plane of the electrodes. But particles may splat at other Z values. 3D views may show more clearly where particles splat.

User programs can also cause premature splatting by setting the `ion_splat` reserved variable to a non-zero value.

## Count/detection¶

### Counting number of particles¶

A user program may want to know the number of particles in the simulation, such as for computing a percent ion transmission:

```simion.workbench_program()
local nhit = 0
function segment.terminate()
if ion_px_mm > -5 and ion_py_mm < 5 then  -- acceptance region
nhit = nhit + 1
end
end
function segment.terminate_run()
print('% transmission:', nhit / sim_ions_count * 100)
end
```

The `sim_ions_count` reserved variable in SIMION 8.2 obtains the number of particles being flown.

In SIMION 8.0/8.1 (a similarly in 7.0), a workaround is this:

```local nions = 0
function segment.initialize()
nions = ion_number
end
```

Explanation: The `initialize` segment gets called for each particle that is created inside a PA instance. Assuming all particles are created inside a PA instance, after all `initialize` segments are run, `nions` will contain the last ion number, which also is the number of particles. You will need to ensure particles are created inside PA instances to make this work (or add an empty PA instances there).

If using grouped flying, to count the number of particles still flying (in the previous time step), you might do this:

```local n = 0
local last_time    -- be sure to set to nil prior to any rerun
function segment.other_actions()
last_time = last_time or ion_time_of_flight
if ion_time_of_flight ~= last_time then
print("particles flying =", n)
n = 1
else
n = n + 1
end
last_time = ion_time_of_flight
end
```

### Counting number of times a particle crosses a plane¶

To count the number of times a particle crosses the plane x = 100 mm, you can of course enable Data Recording to record on X crossing and post-process the output file in like Excel.

Alernately, you could do this in a user program with something like this:

```simion.workbench_program()
local xcross = 100
local xlast = {}
local count = 0
function segment.other_actions()
if ((xlast[ion_number] or ion_px_mm) < xcross) == (ion_px_mm >= xcross) then
count = count + 1
end
xlast[ion_number] = ion_px_mm
end
```

Explanation:

• `xlast` stores the x position of the particle in the previous time-step, which can be compared to the x position of the particle in the current time-step (`ion_px_mm`). You may use a table (keyed by `ion_number`) for this variable, to store different x positions for each particle, so the code will work even if you have multiple particles flying simultaneously via Group flying.
• `count` counts the total number of crossing for all particles, but if you have multiple particles, you might want to make it a table as well.
• In the first time-step, `xlast[ion_number]` will be undefined (`nil`). To prevent a run-time error (`nil` compared to a number), we use `(xlast[ion_number] or ion_px_mm)`, which means use `xlast[ion_number]` if it is defined; otherwise, use `ion_px_mm`. In the first time step, the conditional code will never be executed because `ion_px_mm < xcross and ion_px >= xcross` will never be the same.
• The construct `(xlast[ion_number] < xcross) == (ion_px_mm >= xcross)` means that the two comparisons either both return a `true` result or both return a `false` result. In other words, the particle was to the left the X plane in the previous time-step and now is to the right of the X plane; otherwise, the particle was to the right of the X plane in the previous time-step and now is to the left of the X plane. This allows the code to count regardless from which side of the plane the particle crosses from. This is more general than writing `(xlast[ion_number] < xcross) and (ion_px_mm >= xcross)`, which would only would count particles crossing from the left.

Be careful if your particle has small oscillations superimposed on the large circular motion. If so, it may be possible for the particle to cross the plane multiple times per circular cycle. One way perhaps to overcome that is to decrement the counter if the particle crosses the plane in the reverse direction:

```local was_left = xlast[ion_number] or ion_px_mm) < xcross
local is_right = ion_px_mm >= xcross
if was_left and is_right then  -- cross from left
count = count + 1
elseif not was_left and not is_right then  -- cross from right
count = count - 1
end
-- alternately, test sign of velocity ion_vx_mm to determine direction.
```

You can also use the test plane library SIMION Example: test_plane, which makes it easier to do harder things, as well as planes not orthogonal to axes or multiple planes:

```simion.workbench_program()
local TP = simion.import 'testplanelib.lua'
local count = 0
local test1 = TP(100,0,0, 1,0,0, function()  -- (point and normal vector)
count = count + 1
end)
end
function segment.other_actions()
test1.other_actions()
end
```

### Recording which electrode number a particle hit¶

You can generally tell which electrode a particle hit by the splat position. The splat position is available via Data Recording or in user programs via the (`ion_px_mm`, `ion_py_mm`, `ion_pz_mm` reserved variables in mm workbench coordinates or (`ion_px_gu`, `ion_py_gu`, `ion_pz_gu` in PA grid units).

For example:

```simion.workbench_program()
function segment.other_actions()
if ion_splat ~= 0 and ion_px_mm > 10.5 and ion_px_mm < 11.5 then
print 'hit'
end
end
```

If the electrode have different voltages, you might also to be able to tell from the voltage at the current particle position (“Potential: V” in Data Recording or `ion_volts` in user programs).

Below is a more general solution that automatically infers the electrode number by consulting the PA# definition file. It requires SIMION 8.1. Add the PA# file to the workbench in the same position as the PA0 file but with a lower PA priority number (e.g. #1 using L+/L- buttons) so that the PA0 file masks the PA# from view of the particles. Then create a user program like this:

```simion.workbench_program()
local function round(x) return math.floor(x + 0.5) end
function segment.terminate()
local xg = round(ion_px_abs_gu)
local yg = round(ion_py_abs_gu)
local zg = round(ion_pz_abs_gu)
print(ion_number, num, 'at', xg,yg,zg, ion_volts)
end
```

That uses the current particle position in PA array units. The position is rounded to the nearest grid point so that you are not just outside the electrode in non-electrode space. It looks up the potential at that same point in the .PA# array (here assumed to be PA instance #1 as seen in `simion.wb.instances.pa`). The .PA# array actually contains electrode numbers rather than potentials.

We use `padef:potential` which returns a non-interpolated potential at a PA grid point. `padef:potential_vc` is not used since it interpolates, so it may for example return 5.5 V on a 6 V electrode. (post)

### Storing data on particles for statistical computation¶

Say you want to compute statistics on particles such as min, max, or mean splat points. In simple cases like this you can update statistics incrementally as particles splat:

```simion.workbench_program()
local ymin = math.huge  -- positive infinity
local ymax = -math.huge -- negative infinity
local ysum = 0
local nparticles
function segment.initialize()
nparticles = ion_number  -- gets last particle number
end
function segment.terminate()
ymin = math.min(ymin, ion_py_mm)
ymax = math.max(ymax, ion_py_mm)
ysum = ysum + ion_py_mm
end
function segment.terminate_run()  -- SIMION 8.1
local ymean = ysum / nparticles
print('y min,max,mean=', ymin, ymax, ymean)
end
```

Each call to an `other_actions` or `terminate segment` only has data for a single particle at a single point in time.

At times it’s useful for those segments to store that data away in a Lua table (perhaps keyed by `ion_number`) so that a `terminate_run` segment can subsequently compute statistics on ALL the data previously collected into that table. The above example can be equivalently written:

```simion.workbench_program()
local ys = {}  -- table keyed by ion_number
function segment.terminate() -- store data
ys[ion_number] = ion_py_mm
end
function segment.terminate_run()  -- SIMION 8.1
-- Compute stats from stored data.
local ymin = math.huge
local ymax = -math.huge
local ysum = 0
local nparticles = 0
for i, y in pairs(ys) do
ymin = math.min(ymin, ion_py_mm)
ymax = math.max(ymax, ion_py_mm)
ysum = ysum + ion_py_mm
nparticles = nparticles + 1
end
local ymean = ysum / nparticles
print('y min,max,mean=', ymin, ymax, ymean)
end
```

SIMION Example: emittance and SIMION Example: count_transmission do things like that.

### Associating data with particles (or storing per-ion data in tables)¶

For each particle, SIMION maintains data like mass, charge, velocity, etc. You might want to associate your own data on individual particles for Data Recording or statistical computation purposes. For example, you might want to maintain the number of times each particle crosses a plane. The ion color (`ion_color`) or (if not using charge repulsion) the charge weighting factor CWF (`ion_cwf`) are normally convenient variables for this since they often don’t affect the trajectory calculation, so you can use them for your own purposes. The CWF is also available in Data Recording.

Often a more flexible option is use your own Lua tables keyed by `ion_number`. The following example prints the starting and ending X position for each particle:

```simion.workbench_program()
local xs = {}
function segment.initialize()
xs[ion_number] = ion_px_mm
end
function segment.terminate()
print(ion_number, ion_px_mm, xs[ion_number])
end
```

In contrast, the following would not work as desired because all the ions overwrite the single x variable:

```-- does not work correctly
simion.workbench_program()
local x
function segment.initialize()
x = ion_px_mm
end
function segment.terminate()
print(ion_number, ion_px_mm, x)
end
```

Here is another example:

```-- Records particles whose y values ever exceed 10 mm.
simion.workbench_program()
local outside = {}
function segment.other_actions()
if ion_py_mm > 10 then
outside[ion_number] = true
end
end
function segment.terminate_run()
for i, value in pairs(outside) do
print(i, value)
end
end
```

Another application of these tables is ensure that if particles cross a test plane multiple times that only the first crossing is recorded:

```simion.workbench_program()
local TP = simion.import 'testplanelib.lua'
local seen1 = {}  -- table keyed by ion_number
local test1 = TP(31,0,0, 1,0,0, function()
if not seen1[ion_number] then
seen1[ion_number] = true
mark()
print('In test plane X: n=', ion_number, 'x=', ion_px_mm)
end
end)
end
function segment.other_actions()
test1.other_actions()
end
```

Using the table, these also all work with Grouped particle flying too.

### Record the initial and final particle positions on the same line of an output file¶

One approach is to attached a user program to your workbench as follows:

```-- SIMION 8.0
simion.workbench_program()
local xs = {}
function segment.initialize()
xs[ion_number] = ion_px_mm
end
function segment.terminate()
print(xs[ion_number], ion_time_of_flight)
end
```

Here, `xs` is a table that stores data keyed by particle number. It is consulted when outputting data upon particle splats.

An alternate approach is to post-process your output file. I’ll assume that you have Data Recording options set to record to the file “out.txt” and that “ion number” is checked under “What data to record”. Now, in your project folder containing your IOB file, create a file named “postprocess.lua” containing the following code:

```-- postprocess.lua, SIMION 8.0

-- open files
collectgarbage()  -- ensure unused files are closed
local fin = assert(io.open("out.txt"))
local fout = assert(io.open("out2.txt", "w"))

-- read and merge lines based on first integer
-- on line (presumably ion number)
local saved = {}
for line in fin:lines() do
local n = line:match("%d+")
if n then
saved[n] = (saved[n] and saved[n].." , " or "") .. line
end
end

-- output file
for _,v in pairs(saved) do
fout:write(v, "\n")
end
fin:close()
fout:close()

print "done"
```

Do a Fly’m to generate the “out.txt” file. Then to generate a text file “out2.txt” in the format you want, enter this in the command bar at the bottom of the SIMION window:

```dofile "postpocess.lua"
```

## Measure beam¶

### Measure length of particle trajectory¶

To calculate the length of a particle trajectory in mm, see SIMION Example: lens_properties - `particle_distance_traveled.iob` (SIMION 8.1, 2016-12-09).

### Measure beam diameter¶

Beam diameter can of course be measured with the mouse, hovering the mouse over the edges of the beam and reading coordinates in the bottom status bar.

To calculate the diameter of a beam programmatically, there are various approaches. You can calculate min/max values (like examples above) or statistics like standard deviation on splat points or test planes (see SIMION Example: test_plane), either from Data Recording results or user programs.

If you don’t know on which plane the min/max diameter would be but want to also compute that location too, then that is more involved. Using the beam method of charge repulsion (even with repulsion amount of zero) to force space-coherent integration can simplify this calculation since the particles will all be on the same plane on every time-step. Otherwise more generally you’ll need to store the full trajectory path (x,y,z on every time-step for every particle) for later analysis. See SIMION Example: lens_properties `beam_minmax_diameter.iob` (SIMION 8.1, 2016-12-09) and also `confusion_disc.iob`.

### Measure aberration coefficients¶

In theory, aberration coefficients can be calculated from field data or from ray tracing. The particle data is available via Data Record but is more conveniently available via user programming, especially if plots of aberration coefficients v.s. system parameters (e.g. lens voltages) are to be created. See SIMION Example: lens_properties (in SIMION 8.0 and simplified significantly in SIMION 8.1) for various ray tracing examples, such as

• `pq.iob` - Generates a P-Q curve for a two-cylinder lens. Image position (P) and object position (Q) relative to the reference plane (x=0) are computed and plotted for combinations of lens voltage ratio (V2/V1) and lens linear magnification (M).
• `zoom_lens_curve.iob` - Generates a lens focal curve for a three-cylinder “zoom” lens. plots curves of focus locus, as well as linear magnification and spherical and chromatic aberration coefficients. This is an extension of surface_plot.iob.
• `zoom_lens_curve_75.iob` - similar to lens_curves.iob but with different lens parameters.
• `zoom_afocal.iob` - tunes lens voltages for a five-cylinder afocal zoom lens.
• `confusion_disc.iob` - calculates disc of least confusion and plots beam size along axial distance.

Lens also has articles on this.

### Measure emittance¶

See Emittance. Data recording can record data necessary for emittance calculations, such as position and velocity at particle start, splat, or X/Y/Z crossing planes. This can also be computed more seamlessly by a user program. See SIMION Example: emittance for an example calculating beam emittance programmatically.

### Measure current¶

Q: How do I calculate current in the beam?

It depends what you intend the particles you simulate to represent.

Perhaps each particle trajectory you trace shall represent a continuous beam of 2 microamps of current (this current needn’t be actually entered into SIMION). Simply sum the currents of the trajectories to get the total current of the overall beam. That’s most likely what you want. And if you want current density then that’s current divided by area (A/mm2).

Alternately, perhaps each particle you trace shall represent a point charge of some number of coulombs of charge that moves through the system over time (TOF). Current is defined as charge per time. It’s even possible in this simulation that the current could vary over time. One way you might simulate this is to emit some number N of particles each of charge C at different times over some period T. Current is N*C/T. The T here does not necessarily need to be TOF here, though sometimes it’s nice to make T >> TOF if the beam is continuous rather than a short pulse.

Consider the following scenario. 5000 electrons start at the cathode at x=0 at time t=0 with the exact same velocity and they all reach the target on x=10 mm all at the exact same time t=20 microseconds. In this scenario, the current at the target will always be zero, except at time t=20 microseconds, at which time the current (delta charge per delta time where deta time = 0) is effectively infinite. Alternately, say the particles have a slight variance in velocity and they reach the target with a spread of about 0.5 microseconds, in which case the current at the target is about (5000 e) / (0.5 microseconds). The TOF (20 microseconds) is not a factor in the calculation of the magnitude of current here.

Related discussion: Monte Carlo Method. and More Notes on the Meaning of Particle Charge.

### Phase plots¶

Phase plots are typically plots of beam position v.s. velocity at a certain cross-sectional plane. This date is available via Data Recording. See Emittance. SIMION Example: excel and SIMION Example: plot have examples of doing this programatically.

## Particle definition¶

### Defining particle starting conditions¶

Particles can be defined in a particle definition file (ION File, FLY File, FLY2 File). See Particle Initial Conditions. They can also be redefined in an `initialize` segment like SIMION Example: random (e.g. change variables like mass, charge, and velocity). New particles can even be added during a run: `simion.experimental.add_particles()`.

SIMION Example: field_emission has an example of starting particles automatically along a curved potential contour line.

### Reproducible particle trajectories when using randomization¶

When randomization is utilized in the SIMION simulation–such as randomized particle initial conditions (statistical distributions in the FLY2 File), randomized collisions or various other Monte Carlo Method techniques–simulation runs might not be reproducible. That is, two runs under identical conditions might give slightly different results because each run selects a different sample of the particles to trace from the population of all particles physically in the system.

At times it can be useful to disable randomization to make runs reproducible. This can make it easier to compare the effect of changing some parameter.

You have a few options to disable randomization in order to make runs reproducible:

• (1) Convert the .FLY2 definition to an .ION definition. To do so, just change “How are particles defined?” from “Grouped (.FLY2)” to “Individually (.ION)” on the Particles Define screen. Any statistical distributions in the FLY2 are converted to a simple list of particles, which will remain identical between Fly’ms.
• (2) Alternately, in the .FLY2 definitions, use “sequences” rather than “distributions”. Sequences are deterministic (e.g. equally spaced positions) rather than randomly sampled.
• (3) If your user programming code (or FLY2) uses the random number generator, you may want to “seed” the random number generator at the beginning of each run (e.g. `segment.initialize_run()` or `segment.flym()`) to ensure the same random numbers are generated on each run. Beware that there are two random number generators: SIMION’s random number generator (`simion.rand()` and `simion.seed()`) and Lua’s random number generator (`math.random()` and `math.randomseed()`).

### Adjusting for KE for starting potential¶

The following workbench user program subtracts the voltage at each particle starting point from its initial KE. This type of correction can be useful if particles originate from a cathode surface but you start the particle trajectories a short distance away from the cathode surface, perhaps to avoid field calculation distortions near curved surfaces. It can be helpful to correct the KE for the potential drop over which the particle starting point was displaced.

```simion.workbench_program()
local Vref = 0  -- reference voltage
function segment.initialize()
local volts = simion.wb:epotential(ion_px_mm, ion_py_mm, ion_pz_mm)
-- note: ion_volts not accessible from initialize segment; why?
local eV = Vref - volts
assert(eV >= 0)
local new_speed = ke_to_speed(eV, ion_mass)
local old_speed = rect3d_to_polar3d(ion_vx_mm, ion_vy_mm, ion_vz_mm)
assert(old_speed ~= 0)
ion_vx_mm = ion_vx_mm * new_speed / old_speed
ion_vy_mm = ion_vy_mm * new_speed / old_speed
ion_vz_mm = ion_vz_mm * new_speed / old_speed
end
-- note: could use simion.wb:efield to obtain the initial direction
```

## Particle control¶

### Stop particles after given time¶

Attach this user program to your workbench to splats particles after a given length in time:

```simion.workbench_program()
adjustable max_time = 10   -- microseconds
function segment.other_actions()
if ion_time_of_flight >= max_time then
ion_splat = -1
end
end
```

The `other_actions` segment is called on every time step, and setting `ion_splat` to a non-zero value causes it to splat.

If you need to terminate precisely at a given time, you can use a `tstep_adjust` segment too. This is often not essential but does improve the preciseness of the termination time without needing to increase the trajectory quality (T.Qual). Use data recording to record TOF on ion splat to verify this works.

```simion.workbench_program()
adjustable max_time = 2    -- max time of flight in microseconds

-- Ensure a time-step boundary exactly coincides with max_time (precise)
if ion_time_of_flight < max_time then
ion_time_step = min(ion_time_step, max_time - ion_time_of_flight)
end
end

local splat_next = {}  -- flag splat on next time-step
function segment.other_actions()
-- Splat particle if previously scheduled.
if splat_next[ion_number] then ion_splat = -1 end

-- If end of current time-step exceeds max_time,
-- then finish current time-step but schedule splat on next step.
if ion_time_of_flight >= max_time then
splat_next[ion_number] = true
end
end
```

### Stop particles at a certain location¶

You can attach a user program like this to your workbench (IOB):

```simion.workbench_program()
function segment.other_actions()
if ion_px_mm >= 50 then  -- or "ion_time_of_flight" for time in microseconds
ion_splat = 1
end
end
```

Note however that it’s likely that a time-step straddles the position x = 50 mm, in which case your particle will stop slightly after 50 mm (which might not be ok with you). A simple way to reduce that small skip is to just increase your trajectory quality (TQual) factor. Or you could write a `tstep_adjust` segment to ensure that time steps end exactly at the desired end positions. More simply, you can use the test plane library (see SIMION Example: test_plane) available in 8.1 (or 8.0 Check for Updates), which handles this all for you.

Yet another solution, which involves no user programming and avoids the time-step inaccuracies, is to add a “stopper electrode” to your workbench. To do so, create a small PA file that contains only electrode points (no empty space) and add it to your workbench at the location where you want to splat ions. The PA must have the highest priority number in the PA Instances list on the PAs tab on the View screen; otherwise your particles might never see it. This stopper electrode, thankfully, will not affect your fields in the rest of your system because fields between PA instances do not interact in SIMION, so your particles will (correctly) not observe any effects from the stopper electrode until they actually hit it.

Another possibility is to refine your PA and then use the Modify Find/Replace functions to convert some non-electrode points to electrode points. However, be careful never to refine the array again, even if SIMION prompts you, because this electrode will affect the fields around it if you re-refine. Therefore, this approach may be error prone. Setting `simion.pas pa.refinable` to `false` can suppress the prompts to refine though.

### Get/set kinetic energy of current particle¶

In SIMION 8.2 you can obtain or change the current particle’s kinetic energy in a user program via the `ion_ke` reserved variables. Internally these read/write to the velocity (`ion_vx_mm`, `ion_vy_mm`, `ion_vz_mm`).

In SIMION 8.0/8.1 (a similarly in 7.0), you can use these functions:

```local function get_ke()
local speed = math.sqrt(ion_vx_mm^2 + ion_vy_mm^2 + ion_vz_mm^2)
return speed_to_ke(speed, ion_mass)
end

local function set_ke(ke)
local speed = math.sqrt(ion_vx_mm^2 + ion_vy_mm^2 + ion_vz_mm^2)
assert(speed ~= 0)
local new_speed = ke_to_speed(ke, ion_mass)
local scale = new_speed/speed
ion_vx_mm = ion_vx_mm * scale
ion_vy_mm = ion_vy_mm * scale
ion_vz_mm = ion_vz_mm * scale
end
```

### Flying particles in the reverse direction¶

As mentioned in Chapter 5 of the SIMION User Manual: “Try to back-fly your ions from their termination points (using reversed directions and termination energies). The ability to duplicate ion trajectories bi-directionally (at least approximately) means that the modeling is reasonably self-consistent. However, neither test guarantees that the results are really correct.”

This backflying can be done with Lua code. The ion color can be used to give these different colors and avoid backflying the backflown particles.

```function segment.other_actions()
if ion_splat ~= 0 and ion_color == 1 then
ion_splat = 0
ion_vx_mm = - ion_vx_mm
ion_vy_mm = - ion_vy_mm
ion_vz_mm = - ion_vz_mm
ion_color = 3
end
end
```

## Time step control¶

See Time Step Size on adjustment of time step size (trajectory quality, T.Qual) in the trajectory calculation.

### Data Record periodically, with constant size time steps, possibly for FFT¶

Sometimes you need particle data at periodic intervals. This is especially critical for Fast Fourier Transform (FFT) analysis.

The easiest way is to enable Time Markers (on the View screen Particles tab) for a given period in microseconds and then in Data Recording set When to Record to All Markers. Time markers force time steps to end precisely on multiples of the given period and can also be used as a Data Recording event.

A similar thing can alternately be done with user programming by executing `print` statements at the appropriate times (`ion_time_of_flight`) in an `other_actions` segment or invoking `mark()` within `other_actions` to trigger Data Recording when the Data Recording `All Markers` event is selected.

If using the programming method, care is required to ensure that the data output occurs on the precise times. `other_actions` is only executed at the end of time-steps, so you want to ensure that time-steps end at the precise times in which you want to print data, such as by using a `tstep_adjust` segment (or very large T.Qual values). Note, however, that numerical error (e.g. 0.9999999999999 does not equal 1) may require introducing a small tolerance factor in time calculation in a `tstep_adjust` segment.

The following code ensures that time step boundaries fall on multiples of time `step_size`. For example, 0.001, 0.002, 0.0003, etc. microseconds. SIMION may still insert additional time-steps unless you set the trajectory quality (TQual) to zero:

```-- SIMION 8.0/8.1
simion.workbench_program()
adjustable step_size = 0.001  -- microseconds
ion_time_step = math.min(ion_time_step, step_size)
end
```

## Custom physics¶

### Custom forces¶

Normally SIMION applies forces on particles according to the Lorentz Force Law, but you can customize this as described on Lorentz Force Law.

### Randomly convert ion to electron¶

To convert one particle to another particle (e.g. ion is replaced with an electron), just change the charge and mass (`ion_charge` and `ion_mass` variables). You probably also want to change the velocity (`ion_vx_mm`, `ion_vy_mm`, `ion_vz_mm`) or KE (Get/set kinetic energy of current particle) too.

```simion.workbench_program()

function segment.other_actions()
if ion_px_mm > 10 and ion_px_mm < 20 then
local speed_ion = sqrt(ion_vx_mm^2 + ion_vy_mm^2 + ion_vz_mm^2)
local collision_prob = 1 - exp(- speed_ion * ion_time_step / effective_mean_free_path_mm)
-- Test for collision
if rand() < collision_prob then
ion_charge = -1
ion_mass = 0.00054857990946
-- TODO: also change ion_vx_mm, ion_vy_mm, ion_vz_mm appropriately.
end
end
end
```

The above example performs conversions at random times with a given region. Here, `collision probability` is derived from a mean-free-path, similar to that done in SIMION Example: collsion_hs1 (also Collision Model HS1).

You can kind-of simulate this type of thing without user programs too, if you just do a second Fly’m with the child particles. Note that Data Recording has a “.ION” option which saves particles to its output file in ION format that can be read back (optionally after manual editing) into the Particles screen for a subsequent simulation.

### Fragment particles¶

What if one particle fragments into multiple particles or a particle hitting a surface causes multiple secondary emission particles? Techniques other than the above are required.

In SIMION 8.0, the simplest method is to fly two identical AB++ particles and then at a certain point in time (`other_actions` segment in a user program), change one particle to A+ and the other to B+. You should change not just the mass and charge but also the velocities (momentum) of the particles. For example:

```simion.workbench_program()
function segment.other_actions()
if ion_px_mm > 10 and ion_mass == 100 then
if ion_number == 1 then
ion_mass = 99
ion_charge = 1
ion_vx_mm = -0.1
ion_vy_mm = 0
ion_vz_mm = 0
else
ion_mass = 1
ion_charge = 1
ion_vx_mm = 10
ion_vy_mm = 0
ion_vz_mm = 0
end
end
end
```

SIMION also has an SIMION Example: child_particles though is a little more complicated.

More recently, SIMION 8.2 also has a newer API to create particles during the Fly’m, which is a more flexible way of doing these things: `simion.experimental.add_particles()`.

### Secondary emission¶

Seconary emissions effects is more involved topic with multiple approaches possible for simulating it. Specifically, the surface normal might not be obtained from the jagged PA surface but rather via the nearby (smooth) potential contour, a GEM file defintion, STL file definition, perhaps even surface enhancement info, or just a simple analytical equation provided. There are also various statistics that could applied to secondary energies and angles. See SIMION Example: secondary and SIMION Example: child_particles for some starting points.

### Collisional effects¶

For things like ion-gas collision or other randomized collisions, see Ion-Gas Collisions.

## Multiple runs or flym segment¶

Performing multiple runs is often used for things like Optimization. Examples with multiple runs include SIMION Example: tune (e.g. `tune_series.iob/lua`). In SIMION 8.0, this was done via the `sim_rerun_flym` variable. In SIMION 8.1, there is a new `segment.flym()` segment that makes this nicer.

Here’s a few additional notes on the `flym` segment. When click the `Fly'm` button the `segment.flym` executes. If you didn’t define a `segment.flym`, the default segment.flym is essentially as follows, which does a does a run:

```function segment.flym() run() end
```

But you can have any code you want inside the `segment.flym`. This prints “hello” without doing any runs:

```function segment.flym() print "hello" end
```

This prints “hello” and does two runs:

```function segment.flym()
print "hello"
run() run()
end
```

## Trajectory Display¶

### Clearing trajectories¶

The View and Retain checkboxes on the Particles tab control whether calculated particle trajectories are displayed on the screen and/or retained on disk. On redrawing the screen, particles will be cleared on the screen if they are not retained on disk. View and retain options can also be controlled programmatically via the `sim_trajectory_image_control` reserved variable. Redrawing occurs on window resizing or by issuing a `simion.redraw_screen()` command. Having `Rerun` enabled (or `sim_rerun_flym` reserved variable) also effectively disables retaining, although doing multiple runs instead via the `segment.flym()` segment does not cause that.

### Displaying only some trajectories¶

This will cause only every 10th particle trajectory to be displayed.

```simion.workbench_program()
function segment.other_actions()
if ion_number % 10 == 0 then
sim_trajectory_image_control = 0
else
sim_trajectory_image_control = 3
end
end
function segment.terminate_run()
sim_trajectory_image_control = 1  -- keep viewing
end
```

It uses the `sim_trajectory_image_control` reserved variable to suppress viewing and retaining of most trajectories. This technique may be useful if you are computing statistics on a lot of trajectories but only care to display some of them. It may be equivalently written using the modulus `%` operator:

```simion.workbench_program()
function segment.other_actions()
sim_trajectory_image_control = (ion_number % 10 == 0) and 0 or 3
end
function segment.terminate_run()
sim_trajectory_image_control = 1
end
```

This example preserves only every 10th time-step (for all particles):

```simion.workbench_program()
local n = 0
function segment.other_actions()
n = (n + 1) % 10
sim_trajectory_image_control = (n == 0) and 0 or 3
end
function segment.terminate_run()
sim_trajectory_image_control = 1
end
```

Another option, without programs, is to just do another Fly’m with a different and smaller set of particle definitions intended only for visualization, not statistics.

### Marker Color¶

If you set the marker color to 15 (transparent/white), SIMION will use the particle color for the marker color.

### Particle trails¶

This interactively shows particle traces only for the previous 10 time-steps. The effect is similar to Dots most flying but the particle has a short tail (like a tadpole).

```simion.workbench_program()
function segment.initialize_run()
sim_trajectory_image_control = 1
end
local n = 0
function segment.other_actions()
n = n + 1
if n % 10 == 0 then
simion.sleep(1.0)
simion.redraw_screen()
end
end
```

Explanation: The `sim_trajectory_image_control = 1` turns off trajectory retaining so that `simion.redraw_screen()` results in previous particle trajectories disappearing. This occurs every 10 time-steps (using the modulus operator `%` to compute the remainder of the counter divided by 10). `simion.sleep(0.5)` adds a 0.5 second delay, which can be adjusted.

## Fields seen by particles¶

### Overlapping PA instances¶

A particle in SIMION might not see the fields from all potential array instances in the workbench. A particle only sees the electric fields from the electric potential array instance (if any) with the highest priority number among the electric PA instances containing the current particle position. Likewise, a particle only sees the magnetic fields from the magnetic potential array instance (if any) with the highest priority number among the magnetic PA instances containing the current particle position. The priority numbers are displayed in the PA Instances list on the PAs tab on the View screen. They can be changed with the L- and L+ buttons.

Sometimes, you might not want the fields in the highest priority PA instance to be effective within the entire volume of that PA instance. Perhaps you broke a system into to PA instances that overlap and the fields are not valid on the ends of both PA instances (see SIMION Example: bender_cut). One solution is to crop that PA instance (using the Modify screen Crop function after refining the PA, as shown in SIMION Example: bender_cut). An alternative and simpler solution is to just utilize a user program to apply fields from a different PA instance, like in the example below:

```-- Apply E field from different PA instance (requires SIMION 8.2EA-20170127)
simion.workbench_program()
local pa1 = simion.wb.instances
local pa2 = simion.wb.instances
local painst = (ion_px_mm < 100) and pa1 or pa2
local ex,ey,ez = painst:field(ion_px_mm, ion_py_mm, ion_pz_mm, ion_time_of_flight)
ion_dvoltsx_gu, ion_dvoltsy_gu, ion_dvoltsz_gu
= painst:wb_to_pa_orient(-ex,-ey,-ez)
end
```

8.1 versions do not support time-dependent fields in `painst:field` (using `ion_time_of_flight`) but otherwise can do this if `fast_adjust`/`efield_adjust` segments are not used:

```-- Apply E field from different PA instance (requires SIMION 8.1.2.30 or above)
simion.workbench_program()
local pa1 = simion.wb.instances
local pa2 = simion.wb.instances
local painst = (ion_px_mm < 100) and pa1 or pa2
local ex,ey,ez = painst:field(ion_px_mm, ion_py_mm, ion_pz_mm)
ion_dvoltsx_gu, ion_dvoltsy_gu, ion_dvoltsz_gu
= painst:wb_to_pa_orient(-ex,-ey,-ez)
end
```

The following code superimposes fields from two magnetic fields:

```-- Superimposing two B fields (requires SIMION 8.2EA-20170127)
simion.workbench_program()
local pa1 = simion.wb.instances
local pa2 = simion.wb.instances
local bx1,by1,bz1 = pa1:bfield(ion_px_mm, ion_py_mm, ion_pz_mm, ion_time_of_flight)
local bx2,by2,bz2 = pa2:bfield(ion_px_mm, ion_py_mm, ion_pz_mm, ion_time_of_flight)
local bx = (bx1 or 0) + (bx2 or 0)
local by = (by1 or 0) + (by2 or 0)
local bz = (bz1 or 0) + (bz2 or 0)
ion_bfieldx_gu, ion_bfieldy_gu, ion_bfieldz_gu = wb_to_pa_orient(bx,by,bz)
end
-- B field in gauss at current particle position in workbench
-- coordinates (mm) from both PA instances.  Sums them.
-- Rotates vector back to current PA instance coordinates and sets it.
```

8.1 versions do not support time-dependent fields in `painst:field` (using `ion_time_of_flight`), but the `field` can accept electrode voltages as parameters:

```simion.workbench_program()
local pa1 = simion.wb.instances
local pa2 = simion.wb.instances
local ex,ey,ez
if ion_px_mm < 100 then
ex,ey,ez = pa1:field(ion_px_mm, ion_py_mm, ion_pz_mm, {
 = 200,  = 400
})
ex,ey,ez = pa1:wb_to_pa_orient(-ex,-ey,-ez)
else
ex,ey,ez = pa2:field(ion_px_mm, ion_py_mm, ion_pz_mm, {
 = -100,  = 300 * sin(ion_time_of_flight)
})
ex,ey,ez = pa2:wb_to_pa_orient(-ex,-ey,-ez)
end
ion_dvoltsx_gu, ion_dvoltsy_gu, ion_dvoltsz_gu = ex,ey,ez
end
```

Prior to 8.1.2.30, the `field` function was named `field_wc`.

## Output logging¶

Currently, everything written to the Log window (either via Data Recording GUI or printing in Lua) also gets written to the data recording output file (specified on the Data Recording screen).

Now, you could open your own file handle in Lua and write to it, and that would only write to that file (not the Log window):

```simion.workbench_program()
local fh = assert(io.open("output.txt", "w")
function segment.other_actions()
fh:write(ion_number, " ", ion_px_mm, "\n")
end
function segment.terminate_run()
fh:close()
end
```

This useful function ensures all the arguments are strings. It concatenates them with a tab and appends a new line. So, it works quite like `print` but accepts as a first argument a file handle:

```local function fprint(fh, ...)
local t = {n=select('#', ...), ...}
for i=1,t.n do t[i] = tostring(t[i]) end
local s = table.concat(t, '\t')
fh:write(s, "\n")
end
```

Example:

```fprint(fh, 123.456, "test", 789, math.pi, true)
```

Be sure to close the file (`fh:close()`) or at least flush (`fh:flush()`) when done writing.

Something that might be useful is the `simion.status()` function, which prints a message to the status bar at the bottom of the SIMION window:

```simion.workbench_program()
function segment.other_actions()
simion.status("The time is " .. ion_time_of_flight)
end
```

Another example displaying only periodically:

```simion.workbench_program()
local nions = 0
function segment.initialize() nions = ion_number end
local last_record_time = 0
local record_period = 0.5  -- microseconds
function segment.other_actions()
if math.abs(ion_time_of_flight - last_record_time) > record_period then
last_record_time = ion_time_of_flight
simion.status(("%0.2g usec, ion=%d/%d"):format(ion_time_of_flight, ion_number, nions))
end
end
```

You can also use the GUI Library (simion.experimental.dialog) library to display data in another window. See SIMION Example: gui (e.g. counter.lua in 2016-12-09 or above).

You can also write to Excel or gnuplot. See SIMION Example: excel, SIMION Example: gnuplot and SIMION Example: plot.