calc

tit.calc

Temporal-interference envelope and carrier-field metrics.

Fields are zero-to-peak amplitudes, (N, 3) arrays in V/m for zero-to-peak channel currents, one array per channel (electrodes Ek+ and Ek-). The envelope functions take fields = [E1, E2, ..., E2K] in channel order, paired positionally: channels 2k-1 and 2k form beating pair k (carriers f_k and f_k + df, beating at the common difference frequency df), all K beats phase-aligned (K = 1 is TI, K >= 2 is mTI). In the formulas E_ka and E_kb are pair k's two channel fields. The carrier metrics take the fields of any number >= 1 of channels, one carrier each.

A channel with no beating partner enters the envelope functions with a zero field as its partner, [E1, E2, Eu, np.zeros_like(Eu)]: it adds to the carrier power but not to the beat, so it lowers the envelope. :func:hf_peak and :func:hf_power take it as one more channel field in the list.

• :func:envelope -- envelope amplitude along a direction (TI_normal)
• :func:envelope_max -- its maximum over all directions (TI_max, mTI_max) or within a plane (tangential maximum)
• :func:envelope_mean -- its mean over all directions (TI_avg)
• :func:hf_peak -- peak instantaneous carrier field
• :func:hf_power -- time-averaged mean-square carrier field, E_rms^2

The math and the accuracy table are in docs/wiki/analyzer.md. The K >= 2 sphere maximum runs in the numba kernel :mod:tit._calc_kernel; numba is a hard dependency (the image ships it).

References

Grossman et al. 2017, Cell 169(6):1029-1041, STAR Methods (K = 1 envelope). Botzanowski et al. 2025, Bioelectron Med 11:7, doi:10.1186/s42234-025-00169-6, pp. 4-5 (K >= 2 envelope, extending Grossman's: square, low-pass, square root). Hirata et al. 2024, Comput Biol Med 178:108697, Eq. 2 (K = 1 maximum). Cassarà et al. 2025, Bioelectromagnetics, Part I 46(2):e22542 and Part II 46(1):e22536 (carrier metrics); IEEE C95.1-2019 p. 35 (SAR).

The K >= 2 envelope model and the Fibonacci direction grid behind _hemisphere_directions are ported from Larissa Albantakis's branch alba/mTI_testing (the solvers here are new); the unpaired-channel convention follows Michael Sutherland's feature/uni_HF work.

envelope

envelope(fields: Sequence[ArrayLike], direction: ArrayLike) -> NDArray[float64]

Envelope modulation amplitude along a direction.

E_AM(n) = | |S n| - |D n| | = 4 |sum_k a_k b_k| / (|S n| + |D n|), with a_k = E_ka . n, b_k = E_kb . n and S, D the K x 3 matrices of rows E_ka + E_kb, E_ka - E_kb; the right-hand form (no cancellation) is evaluated. K = 1 gives 2 min(|a|, |b|).

Method: closed form, exact for every K. Literature: Grossman 2017's "envelope modulation amplitude along the direction n" (SimNIBS get_dirTI); along the cortical normal it is TI_normal.

Defined edge cases: a zero direction row and a null field (every projection 0) give 0.

Parameters:

Name	Type	Description	Default
fields	Sequence[ArrayLike]	[E1, E2, ..., E2K], one (N, 3) array per channel, in V/m.	required
direction	ArrayLike	(3,) or (N, 3); normalised internally, a zero row gives 0.	required

Returns:

Type	Description
NDArray[float64]	(N,) peak-to-trough envelope amplitude in V/m.

Raises:

Type	Description
ValueError	odd or fewer than two fields, shapes not all one (N, 3), NaN or inf, or direction not (3,) or (N, 3).

Source code in tit/calc.py

def envelope(fields: Sequence[ArrayLike], direction: ArrayLike) -> NDArray[np.float64]:
    """Envelope modulation amplitude along a direction.

    ``E_AM(n) = | |S n| - |D n| | = 4 |sum_k a_k b_k| / (|S n| + |D n|)``, with
    ``a_k = E_ka . n``, ``b_k = E_kb . n`` and ``S``, ``D`` the K x 3 matrices
    of rows ``E_ka + E_kb``, ``E_ka - E_kb``; the right-hand form (no
    cancellation) is evaluated.  K = 1 gives ``2 min(|a|, |b|)``.

    Method: closed form, exact for every K.  Literature: Grossman 2017's
    "envelope modulation amplitude along the direction n" (SimNIBS
    ``get_dirTI``); along the cortical normal it is ``TI_normal``.

    Defined edge cases: a zero ``direction`` row and a null field (every
    projection 0) give 0.

    Args:
        fields: ``[E1, E2, ..., E2K]``, one (N, 3) array per channel, in V/m.
        direction: (3,) or (N, 3); normalised internally, a zero row gives 0.

    Returns:
        (N,) peak-to-trough envelope amplitude in V/m.

    Raises:
        ValueError: odd or fewer than two fields, shapes not all one (N, 3),
            NaN or inf, or ``direction`` not (3,) or (N, 3).
    """
    E = _stack(fields, paired=True)
    n = _unit_rows(direction, E.shape[1], "direction")
    p = np.einsum("fnx,nx->fn", E, n)
    a, b = p[0::2], p[1::2]
    # w = u - v = 2 sum a b straight from the projections: _uvw's (s + d)(s - d) would lose
    # the weaker field of a pair to cancellation when the other dominates.
    u = 0.5 * np.sum((a + b) ** 2, axis=0)
    v = 0.5 * np.sum((a - b) ** 2, axis=0)
    return _env_uvw(u, v, 2.0 * np.sum(a * b, axis=0))

envelope_max

envelope_max(fields: Sequence[ArrayLike], *, perpendicular_to: ArrayLike | None = None, return_direction: bool = False) -> NDArray[float64] | tuple[NDArray[float64], NDArray[float64]]

Maximum of the envelope over all directions, or within a plane.

max_n E_AM(n) over unit n (:func:envelope), or over the unit n perpendicular to perpendicular_to, which equals the sphere maximum of the fields projected onto that plane.

Method: K = 1, sphere or plane, is an exact closed form (SimNIBS get_maxTI to ~1e-15; 0 on null fields). K >= 2 is exact in a plane (polynomial roots); on the sphere, an eigen-branch root-find plus Newton polish, never measured more than 1e-6 relative below a dense direction search. Literature: Grossman 2017's "maximal modulation amplitude along any orientation" (SimNIBS TImax), K = 1 in the form of Hirata 2024 Eq. 2. Outputs TI_max/mTI_max; the plane form is the tangential maximum.

Parameters:

Name	Type	Description	Default
fields	Sequence[ArrayLike]	[E1, E2, ..., E2K], one (N, 3) array per channel, in V/m.	required
perpendicular_to	ArrayLike | None	plane normal(s), (3,) or (N, 3); None maximises over the sphere. Normalised internally; a zero row gives 0.	None
return_direction	bool	also return the maximising unit direction.	False

Returns:

Type	Description
NDArray[float64] | tuple[NDArray[float64], NDArray[float64]]	amp (N,) in V/m, or (amp, n_hat) with n_hat (N, 3) unit.
NDArray[float64] | tuple[NDArray[float64], NDArray[float64]]	-n_hat is equally a maximiser; the sign gives n_hat . E_L >= 0
NDArray[float64] | tuple[NDArray[float64], NDArray[float64]]	for E_L the largest (projected) field, so at K = 1 amp * n_hat
NDArray[float64] | tuple[NDArray[float64], NDArray[float64]]	is the Grossman TI vector. Where amp == 0, n_hat is (1, 0, 0),
NDArray[float64] | tuple[NDArray[float64], NDArray[float64]]	or an in-plane unit vector for a plane.

Raises:

Type	Description
ValueError	the :func:envelope errors for fields, or perpendicular_to not (3,) or (N, 3) or containing NaN or inf.
RuntimeError	K >= 2 on the sphere, a 3x3 eigen-solve in the kernel did not converge (seen only for fields that overflow float64 once squared).

Source code in tit/calc.py

def envelope_max(
    fields: Sequence[ArrayLike],
    *,
    perpendicular_to: ArrayLike | None = None,
    return_direction: bool = False,
) -> NDArray[np.float64] | tuple[NDArray[np.float64], NDArray[np.float64]]:
    """Maximum of the envelope over all directions, or within a plane.

    ``max_n E_AM(n)`` over unit ``n`` (:func:`envelope`), or over the unit ``n``
    perpendicular to ``perpendicular_to``, which equals the sphere maximum of
    the fields projected onto that plane.

    Method: K = 1, sphere or plane, is an exact closed form (SimNIBS
    ``get_maxTI`` to ~1e-15; 0 on null fields).  K >= 2 is exact in a plane
    (polynomial roots); on the sphere, an eigen-branch root-find plus Newton
    polish, never measured more than 1e-6 relative below a dense direction
    search.  Literature: Grossman 2017's "maximal modulation
    amplitude along any orientation" (SimNIBS ``TImax``), K = 1 in the form of
    Hirata 2024 Eq. 2.  Outputs ``TI_max``/``mTI_max``; the plane form is the
    tangential maximum.

    Args:
        fields: ``[E1, E2, ..., E2K]``, one (N, 3) array per channel, in V/m.
        perpendicular_to: plane normal(s), (3,) or (N, 3); ``None`` maximises
            over the sphere.  Normalised internally; a zero row gives 0.
        return_direction: also return the maximising unit direction.

    Returns:
        ``amp`` (N,) in V/m, or ``(amp, n_hat)`` with ``n_hat`` (N, 3) unit.
        ``-n_hat`` is equally a maximiser; the sign gives ``n_hat . E_L >= 0``
        for ``E_L`` the largest (projected) field, so at K = 1 ``amp * n_hat``
        is the Grossman TI vector.  Where ``amp == 0``, ``n_hat`` is (1, 0, 0),
        or an in-plane unit vector for a plane.

    Raises:
        ValueError: the :func:`envelope` errors for ``fields``, or
            ``perpendicular_to`` not (3,) or (N, 3) or containing NaN or inf.
        RuntimeError: K >= 2 on the sphere, a 3x3 eigen-solve in the kernel did
            not converge (seen only for fields that overflow float64 once squared).
    """
    E = _stack(fields, paired=True)
    n_el = E.shape[1]
    plane = perpendicular_to is not None
    if plane:
        m = _unit_rows(perpendicular_to, n_el, "perpendicular_to")
        E = E - np.einsum("fnx,nx->fn", E, m)[..., None] * m
        E[:, ~m.any(axis=1)] = 0.0  # undefined plane -> 0, like a zero direction
    k = E.shape[0] // 2

    amp = np.empty(n_el)
    n_hat = np.empty((n_el, 3))
    for lo in range(0, n_el, _CHUNK):
        sl = slice(lo, min(lo + _CHUNK, n_el))
        if k == 1:
            amp[sl], n_hat[sl] = _max_k1(E[0, sl], E[1, sl])
        elif plane:
            amp[sl], n_hat[sl] = _max_plane(E[:, sl], m[sl])
        else:
            amp[sl], n_hat[sl] = _max_sphere(E[:, sl])
    if not return_direction:
        return amp

    if k >= 2:
        norms = np.linalg.norm(E, axis=2)
        e_big = E[norms.argmax(axis=0), np.arange(n_el)]
        n_hat[np.einsum("nx,nx->n", n_hat, e_big) < 0] *= -1.0
    null = amp <= 0.0
    if null.any():
        n_hat[null] = _frame(m[null])[1] if plane else (1.0, 0.0, 0.0)
    return amp, n_hat

envelope_mean

envelope_mean(fields: Sequence[ArrayLike]) -> NDArray[float64]

Mean of the envelope over all directions (TI_avg).

The mean of E_AM(n) over the unit sphere; at K = 1 exactly |E1| + |E2| - (|E1 + E2| + |E1 - E2|) / 2. Never exceeds :func:envelope_max.

Method: K = 1 closed form, exact. K >= 2: equal-area quadrature on 2048 hemisphere directions (E_AM(n) = E_AM(-n), so a 4096-point sphere rule), ~7e-4 relative error at worst. Literature: none, TI_avg is defined by this toolbox.

Parameters:

Name	Type	Description	Default
fields	Sequence[ArrayLike]	[E1, E2, ..., E2K], one (N, 3) array per channel, in V/m.	required

Returns:

Type	Description
NDArray[float64]	(N,) direction-averaged envelope amplitude in V/m.

Raises:

Type	Description
ValueError	as :func:envelope.

Source code in tit/calc.py

def envelope_mean(fields: Sequence[ArrayLike]) -> NDArray[np.float64]:
    """Mean of the envelope over all directions (``TI_avg``).

    The mean of ``E_AM(n)`` over the unit sphere; at K = 1 exactly
    ``|E1| + |E2| - (|E1 + E2| + |E1 - E2|) / 2``.  Never exceeds
    :func:`envelope_max`.

    Method: K = 1 closed form, exact.  K >= 2: equal-area quadrature on 2048
    hemisphere directions (``E_AM(n) = E_AM(-n)``, so a 4096-point sphere
    rule), ~7e-4 relative error at worst.  Literature: none, ``TI_avg`` is
    defined by this toolbox.

    Args:
        fields: ``[E1, E2, ..., E2K]``, one (N, 3) array per channel, in V/m.

    Returns:
        (N,) direction-averaged envelope amplitude in V/m.

    Raises:
        ValueError: as :func:`envelope`.
    """
    E = _stack(fields, paired=True)
    if E.shape[0] == 2:

        def norm(x):
            """Row norms."""
            return np.linalg.norm(x, axis=1)

        E1, E2 = E
        return norm(E1) + norm(E2) - 0.5 * (norm(E1 + E2) + norm(E1 - E2))

    dirs = _hemisphere_directions(2048)  # the ~7e-4 rule of the docstring
    x, y, z = dirs.T
    monomials = np.stack([x * x, y * y, z * z, 2 * x * y, 2 * x * z, 2 * y * z])

    def pack(C):
        """The 6 unique entries of each symmetric ``C``, in the order of ``monomials``."""
        return np.stack(
            [C[:, 0, 0], C[:, 1, 1], C[:, 2, 2], C[:, 0, 1], C[:, 0, 2], C[:, 1, 2]], 1
        )

    n_el = E.shape[1]
    out = np.empty(n_el)
    for lo in range(0, n_el, _MEAN_CHUNK):
        sl = slice(lo, min(lo + _MEAN_CHUNK, n_el))
        u, v, w = (pack(C) @ monomials for C in _forms(E[:, sl]))
        out[sl] = _env_uvw(u, v, w).mean(axis=1)
    return out

hf_peak

hf_peak(fields: Sequence[ArrayLike], *, coherent: bool = False) -> NDArray[float64]

Peak instantaneous carrier field, max_s |sum_c s_c E_c| over s_c = +-1.

Equal to max_n sum_c |E_c . n|; two channels give max(|E1 + E2|, |E1 - E2|). The default is the worst case over channel phases (reached by incommensurate frequencies), the safety-relevant value. coherent=True gives |sum_c E_c|, the phase-aligned sum some scripts call "peak HF".

Method: exact for any number of channels; closed form up to two, every sign pattern from three to eight (2^(F-1) of them), zonotope vertices from nine (cost ~F^3 per element for F fields, 3.5 ms at 52). Literature: Cassarà et al. 2025 Part I Eq. 3.

Parameters:

Name	Type	Description	Default
fields	Sequence[ArrayLike]	one or more channel fields (N, 3) in V/m, one per channel.	required
coherent	bool	carrier-only control (all channels at one frequency, in phase, ΔF = 0): the fields add as vectors (Cassarà 2025 Part II p. 8).	False

Returns:

Type	Description
NDArray[float64]	(N,) peak carrier field magnitude in V/m (zero-to-peak).

Raises:

Type	Description
ValueError	no fields, shapes not all one (N, 3), or NaN or inf.

Source code in tit/calc.py

def hf_peak(
    fields: Sequence[ArrayLike], *, coherent: bool = False
) -> NDArray[np.float64]:
    """Peak instantaneous carrier field, ``max_s |sum_c s_c E_c|`` over ``s_c = +-1``.

    Equal to ``max_n sum_c |E_c . n|``; two channels give ``max(|E1 + E2|,
    |E1 - E2|)``.  The default is the worst case over channel phases (reached
    by incommensurate frequencies), the safety-relevant value.  ``coherent=True``
    gives ``|sum_c E_c|``, the phase-aligned sum some scripts call "peak HF".

    Method: exact for any number of channels; closed form up to two, every
    sign pattern from three to eight (``2^(F-1)`` of them), zonotope vertices
    from nine (cost ~F^3 per element for F fields, 3.5 ms at 52).
    Literature: Cassarà et al. 2025 Part I Eq. 3.

    Args:
        fields: one or more channel fields (N, 3) in V/m, one per channel.
        coherent: carrier-only control (all channels at one frequency, in phase,
            ΔF = 0): the fields add as vectors (Cassarà 2025 Part II p. 8).

    Returns:
        (N,) peak carrier field magnitude in V/m (zero-to-peak).

    Raises:
        ValueError: no fields, shapes not all one (N, 3), or NaN or inf.
    """
    E = _stack(fields, paired=False)
    if coherent or len(E) == 1:
        return np.linalg.norm(E.sum(axis=0), axis=-1)
    if len(E) == 2:
        return np.maximum(
            np.linalg.norm(E[0] + E[1], axis=-1), np.linalg.norm(E[0] - E[1], axis=-1)
        )
    exact = _hf_peak_signs if len(E) <= _SIGN_ENUMERATION_MAX else _hf_peak_vertices
    n_el = E.shape[1]
    out = np.empty(n_el)
    for lo in range(0, n_el, _CHUNK):
        sl = slice(lo, min(lo + _CHUNK, n_el))
        out[sl] = exact(E[:, sl])
    return out

hf_power

hf_power(fields: Sequence[ArrayLike], *, coherent: bool = False) -> NDArray[float64]

Time-averaged mean-square carrier field, (1/2) sum_c |E_c|^2.

Channel fields at distinct frequencies add in power (their cross terms average out). The saved hf_rms field is sqrt(hf_power), and sigma * hf_power (isotropic sigma) the absorbed (Joule) power density hf_power_density in W/m^3. Unlike every other output this is a time average (an RMS squared), not a zero-to-peak amplitude.

Method: closed form, exact for any number of channels. Literature: Cassarà et al. 2025 Part I p. 11 and Part II Eq. 1; IEEE C95.1-2019 p. 35.

Parameters:

Name	Type	Description	Default
fields	Sequence[ArrayLike]	one or more zero-to-peak channel fields (N, 3) in V/m.	required
coherent	bool	carrier-only control (all channels at one frequency, in phase, ΔF = 0): one sinusoid of amplitude sum_c E_c, so (1/2) |sum_c E_c|^2.	False

Returns:

Type	Description
NDArray[float64]	(N,) mean-square carrier field in (V/m)^2.

Raises:

Type	Description
ValueError	as :func:hf_peak.

Source code in tit/calc.py

def hf_power(
    fields: Sequence[ArrayLike], *, coherent: bool = False
) -> NDArray[np.float64]:
    """Time-averaged mean-square carrier field, ``(1/2) sum_c |E_c|^2``.

    Channel fields at distinct frequencies add in power (their cross terms
    average out).  The saved ``hf_rms`` field is ``sqrt(hf_power)``, and
    ``sigma * hf_power`` (isotropic ``sigma``) the absorbed (Joule) power density
    ``hf_power_density`` in W/m^3.  Unlike every other output this is a time
    average (an RMS squared), not a zero-to-peak amplitude.

    Method: closed form, exact for any number of channels.  Literature:
    Cassarà et al. 2025 Part I p. 11 and Part II Eq. 1; IEEE C95.1-2019 p. 35.

    Args:
        fields: one or more zero-to-peak channel fields (N, 3) in V/m.
        coherent: carrier-only control (all channels at one frequency, in phase,
            ΔF = 0): one sinusoid of amplitude ``sum_c E_c``, so
            ``(1/2) |sum_c E_c|^2``.

    Returns:
        (N,) mean-square carrier field in (V/m)^2.

    Raises:
        ValueError: as :func:`hf_peak`.
    """
    E = _stack(fields, paired=False)
    if coherent:
        E = E.sum(axis=0, keepdims=True)
    return 0.5 * np.einsum("fnx,fnx->n", E, E)