Beam markers

Download notebook View on GitHubRun when the site was built.

Describe a neutral beam as markers: streams of fast atoms with a position, a velocity and a weight in atoms per second. Aim the beam with tangential_injection, give it a width, a divergence and a halo, and combine energy components.

A NeutralBeam holds its markers in a plasmacoll ParticleArrays (positions and velocities (N, 3) in Cartesian coordinates with z vertical, weights (N,)) and the isotope’s mass number. Its power is sum(weight * energy).

import matplotlib.pyplot as plt
import numpy as np

import freeofcharge as foc

origin, direction = foc.tangential_injection(5.3, 12.0)
closest = origin[:2] - (origin[:2] @ direction[:2]) * direction[:2]
print("start", origin.round(2), "direction", direction.round(3))
assert abs(np.hypot(*closest) - 5.3) < 1e-12
assert abs(np.hypot(*origin[:2]) - 12.0) < 1e-12

start [  5.3  -10.77   0.  ] direction [-0.  1.  0.]

tangential_injection(tangency_radius, start_radius) returns the start point and direction of a horizontal axis that starts at R = start_radius and passes closest to the machine’s axis at R = tangency_radius: the tangency radius is the injector’s main aiming choice. sense=-1 injects against +phi (counter-injection), height and toroidal_angle move the axis.

equilibrium = foc.Equilibrium.solovev(6.2, 1.8, 1.7, 5.3, wall_psi_n=1.3)
r_in, r_out = (
    equilibrium.boundary.vertices[:, 0].min(),
    equilibrium.boundary.vertices[:, 0].max(),
)
wall_out = equilibrium.wall.vertices[:, 0].max()

phi = np.linspace(0, 2 * np.pi, 400)
fig, ax = plt.subplots(figsize=(7.5, 5.5), layout="constrained")
ax.fill(
    r_out * np.cos(phi),
    r_out * np.sin(phi),
    color="C1",
    alpha=0.25,
    label="plasma (midplane)",
)
ax.fill(r_in * np.cos(phi), r_in * np.sin(phi), color="white")
ax.plot(wall_out * np.cos(phi), wall_out * np.sin(phi), color="0.5", label="outer wall")
ax.plot(5.3 * np.cos(phi), 5.3 * np.sin(phi), ":", color="0.4", label="tangency radius")
for tangency in (4.5, 5.3, 6.0):
    start, axis = foc.tangential_injection(tangency, 12.0)
    path = start[:, None] + np.linspace(0, 20, 2) * axis[:, None]
    ax.plot(path[0], path[1], label=f"axis, tangent to R = {tangency} m")
ax.set_aspect("equal")
ax.set_xlabel("x [m]")
ax.set_ylabel("y [m]")
ax.legend(loc="upper left", bbox_to_anchor=(1.02, 1.0), fontsize=8)
ax.set_title("Top view")
plt.show()

Width, divergence and halo

NeutralBeam.gaussian spreads the markers’ start points about the axis with a standard deviation width [m], and their directions with a standard deviation divergence [rad], each a scalar or a (horizontal, vertical) pair. A share halo_fraction of the markers gets the larger halo_divergence instead: the core and halo of a beamlet, as BBNBI describes it. All markers carry the same weight, power / (energy e) / markers atoms per second.

def angles(beam, axis):
    """The horizontal and vertical angles of the markers' directions to the axis."""
    horizontal = np.cross(axis, [0.0, 0.0, 1.0])
    horizontal /= np.linalg.norm(horizontal)
    vertical = np.cross(horizontal, axis)
    d = np.asarray(beam.directions)
    along = d @ axis
    return np.arctan(d @ horizontal / along), np.arctan(d @ vertical / along)


core = foc.NeutralBeam.gaussian(
    origin,
    direction,
    1.0e6,
    16.5e6,
    width=0.3,
    divergence=(0.005, 0.01),
    markers=20000,
    seed=1,
)
horizontal, vertical = angles(core, direction)
print(
    f"core divergence: {horizontal.std() * 1e3:.2f} mrad horizontal, {vertical.std() * 1e3:.2f} mrad vertical"
)
print(
    f"power {core.power() / 1e6:.2f} MW, {core.size} markers of {np.asarray(core.markers.weights)[0]:.3e} atoms/s"
)
assert (
    abs(horizontal.std() / 0.005 - 1) < 0.03 and abs(vertical.std() / 0.01 - 1) < 0.03
)
assert abs(core.power() / 16.5e6 - 1) < 1e-12

with_halo = foc.NeutralBeam.gaussian(
    origin,
    direction,
    1.0e6,
    16.5e6,
    divergence=0.005,
    halo_fraction=0.15,
    halo_divergence=0.03,
    markers=20000,
    seed=1,
)
halo_horizontal, _ = angles(with_halo, direction)
far = np.mean(np.abs(halo_horizontal) > 0.02)  # 4 core standard deviations
print(f"{100 * far:.1f} % of the markers beyond 20 mrad")
assert 0.15 * 0.45 < far < 0.15 * 0.60  # the halo's share beyond 2/3 of its sigma

bins = np.linspace(-0.08, 0.08, 81)
fig, ax = plt.subplots(figsize=(6, 3.5), layout="constrained")
ax.hist(horizontal, bins=bins, histtype="step", label="core only")
ax.hist(halo_horizontal, bins=bins, histtype="step", label="core and 15 % halo")
ax.set_yscale("log")
ax.set_xlabel("horizontal angle to the axis [rad]")
ax.set_ylabel("markers")
ax.legend()
plt.show()

core divergence: 5.00 mrad horizontal, 9.93 mrad vertical
power 16.50 MW, 20000 markers of 5.149e+15 atoms/s
8.0 % of the markers beyond 20 mrad

Energy components

A positive-ion source also extracts molecular ions, which leave atoms at a half and a third of the full energy. Build one beam per component and join them with NeutralBeam.concatenate (one isotope only). Here 70 %, 20 % and 10 % of the power are at the full, half and third energy:

full_energy, power = 120e3, 16.5e6
components = [
    foc.NeutralBeam.gaussian(
        origin,
        direction,
        full_energy / k,
        share * power,
        divergence=0.01,
        markers=2000,
        seed=k,
    )
    for k, share in ((1, 0.7), (2, 0.2), (3, 0.1))
]
beam = foc.NeutralBeam.concatenate(components)
energies = np.asarray(beam.energies)
for k, component in enumerate(components, start=1):
    print(
        f"E/{k}: {component.power() / 1e6:.2f} MW, {np.asarray(component.markers.weights).sum():.2e} atoms/s"
    )
assert abs(beam.power() / power - 1) < 1e-12
assert sorted(np.unique(energies.round(3))) == sorted(
    [full_energy, full_energy / 2, round(full_energy / 3, 3)]
)

E/1: 11.55 MW, 6.01e+20 atoms/s
E/2: 3.30 MW, 3.43e+20 atoms/s
E/3: 1.65 MW, 2.57e+20 atoms/s

The lower-energy components carry more atoms per watt and are stopped closer to the edge, as the next tutorial shows. With MPI, pass comm= to gaussian: each rank draws its own share of the markers and the deposition sums their profiles.