# 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](configuration.md) 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.

```python
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](#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:

```python
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

```python
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).

```python
# 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`:

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

See [Simulation Engine](backends.md) 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()`:

```python
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`. |

```python
# 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):

```python
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

```python
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()
```