Simulation¶

PyExtal computes diffraction intensities using the Bloch wave method: for each incident beam direction (tilt) the structure-factor matrix is diagonalised, and intensities are assembled from the resulting eigenvalues and eigenvectors as a function of sample thickness.

Two types of simulation are supported:

CBED — convergent beam electron diffraction: a 2D image where each pixel corresponds to one incident angle inside a cone.
LARBED — large angle rocking beam: a stack of 2D maps, one per diffraction disk, where each pixel is the intensity in that disk at a given sample tilt.

Both start from a .dat configuration file. See Configuration for the full reference.

CBED simulation¶

Quick path — `simulate()`¶

simulate() runs the full pipeline (parse → eigendecompose → render) and returns a 2D intensity image alongside the SimParams object.

import numpy as np
import matplotlib.pyplot as plt
import pyextal.callBloch as cb

# displayParam: [x0, y0, gl, xs, ys, theta]
#   x0, y0  — pattern centre in pixels
#   gl      — geometric scale (pixel size in reciprocal space, Å⁻¹/pixel)
#   xs, ys  — output image size in pixels
#   theta   — pattern rotation angle (degrees)
displayParam = [512, 512, 0.05, 1024, 1024, 0.0]
thickness = 300  # Å

image, param = cb.simulate(
    'path/to/crystal.dat',
    thickness,
    displayParam,
)

plt.imshow(image, cmap='inferno')
plt.colorbar()
plt.show()

simulate() uses subaper=0 (full square tilt grid) internally. Do not pass subaper=1 — see Aperture modes below.

Step-by-step path — `bloch()` then assemble manually¶

For more control (e.g., to inspect per-tilt intensities or use a non-standard aperture), call bloch() and compute intensities yourself:

import numpy as np
import pyextal.callBloch as cb

param = cb.bloch('path/to/crystal.dat')

# param.tilts is a list of Data objects, one per tilt
# For each tilt, compute intensity from eigenvalues (gamma) and eigenvectors:
#   psi[g] = Σ_k  C[g,k] · cinv[k] · exp(i π γ_k t / K_n)
#   I[g]   = |psi[g]|²

thickness = 300  # Å
K_n = param.bigK * np.cos(np.radians(param.alpha))

intensities = []
for data in param.tilts:
    phase = np.exp(1j * np.pi * data.gamma * thickness / K_n)
    psi   = data.eigenvector @ (data.cinv * phase)
    I     = np.abs(psi) ** 2
    intensities.append(I)   # shape: (nout,)

Each Data object in param.tilts contains:

Field	Shape	Description
`ndiag`	scalar	Number of strongly excited beams for this tilt
`gamma`	`(ndiag,)`	Complex eigenvalues
`eigenvector`	`(nout, ndiag)`	Rows are output beams; columns are Bloch modes
`cinv`	`(ndiag,)`	First column of C⁻¹ — initial excitation amplitudes
`ix`, `iy`	scalar	Grid position of this tilt

LARBED simulation¶

LARBED computes the intensity of each diffraction disk as a 2D map over the tilt grid. The two-step structure (eigendecompose once, evaluate at many thicknesses) makes it efficient: call bloch() once, then call LARBED() or LARBED_tilt() at any number of thicknesses without repeating the eigensolve.

`LARBED()` — render the full 2D map¶

import numpy as np
import matplotlib.pyplot as plt
import pyextal.callBloch as cb

param = cb.bloch('path/to/crystal.dat')

thickness = 400  # Å
disks = cb.LARBED(param, thickness)   # shape: (nout, H, W)

fig, axes = plt.subplots(1, param.nout, figsize=(4 * param.nout, 4))
for i, ax in enumerate(axes):
    ax.imshow(disks[i], cmap='inferno')
    h, k, l = param.hklout[:, i]
    ax.set_title(f'({h} {k} {l})')
plt.tight_layout()
plt.show()

LARBED() renders the tilt grid of the param it is given, from the eigen-solutions stored by bloch() / bloch_run().

`LARBED_tilt()` — per-disk intensity array¶

LARBED_tilt() returns a (nout, ntilt) array — the intensity of each output disk at each tilt. This is the function used during refinement because it works on any tilt list, including one passed to tilt_run(), and supports thickness arrays (one value per tilt).

# Single thickness
intensity = cb.LARBED_tilt(param, thickness=400, nout=param.nout)
# shape: (nout, ntilt)

# Per-tilt thickness (e.g., wedge-shaped sample)
thickness_map = np.linspace(300, 500, param.ntilt)
intensity = cb.LARBED_tilt(param, thickness=thickness_map, nout=param.nout)

Number of cores¶

The eigen-decomposition runs on a thread pool; ncores sets its size (ncores=0 uses all physical cores). bloch(), bloch_run(), tilt_run() and simulate() default to ncores=1:

param = cb.bloch('crystal.dat', ncores=8)

See Simulation Engine for how the eigensolver works and how it scales.

Adjusting structure factors¶

To modify structure factors before simulation — for example to test a refinement step — use updateUgh():

from pyextal.callBloch import bloch_parse, bloch_run, updateUgh

param = bloch_parse('crystal.dat')

# reflections: list of (h, k, l) tuples
# values: (nref, 4) array — columns are [|Ugh|, phase°, |U'gh|, phase°]
reflections = [(1, 1, 1), (2, 2, 0)]
values = np.array([
    [0.0312, 180.0, 0.0015, 0.0],   # (1 1 1)
    [0.0210,   0.0, 0.0010, 0.0],   # (2 2 0)
])
updateUgh(reflections, values)

param = bloch_run(param)

updateUgh() must be called after bloch_parse() and before bloch_run(). Calling bloch_parse() again resets structure factors to the values in the .dat file.

Aperture modes¶

The subaper parameter controls which tilts are included in the tilt grid:

`subaper`	Description
`0`	Full square grid of side `2·nx+1`. Required by `simulate()`.
`1`	Circular aperture of radius `nx` (standard CBED geometry).
`2`	Custom circular aperture; specify radius via `subnx`.
`3`	Custom rectangular aperture; specify bounds via `subnx` and `subny`.

# Circular aperture
param = cb.bloch('crystal.dat', subaper=1)

# Custom rectangle: x in [-20, 30], y in [-10, 15]
param = cb.bloch(
    'crystal.dat', subaper=3,
    subnx=np.array([-20, 30], dtype=np.int32),
    subny=np.array([-10, 15], dtype=np.int32),
)

simulate() calls cbedp internally, which assumes a full square grid. Passing subaper=1 to simulate() will produce a shape mismatch and a crash. Use bloch() + manual rendering for non-square apertures.

Sample tilt¶

A sample tilt away from the zone axis can be set via t_x and t_y (degrees):

param = cb.bloch('crystal.dat', t_x=0.5, t_y=-0.3)

These shift the centre of the tilt grid relative to the zone axis direction defined in the .dat file.

Complete LARBED example¶

import numpy as np
import matplotlib.pyplot as plt
import pyextal.callBloch as cb

# 1. Eigendecompose once
param = cb.bloch('path/to/crystal.dat')
print(f'{param.ntilt} tilts, {param.nout} output disks')

# 2. Render disk maps at several thicknesses
thicknesses = np.arange(100, 600, 50)   # Å

fig, axes = plt.subplots(
    len(thicknesses), param.nout,
    figsize=(3 * param.nout, 3 * len(thicknesses)),
)

for row, t in enumerate(thicknesses):
    disks = cb.LARBED(param, thickness=t)   # shape: (nout, H, W)
    for col in range(param.nout):
        ax = axes[row, col] if param.nout > 1 else axes[row]
        ax.imshow(disks[col], cmap='inferno', origin='lower')
        if row == 0:
            h, k, l = param.hklout[:, col]
            ax.set_title(f'({h} {k} {l})')
        ax.set_ylabel(f't = {t} Å' if col == 0 else '')
        ax.axis('off')

plt.tight_layout()
plt.show()