Metadata-Version: 2.1
Name: treaaams
Version: 0.4.6rc6
Summary: Derivative of treams with complex frequency support
License: MIT
Requires-Python: >=3.10
Description-Content-Type: text/markdown
Requires-Dist: numpy >=1.23
Requires-Dist: scipy <1.17,>=1.9
Requires-Dist: tbump >=6.11.0
Requires-Dist: threadpoolctl >=3.6.0

# treaaams v0.4.6c6

A derivative of [treams](https://github.com/tfp-photonics/treams) with targeted
source-level patches to enable complex frequency evaluation.

## What treaaams changes in treams

treams uses the time convention $\mathrm{Re}(\boldsymbol E(\boldsymbol r, k_0)\,\mathrm e^{-\mathrm i c k_0 t})$
([treams paper][treams], Sec. 2.1), so resonances (poles of T- and S-matrices) lie
at $\mathrm{Im}\,k_0 < 0$. To reach them, every quantity has to be the **analytic
continuation** of its value at real $k_0$. treams instead picks branches pointwise,
e.g. "the root with $\mathrm{Im} \geq 0$", which is correct on the real axis but
jumps when $k_0$ crosses it. treaaams makes the changes below. At real $k_0$ the
physical results equal those of treams. The upstream test suite passes except for
two tests that evaluate $\Gamma(n, z)$ directly on its moved branch cut (see
[Tests](#tests)).

**1. Complex $k_0$ is accepted everywhere.** `TMatrix`, `TMatrixC`, `SMatrix` and
the `k0` annotation accept complex values. The basis `CylindricalWaveBasis` accepts
complex $k_z$. The Mie coefficients of spheres take complex size parameters
$x = k_0 r$, and those of cylinders take complex $k_0$ and $k_z$
(patches `020`–`024`, `050`).

**Branch cut convention.** Square roots, powers and $\Gamma(n, \cdot)$ are
multivalued. treams evaluates them on their principal branches, with the cut along
the negative real axis, and in places adjusts the sign afterwards (see 2.) ([DLMF 4.2.5](https://dlmf.nist.gov/4.2.E5),
[4.2.28](https://dlmf.nist.gov/4.2.E28)). treaaams places the cut of each of them
along a configurable ray $\arg\zeta = \varphi$, with $\varphi \in [-\pi, \pi]$. Of
the two branches bounded by that ray, it takes the one containing the positive
real axis: $\arg\zeta$ is taken in $(\varphi - 2\pi, \varphi]$ if $\varphi \geq 0$,
and in $(\varphi, \varphi + 2\pi]$ if $\varphi < 0$. Thus $\varphi = \pm\pi$ is the
principal branch, and all functions agree with treams for real positive arguments.
The settings are listed under [Configuration](#configuration).

**2. Square roots of the normal wave vector component are analytic.** Plane waves
use $k_z = \sqrt{\zeta}$ with $\zeta = k^2 - k_x^2 - k_y^2$, and cylindrical waves
use $k_\rho = \sqrt{k^2 - k_z^2}$. treams takes the root with $\mathrm{Im} \geq 0$
(see [treams paper][treams], text after Eq. (E.1), and [lattice paper][lattice],
text after Eq. (19)). treaaams uses the cut $\varphi =$ `BRANCH_CUT_SQRT_KZ`
$= -\pi/2$ instead (`misc.angled_sqrt`, patches `023`, `030`–`032`, `041`):

$$
\sqrt{\zeta} =
\begin{cases}
-\sqrt{\zeta}_{\mathrm{p}} & \arg \zeta \leq -\pi/2 \\
\phantom{-}\sqrt{\zeta}_{\mathrm{p}} & \text{otherwise}
\end{cases}
\qquad\Rightarrow\qquad \arg\sqrt{\zeta} \in \left(-\tfrac{\pi}{4}, \tfrac{3\pi}{4}\right],
$$

where $\sqrt{\zeta}_{\mathrm p}$ is the principal root
([DLMF 4.2.28](https://dlmf.nist.gov/4.2.E28)). For real $\zeta$ this gives the
same root as treams. Unlike treams, though, propagating orders ($\zeta$ close to
the positive real axis) and evanescent orders ($\zeta$ close to the negative real
axis) are continued smoothly into $\mathrm{Im}\,k_0 < 0$. There, the propagating
orders grow away from the scatterer, as quasi-normal modes do. All later sign
adjustments of $k_z$ and $k_\rho$ (in `sw`, `cw`, the field evaluations and
`cw.translate_periodic`) are removed.

**3. The incomplete gamma function of the Ewald lattice sums is continued.**

*Where it appears.* The reciprocal-space parts of the lattice sums contain the
upper incomplete gamma function $\Gamma(n, z) = \int_z^\infty t^{n-1}\mathrm e^{-t}\,\mathrm dt$
([DLMF 8.2.2](https://dlmf.nist.gov/8.2.E2)) at $z = -\gamma^2/(2\eta^2)$, with
$\gamma = \sqrt{1 - \beta^2}$, $\beta = |\boldsymbol k_\parallel + \boldsymbol G|/k$
and the Ewald parameter $\eta$:

- $\Gamma(1, z)$: [lattice paper][lattice], Eq. (20)
- $\Gamma(\tfrac12 - n, z)$: [lattice paper][lattice], Eqs. (25) and (35)
- $\Gamma(-n, z)$: [lattice paper][lattice], Eq. (30)
- origin terms $\Gamma(0, -\tfrac{1}{2\eta^2})$ and $\Gamma(-\tfrac12, -\tfrac{1}{2\eta^2})$:
  [lattice paper][lattice], Eq. (14)

*The problem.* For a propagating order at real $k_0$, $\beta < 1$, so $z$ is
negative real. That is exactly the cut of the principal $\Gamma(n, z)$, and the
lattice paper prescribes the value from below, $z - \mathrm i 0$ (text after
Eq. (14)). treams does the same. Once $k_0$ is complex, $z$ moves off the axis.
Below it, the principal value continues the value from below, but above it the
principal value belongs to the other side of the cut and jumps.

*The fix.* treaaams moves the cut out of the way, to the positive imaginary axis
($\varphi =$ `BRANCH_CUT_INCGAMMA` $= \pi/2$, patches `011`–`014`). Points between
the new and the old cut ($\pi/2 < \arg z \leq \pi$, second quadrant) are no longer
evaluated on the principal sheet. They are evaluated on the sheet reached by
continuing from below across the negative real axis. Following DLMF, that is
written $\Gamma(n, z\,\mathrm e^{-2\pi\mathrm i})$: the same point $z$, reached by
one clockwise turn around the origin. Everywhere else the principal value is
kept. The function with the original cut is still available as
`special.unshifted_incgamma`.

*The value on the neighbouring sheet* depends on the order $n$; with $m = -1$ in
this case:

| order $n$ | $\Gamma(n, z\,\mathrm e^{2\pi\mathrm i m})$ | why |
| --- | --- | --- |
| $1, 2, 3, \dots$ | $\Gamma(n, z)$ | the integrand $t^{n-1}\mathrm e^{-t}$ has no branch point, so there is only one sheet |
| $\pm\tfrac12, \pm\tfrac32, \dots$ | $2\,\Gamma(n) - \Gamma(n, z)$ | [DLMF 8.2.9](https://dlmf.nist.gov/8.2.E9) with $\mathrm e^{2\pi\mathrm i m n} = -1$ |
| $0, -1, -2, \dots$ | $\Gamma(n, z) - 2\pi\mathrm i\, m\, \dfrac{(-1)^{n}}{(-n)!}$ | derived below |

For the half-integer orders, the general continuation formula
[DLMF 8.2.9](https://dlmf.nist.gov/8.2.E9),

$$
\Gamma(n, z\,\mathrm e^{2\pi\mathrm i m}) = \mathrm e^{2\pi\mathrm i m n}\,\Gamma(n, z) + \left(1 - \mathrm e^{2\pi\mathrm i m n}\right)\Gamma(n),
$$

applies directly. For $n = -N$ with $N = 0, 1, 2, \dots$ it can't be used, since
$\Gamma(-N)$ is infinite while $1 - \mathrm e^{2\pi\mathrm i m n} = 0$. Instead, the
last row follows in three steps:

1. *Split off the multivalued part.* [DLMF 8.4.15](https://dlmf.nist.gov/8.4.E15)
   writes
   $$\Gamma(-N, z) = \frac{(-1)^N}{N!}\left(E_1(z) - \mathrm e^{-z}\sum_{k=0}^{N-1}\frac{(-1)^k k!}{z^{k+1}}\right).$$
   The sum contains only integer powers of $z$, which are the same on every
   sheet. All of the multivaluedness sits in the exponential integral
   $E_1(z) = \Gamma(0, z)$ ([DLMF 8.4.4](https://dlmf.nist.gov/8.4.E4)).
2. *Continue $E_1$.* $E_1$ has a logarithmic branch point at $0$, like
   $-\ln z$. Each turn around the origin therefore shifts it by a constant:
   $E_1(z\,\mathrm e^{2\pi\mathrm i m}) = E_1(z) - 2\pi\mathrm i m$
   ([DLMF 6.4.2](https://dlmf.nist.gov/6.4.E2)).
3. *Insert.* Only the $E_1$ term changes, so
   $$\Gamma(-N, z\,\mathrm e^{2\pi\mathrm i m}) = \Gamma(-N, z) - 2\pi\mathrm i\, m\,\frac{(-1)^N}{N!}.$$

As a check, for $N = 0$ and $m = -1$ the shift is $+2\pi\mathrm i$. That is
exactly the jump of $E_1$ across the negative real axis, from
$E_1(-x + \mathrm i 0)$ to $E_1(-x - \mathrm i 0)$
([DLMF 6.4.3](https://dlmf.nist.gov/6.4.E3)), which the moved cut removes.

*The accompanying power.* The sums use the reduced function
$\Gamma(n, z)/(-z)^n$ (`_esum._redincgamma`, patch `015`). Its power
$(-z)^n = \exp(n \operatorname{Ln}(-z))$ ([DLMF 4.2.26](https://dlmf.nist.gov/4.2.E26))
uses the cut $\arg(-z) =$ `BRANCH_CUT_CPOW` $= -\pi/2$, which is the same ray as
$\arg z = \pi/2$. The regularization near $z = 0$ can be configured through
`SINGULARITY_REDINCGAMMA` and `SINGULARITY_THRESH_REDINCGAMMA`.

**4. The direction of plane waves is given explicitly in the lattice expansions.**
Periodic arrays of spherical waves on a 2D lattice expand into plane waves with
the factor $\left(A k^2 \sqrt{1 - (\boldsymbol k_\parallel + \boldsymbol G)^2/k^2}\right)^{-1} = (A k k_z)^{-1}$
([treams paper][treams], Eq. (35)). Arrays of cylindrical waves expand with
$\propto (a k_y)^{-1}$ (Eq. (36)). The sign of the normal component follows the
side of the lattice, $k_z = \pm k\gamma$ for $z \gtrless 0$ (text after Eq. (E.1)).
treams applies that sign through a sign adjustment to $\mathrm{Im}\,k_z \geq 0$,
which is not analytic. In treaaams, `sw.periodic_to_pw` and `cw.periodic_to_pw` therefore
use $d\,k_z$ (or $d\,k_y$) with an explicit `direction` $d = \pm 1$
(patches `030`, `040`, `042`, `043`). The operators pass $d$ from the basis
($+1$ for `up` and $-1$ for `down`). For bases given by unit vectors, $d$ is the
sign of the unit vector's normal component. The keyword `direction` is optional.
Without it, $d$ is the treams sign of $k_z$ (or $k_y$), so results match treams,
but they are then not analytic.

**5. Mie coefficients.** The refractive index $n = \sqrt{\epsilon\mu}$ and the
relative impedance $Z = \sqrt{\mu/\epsilon}$ ([treams paper][treams], text after
Eq. (4)) take their roots with configurable cuts, `BRANCH_CUT_SQRT_MIE_N` and
`BRANCH_CUT_SQRT_MIE_Z` (default $-\pi$, i.e. the principal root as in treams).
The same root gives $k_\rho$ in the cylinder coefficients (patch `050`).

**Other changes.** The package is renamed to `treaaams`, so it can be installed
next to treams. BLAS is limited to one thread on import (`010`). The warning for
mismatched positions and T-matrices now shows both counts (`090`).

### Configuration

Every cut follows the [branch cut convention](#what-treaaams-changes-in-treams) and
can be read and set in `treaaams.config` with `get_<NAME>()` and
`set_<NAME>(value)`:

| name | default | cut along |
| --- | --- | --- |
| `BRANCH_CUT_SQRT_KZ` | $-\pi/2$ | $\arg(k^2 - k_\parallel^2) =$ `value` |
| `BRANCH_CUT_INCGAMMA` | $\pi/2$ | $\arg z =$ `value` |
| `BRANCH_CUT_CPOW` | $-\pi/2$ | $\arg(-z) =$ `value` |
| `BRANCH_CUT_SQRT_MIE_N`, `BRANCH_CUT_SQRT_MIE_Z` | $-\pi$ | $\arg(\epsilon\mu) =$ `value`, $\arg(\mu/\epsilon) =$ `value` |

[treams]: https://arxiv.org/abs/2309.03182 "D. Beutel, I. Fernandez-Corbaton, C. Rockstuhl, treams – A T-matrix scattering code for nanophotonic computations, arXiv:2309.03182 (2023)"
[lattice]: https://arxiv.org/abs/2207.06930 "D. Beutel, I. Fernandez-Corbaton, C. Rockstuhl, Unified lattice sums accommodating multiple sublattices for solutions of the Helmholtz equation in two and three dimensions, Phys. Rev. A 107, 013508 (2023)"

References:

- **[treams paper]**: D. Beutel, I. Fernandez-Corbaton, C. Rockstuhl, *treams – A
  T-matrix scattering code for nanophotonic computations*,
  [arXiv:2309.03182](https://arxiv.org/abs/2309.03182) (2023). Equation numbers
  refer to v1.
- **[lattice paper]**: D. Beutel, I. Fernandez-Corbaton, C. Rockstuhl, *Unified
  lattice sums accommodating multiple sublattices for solutions of the Helmholtz
  equation in two and three dimensions*, Phys. Rev. A **107**, 013508 (2023),
  [doi:10.1103/PhysRevA.107.013508](https://doi.org/10.1103/PhysRevA.107.013508),
  [arXiv:2207.06930](https://arxiv.org/abs/2207.06930). Equation numbers refer to
  the arXiv version.
- **DLMF**: [NIST Digital Library of Mathematical Functions](https://dlmf.nist.gov),
  §§ [4.2](https://dlmf.nist.gov/4.2), [6.4](https://dlmf.nist.gov/6.4),
  [8.2](https://dlmf.nist.gov/8.2), [8.4](https://dlmf.nist.gov/8.4).

## How it works

```
GitHub (treams vX.Y.Z)
        │  tarball download
        ▼
_build/treams_src/          ← clean upstream source
        │  patch -p1 (each file in patches/ in order)
        ▼
_build/treams_src/          ← patched source
        │  build (and publish) wheel
        ▼
install treaaams via pip/uv
```

## Prerequisites

- Python ≥ 3.10
- [uv](https://github.com/astral-sh/uv) (recommended) **or** pip
- A C compiler + `patch` utility (`apt install patch` / `brew install gpatch`)
- Cython and NumPy (installed automatically by the build backend)

## Quick start: editable install

```bash
# 1. Create a virtual environment (uv)
uv sync

# 2. Fetch upstream, apply patches, and install in editable mode
make editable

# 3. Verify
python -c "import treaaams; print(treaaams.__version__)"
```

## Tests

```bash
make test           # tests in tests/
make test-upstream  # upstream treams unit tests run against treaaams
```

Upstream tests that intentionally differ on treaaams are marked as strict `xfail`
in `tests/upstream_xfail.py`, together with the reason. CI
(`.github/workflows/tests.yml`) runs both on every push to `main` and on every pull
request for Python 3.10–3.13, and publishing a tag only proceeds when they pass.

## Adding or modifying a patch

### Workflow

```bash
# 1. Make sure you have a clean source tree
python scripts/build_patched.py --no-fetch   # skips download, resets patches

# 2. Edit the file you want to change
$EDITOR _build/treaaams_src/src/treams/_special/sph.pyx

# 3. Capture the diff as a numbered patch file
python scripts/make_patch.py src/treaaams/_special/sph.pyx \
    --name 010-custom-sph

# 4. Inspect it
cat patches/010-custom-sph.patch

# 5. Dry-run all patches to verify
make patch-dry-run

# 6. Rebuild
make editable
```

### Patch naming convention

Patches are applied **lexicographically**, so number them explicitly:

| File | Applied |
|---|---|
| `000-fix-typo.patch` | first |
| `010-override-ufunc.patch` | second |
| `020-custom-lattice.patch` | third |

## Tracking a new upstream version

Run for example
```
uv run tbump 0.4.7c1
```
To apply the patches to v0.4.7 of treams

If a patch no longer applies cleanly, `patch` will create a `.rej` file with
the rejected hunks.  Edit the patch accordingly and re-run.
