Python API

Signatures and docstrings of the public functions and classes in tratt at main (a4354eb), read from the source without importing it. The guides on this site say how the pieces fit together; this page is for looking up an argument.

tratt.core

Two-beam Takagi-Taupin for one frozen phase profile -> I(eta).

Everything is dimensionless:

- depth      ``zeta = z / L`` in [0, 1]
- coupling   ``alpha = L / Lambda_ext``  (Pendelloesung periods across the slab)
- rocking    ``eta``, the deviation parameter already multiplied by L

The common forward susceptibility ``chi_0`` is gauged out exactly: in symmetric
geometry it enters both equations with the same coefficient, so a shared phase
substitution removes it. Complex ``chi_0`` only leaves a common ``exp(-mu z)``,
which cancels in every normalized observable.

Absorption. In Laue a complex ``alpha`` carries the Borrmann (anomalous) part; in
Bragg normal absorption enters as ``eta -> eta + i kappa`` (see ``solve_bragg``).
A real ``alpha`` and ``kappa = 0`` are the non-absorbing case. See DESIGN.md for
the full assumption list.

Cost. Each RK4 step needs the phase factor ``e^{i phi}`` at the start, middle and
end of the step. The end of one step is the start of the next, so the factor is
evaluated twice per step, not four times.

Laue (transmission). Both beams travel +z; entrance ``D0(0)=1, Dh(0)=0``::

    dD0/dzeta = i pi alpha e^{+i phi} Dh
    dDh/dzeta = i pi alpha e^{-i phi} D0 + i eta Dh

Bragg (reflection). ``Dh`` travels -z, so ``d/ds_h = -d/dzeta``, and the boundary
conditions sit at opposite faces (``D0(0)=1``, ``Dh(1)=0``) -- a two-point BVP.
Rather than shoot, integrate the Riccati variable ``R = Dh / D0`` backward from
``R(1) = 0``, which turns it into an IVP in the stable direction::

    dR/dzeta = -i pi alpha e^{-i phi} - i eta R - i pi alpha e^{+i phi} R^2

Analytic ground truth, asserted in ``tests/test_core.py``:

- Laue, ``phi=0``, ``eta=0``  ->  ``|Dh(1)|^2 = sin^2(pi alpha)``
- Laue, real ``eta``          ->  ``|D0|^2 + |Dh|^2 = 1`` (lossless)
- Bragg, ``phi=0``            ->  total-reflection plateau for ``|eta| < 2 pi alpha``

PhaseField

class PhaseField(Protocol):
    ...

Frozen phase profile ``phi(zeta)``, vectorized over depth.

ZERO_PHASE

def ZERO_PHASE(zeta)

Unmodulated reference.

required_nz

def required_nz(
    alpha: float,
    alpha_ac: float,
    per_cycle: int = 20,
    S: float = 0.0,
    eta_max: float = 0.0,
) -> int

Depth steps needed to resolve the fastest oscillation in the integrand.

The grating turns over ``alpha_ac`` times across the slab and the
Pendelloesung solution ``alpha`` times. But the integrand also carries
``e^{-i phi}``, whose phase sweeps up to ``S q`` per unit depth, and
``e^{i eta zeta}``, which turns ``eta / 2 pi`` times; at large ``S`` and
high channels these dominate. Pass ``S`` and the largest ``|eta|`` to
include them: e.g. ``S = 6`` read at ``eta = 5 q`` needs ~4x the steps the
``alpha_ac`` rule gives (notebook Entry 31).

solve_laue

def solve_laue(
    alpha: complex,
    eta: ArrayLike = 0.0,
    phi: PhaseField | None = None,
    n_z: int = 512,
) -> tuple[NDArray[np.complex128], NDArray[np.complex128]]

Integrate the Laue system from ``zeta=0`` to ``1`` (fixed-step RK4).

Vectorized over ``eta``: the phase profile is shared, so a whole rocking
curve costs one sweep. Returns the exit amplitudes ``(D0, Dh)``.

A complex ``alpha`` models the Borrmann part of absorption in a
centrosymmetric crystal (``chi_h = chi_hbar``, so the same coupling enters
both equations). The common ``exp(-mu z)`` of ``chi_0''`` is still gauged out.

solve_bragg

def solve_bragg(
    alpha: float,
    eta: ArrayLike = 0.0,
    phi: PhaseField | None = None,
    n_z: int = 512,
    kappa: float = 0.0,
) -> NDArray[np.complex128]

Integrate the Bragg Riccati equation backward from ``R(1)=0`` to ``zeta=0``.

Returns the complex amplitude reflectivity ``R(0) = Dh(0)/D0(0)`` at the
entrance surface; the measured reflectivity is ``|R|^2``.

``kappa >= 0`` is normal absorption (chi_0''), which in Bragg does *not* gauge
out: the two beams travel in opposite directions, so it enters ``R`` as
``eta -> eta + i kappa``. A layer at depth ``zeta`` then contributes with weight
``exp(-kappa zeta)``, i.e. ``kappa`` is the round-trip amplitude attenuation per
unit depth (``mu L (1/gamma_0 + 1/gamma_h) / 2``). chi_h'' (anomalous absorption)
is not included. Thick slab, kappa > 0 -> the Darwin-Prins root of
``c R^2 + i (eta + i kappa) R + c = 0`` (tests/test_core.py).

amplitude

def amplitude(
    alpha: float,
    eta: ArrayLike,
    phi: PhaseField | None = None,
    geometry: str = 'laue',
    n_z: int = 512,
) -> NDArray[np.complex128]

Complex diffracted amplitude whose squared modulus is the intensity.

``geometry='laue'`` gives ``Dh(L)``; ``'bragg'`` gives ``R(0)``.

The phase matters for more than bookkeeping. Subtracting the unmodulated
solution at the *amplitude* level, ``Dh(S) - Dh(0)``, cancels the part of the
crystal the acoustic wave has not reached: in the kinematic limit the
integrand becomes ``(e^{-i phi} - 1)``, which is identically zero ahead of the
wavefront, so the difference is supported only on the filled depth. The same
subtraction done on intensities does not cancel anything.

rocking

def rocking(
    alpha: float,
    eta: ArrayLike,
    phi: PhaseField | None = None,
    geometry: str = 'laue',
    n_z: int = 512,
) -> NDArray[np.float64]

Diffracted intensity ``I(eta)`` for one frozen snapshot.

``geometry='laue'`` returns ``|Dh(L)|^2``; ``'bragg'`` returns ``|R(0)|^2``.

tratt.cube

One parameter point -> the time-resolved rocking curve ``I(eta, t)``.

Two cubes come out of every condition and they must share the eta grid, the
depth grid and the geometry, otherwise the residual between them mixes physics
with discretization error:

``front``  the real traveling front, ``phi = S A(t - zeta) sin(q zeta - q t)``
``full``   the same-phase steady reference, ``phi = S sin(q zeta - q t)``

Both are stroboscopic: the deformation is frozen at each t and the static
two-beam problem is solved in depth. There is no time derivative of the X-ray
field and no memory of earlier snapshots, so "formation" here means the approach
of the frozen-snapshot response to its periodic form, nothing more.

rocking_cube

def rocking_cube(
    alpha: float,
    drive: Drive,
    eta: NDArray[np.float64],
    times: NDArray[np.float64],
    geometry: str = 'laue',
    n_z: int | None = None,
    steady: bool = False,
) -> NDArray[np.float64]

``I(eta, t)`` as an ``(n_times, n_eta)`` array.

``steady=True`` swaps the traveling front for the same-phase full-profile
reference. Everything else is held fixed.

front_and_full

def front_and_full(
    alpha: float,
    drive: Drive,
    eta: NDArray[np.float64],
    times: NDArray[np.float64],
    geometry: str = 'laue',
    n_z: int | None = None,
) -> tuple[NDArray[np.float64], NDArray[np.float64]]

The two cubes a formation-time analysis needs, on identical grids.

satellite_cube

def satellite_cube(
    alpha: float,
    drive: Drive,
    eta: NDArray[np.float64],
    times: NDArray[np.float64],
    geometry: str = 'laue',
    n_z: int | None = None,
    steady: bool = False,
) -> NDArray[np.float64]

``|D_h(S) - D_h(0)|^2``: the diffraction the acoustic wave is responsible for.

A channel window is not a satellite. The crystal keeps its full thickness at
every instant, so the part the wavefront has not reached goes on producing
full-thickness fringes, and those land inside the satellite windows: at
alpha_ac = 8 the n = +-1 window is mostly them at early times. Measuring a
width there returns the window width or a fringe width, not the satellite.

Subtracting the unmodulated amplitude removes exactly that contribution. In
the kinematic limit the integrand of the difference carries a factor
``e^{-i phi} - 1``, which vanishes identically ahead of the wavefront, so the
difference is supported on the filled depth ``zeta < t`` alone and the
satellite it leaves has the finite-size width of that depth.

This is a constructed observable, not the measured intensity: an experiment
sees ``|D_h(S)|^2``. Report it as the acoustically induced amplitude
difference, and read widths off it rather than off the raw window.

tratt.drive

Traveling acoustic modulation -> frozen phase profile ``phi(zeta)``.

Slow time ``t = t_real / t_tr`` with wavefront transit ``t_tr = L / v_s``.
The drive of the abstract is::

    phi(zeta; t) = S * A(t - zeta) * sin(q*zeta - w*t)

with ``S = h.u0`` the modulation index, ``q = K*L``, ``w = Omega*t_tr``, and ``A``
a causal envelope of retarded time ``s = t - zeta`` (``A(s<0)=0`` reproduces the
wavefront cutoff ``z_max(t) = v_s t``).

**q and w are not independent.** For a non-dispersive traveling wave
``Omega = K*v_s``, so

    w = Omega*L/v_s = K*L = q = 2*pi*alpha_ac,     alpha_ac = L / Lambda_ac

leaving ``S`` and ``alpha_ac`` as the only drive parameters. With ``q = w`` the
phase field collapses to a pure function of retarded time,

    phi = -S * A(s) * sin(q*s),      s = t - zeta

i.e. for ``t >= 1`` the grating is a rigid profile sliding at ``v_s``, already in
its final form. That is why the formation time cannot exceed the transit time
under a hard step onset -- see DESIGN.md, fact 5.

The only knob that can push formation past the transit is a finite transducer
rise time ``rise = T_r / t_tr``, so it is an explicit parameter, not a smoothing
convenience.

Drive

@dataclass(frozen=True)
class Drive:
    S: float
    alpha_ac: float
    rise: float = 0.0

Non-dispersive traveling acoustic modulation in dimensionless groups.

Parameters
----------
S:
    Modulation index ``h . u0``.
alpha_ac:
    Acoustic thickness ``L / Lambda_ac``. Sets ``q = w = 2*pi*alpha_ac``, and
    hence the satellite spacing in ``eta``.
rise:
    Transducer rise time in units of ``t_tr``. ``0`` = hard step onset
    (strictly causal); ``>0`` = C1 smoothstep over ``s in [0, rise]``.

Drive.q

@property
def q(self) -> float

Carrier spatial phase across the slab, ``K*L``. Equals ``w``.

Drive.w

@property
def w(self) -> float

Acoustic phase advanced per wavefront transit, ``Omega*t_tr``.

Drive.period

@property
def period(self) -> float

Acoustic period in slow-time units: ``2*pi/w = 1/alpha_ac``.

Drive.envelope

def envelope(self, s: NDArray[np.float64]) -> NDArray[np.float64]

Causal envelope of retarded time ``s`` (units of ``t_tr``).

Drive.phi

def phi(self, t: float)

Frozen phase profile ``phi(zeta)`` at slow time ``t``.

Drive.steady_phi

def steady_phi(self, t: float)

Same-phase steady reference profile, the envelope removed.

``phi_full(zeta; t) = S sin(q zeta - q t)``: the diffraction the crystal
would show if the acoustic wave had *already* filled it, evaluated at the
same drive phase as ``phi(t)``. This is the baseline the formation time is
measured against.

Comparing against a period-*averaged* steady value is wrong: the steady
state itself oscillates, so that residual never settles. Comparing at
equal phase does settle -- and for ``rise = 0`` it settles exactly, since
``A(t - zeta) = 1`` for all ``zeta`` once ``t >= 1``, making the two
profiles identical. ``rise`` deliberately does not enter here.

Drive.eta_channel

def eta_channel(self, n: int) -> float

Center of diffraction channel ``n`` on the ``eta`` axis: ``eta_n = n*q``.

In the kinematic limit ``Dh(1) = i pi alpha int_0^1 e^{-i phi(zeta)}
e^{i eta (1-zeta)} dzeta``. Expanding
``e^{-i phi} = sum_n (-1)^n J_n(S) e^{i n (q zeta - q t)}`` leaves
``int_0^1 e^{i (n q - eta) zeta} dzeta``, which peaks at ``eta = n*q``.
Positions are analytic; nothing is peak-fitted.

``n`` is the harmonic index of the traveling grating. Labeling the
channels Stokes / anti-Stokes needs the frequency-resolved treatment --
the stroboscopic model carries no photon-energy axis, so that labeling is
deliberately not asserted here.

For integer ``alpha_ac`` the finite-slab sinc of one channel vanishes
exactly at every other channel center, which is why the Bessel ladder can
be read off the peak heights with no leakage correction.

Drive.times

def times(self, t_max: float=3.0, per_period: int=32) -> NDArray[np.float64]

Slow-time grid from 0 to ``t_max``, ``per_period`` samples per acoustic period.

tratt.floquet

Exact harmonic (Floquet) solution of Laue two-beam TT for a uniform travelling grating.

For phi(zeta) = S sin(q zeta - q t) over the whole slab (crystal fully filled), expand

    D0 = sum_k a_k(zeta) e^{+i k q zeta},      Dh = sum_m b_m(zeta) e^{-i m q zeta}.

With e^{-i phi} = sum_p J_p(S) e^{-i p (q zeta - q t)} the TT equations become a
*constant-coefficient* system, x' = i H x, x = (a_k, b_m):

    H[a_k, a_k] = -k q                 H[b_m, b_m] = eta + m q
    H[a_k, b_m] = pi alpha J_{k+m}(S) e^{-i (k+m) q t}
    H[b_m, a_k] = conj(H[a_k, b_m])    (H is Hermitian)

The expansion is redundant at zeta = 0; the choice a_k(0) = delta_k0, b_m(0) = 0 satisfies
D0(0) = 1, Dh(0) = 0, so by ODE uniqueness the series solves the original problem exactly.
Truncation to |k|, |m| <= K is the only approximation (convergence is tested).

Channel eta = n q is the resonance of b_{-n} with a_0 (diagonals 0 and eta - n q), with
direct coupling J_{-n}; every other a_k, b_m is off resonance by a multiple of q.
The eigenvalues of H are the Floquet propagation constants (dispersion-surface branches).

`mask(k, m) -> bool` removes individual couplings a_k <-> b_m (path-removal diagnostics).

auto_K

def auto_K(S: float, eta_max_orders: float=0.0, tol: float=1e-10, extra: int=3) -> int

Smallest K with sum_{|p| > K} J_p(S)^2 < tol, widened to cover the eta window and a margin.

FloquetLaue

@dataclass
class FloquetLaue:
    alpha: float
    alpha_ac: float
    S: float
    K: int
    t: float = 0.0
    mask: Callable | None = None

FloquetLaue.H

def H(self, eta)

Stack of H(eta), shape (n_eta, 2n, 2n).

FloquetLaue.propagate

def propagate(self, eta, zeta=1.0, chunk=2048)

Harmonic amplitudes x(zeta) for every eta: shape (n_eta, 2n).

FloquetLaue.fields

def fields(self, eta, zeta=1.0)

Physical exit amplitudes (D0, Dh) at depth zeta.

FloquetLaue.rocking

def rocking(self, eta)

FloquetLaue.eigen

def eigen(self, eta)

Floquet propagation constants and eigenvectors at each eta.

tratt.observe

``I(eta, t)`` -> channel observables -> formation times.

Channel ``n`` is the window ``W_n = {(n - 1/2) q < eta < (n + 1/2) q}`` around the
analytically known centre ``eta_n = n q``. Inside it,

    p_n(t)      = int_{W_n} I(eta, t) d eta            (integrated intensity)
    eta_c,n(t)  = int_{W_n} eta I d eta / p_n          (first moment)
    sigma_n(t)^2= int_{W_n} (eta - eta_c)^2 I d eta / p_n

``p_n`` is an integrated intensity, not a probability and not an occupancy: it is
never normalized by the sum over channels. ``eta_n`` is where the window is
anchored, not a claim about where the true maximum lies at strong coupling.

Formation time. With ``p_n^front`` the traveling-front result and ``p_n^full`` the
same-phase steady reference (see cube.front_and_full),

    r_n(t) = p_n^front(t) - p_n^full(t)
    a_n    = period average of p_n^full          (scale = window content)
           | <p_n^full> - p_n^static |           (scale = induced change)
    e_n(t) = |r_n(t)| / a_n

The second scale exists because a channel window is not a pure satellite: the
unmodulated rocking curve decays only as 1/eta^2, so its tail can dominate the
outer windows. It cancels in r_n but not in the normalizer.

and ``t_f`` is the earliest time such that ``e_n <= eps`` holds at that sample
*and at every later one*. Two cases are reported rather than faked:

- ``e_n`` never exceeds eps          -> t_f = times[0], "no formation detected"
- ``e_n`` still above eps at the end -> t_f = nan, "not converged in the window"

Returning the last sample in the second case would silently manufacture a
formation time equal to the observation length.

ChannelStats

@dataclass(frozen=True)
class ChannelStats:
    n: int
    p: NDArray[np.float64]
    eta_c: NDArray[np.float64]
    sigma: NDArray[np.float64]

Per-time observables of one channel. Arrays have the shape of ``times``.

Formation

@dataclass(frozen=True)
class Formation:
    t_f: float
    status: str
    eps: float
    a: float
    e: NDArray[np.float64]

Formation-time verdict for one channel, with the evidence attached.

Formation.ok

@property
def ok(self) -> bool

channel_window

def channel_window(eta: NDArray[np.float64], drive: Drive, n: int) -> NDArray[np.bool_]

Boolean mask of ``|eta - n q| < q / 2``.

channel_stats

def channel_stats(
    eta: NDArray[np.float64],
    I: NDArray[np.float64],
    drive: Drive,
    n: int,
) -> ChannelStats

Integrated intensity, centroid and width of channel ``n``.

``I`` is ``(n_times, n_eta)``. Where ``p_n`` underflows, the centroid and the
width are NaN rather than a ratio of two small numbers.

formation_time

def formation_time(
    times: NDArray[np.float64],
    p_front: NDArray[np.float64],
    p_full: NDArray[np.float64],
    period: float,
    eps: float = 0.05,
    a_min: float = 0.0,
    p_static: float | None = None,
) -> Formation

Earliest time after which ``|p_front - p_full| / a <= eps`` holds for good.

The residual ``p_front - p_full`` is clean: both cubes carry the same static
rocking-curve tail inside the window, so it cancels. The normalizer is what
needs care. The window content itself, ``<p_full>``, is not the satellite
strength -- at alpha = 1.3, alpha_ac = 5 the n = +-2 window holds more
unmodulated tail than satellite, which inflates ``a`` and biases ``t_f`` low.

Pass ``p_static``, the window integral of the S = 0 control, to normalize by
the acoustically *induced* change ``|<p_full> - p_static|`` instead. The
absolute value is required: modulation depletes the n = 0 window rather than
filling it, so the induced change is negative there.

Neither choice separates satellite from tail rigorously -- in a coherent
calculation the two interfere and do not add. ``a`` is a stated scale for the
residual, not a decomposition, so every reported ``t_f`` must say which one
it used.

``a_min`` rejects channels whose scale is too weak for the ratio to mean
anything; the default of 0 accepts everything and leaves that judgement to
the caller, which must then look at ``Formation.a``.

peak_fwhm

def peak_fwhm(
    eta: NDArray[np.float64],
    y: NDArray[np.float64],
    centre: float,
    half_width: float,
) -> tuple[float, float, float]

Full width at half maximum of the dominant peak inside one window.

Returns ``(fwhm, peak_height, peak_position)``; ``fwhm`` is NaN when the curve
does not fall to half maximum inside the window, which is the honest answer
for a feature wider than the window rather than a clipped number.

Measured on the highest point in the window, walking outward to the
half-maximum crossing and interpolating linearly. On a raw rocking curve this
picks up whatever fringe happens to be tallest; use it on the amplitude-
subtracted curve of ``cube.satellite_cube``, where the unmodulated crystal's
fringes are gone.

tratt.sweep

Parameter grid -> formation-time table.

One row per (alpha, S, channel). The raw cubes are not kept: a grid of them is
tens of GB, and every observable here is reproducible from the recorded metadata
by rerunning ``condition``.

Each condition needs three solves on one shared eta grid and one shared depth
grid: the traveling front, the same-phase steady reference, and the S = 0 static
control that sets the induced normalizer. Mixing grids between them would put
discretization error into the residual.

Grid

@dataclass(frozen=True)
class Grid:
    n_eta: int = 1201
    eta_span: float = 3.0
    t_max: float = 2.0
    per_period: int = 32
    n_z: int | None = None

Numerical resolution, recorded with every result so a row is reproducible.

condition

def condition(
    alpha: float,
    S: float,
    alpha_ac: float = 5.0,
    rise: float = 0.0,
    geometry: str = 'laue',
    channels: tuple[int, ...] = (0, 1, -1, 2, -2),
    eps: tuple[float, ...] = (0.02, 0.05, 0.1),
    grid: Grid = DEFAULT_GRID,
) -> list[dict]

Formation times for one parameter point, one row per channel.

sweep

def sweep(
    alphas: NDArray[np.float64],
    Ss: NDArray[np.float64],
    progress: bool = True,
    **kw,
) -> list[dict]

Cartesian (alpha, S) grid. Returns a flat list of rows.