diff --git a/CITATION.bib b/CITATION.bib
index 2b980d2..2843bd5 100644
--- a/CITATION.bib
+++ b/CITATION.bib
@@ -2,7 +2,7 @@ @misc{TransitionMatrices.jl
 	author  = {Gabriel Wu <wuzihua@pku.edu.cn> and contributors},
 	title   = {TransitionMatrices.jl},
 	url     = {https://github.com/JuliaRemoteSensing/TransitionMatrices.jl},
-	version = {v0.1.0},
-	year    = {2023},
-	month   = {2}
+	version = {v0.5.0},
+	year    = {2026},
+	month   = {6}
 }
diff --git a/README.md b/README.md
index 037cf94..fa239c0 100644
--- a/README.md
+++ b/README.md
@@ -4,7 +4,43 @@
 [![Build Status](https://github.com/JuliaRemoteSensing/TransitionMatrices.jl/actions/workflows/CI.yml/badge.svg?branch=main)](https://github.com/JuliaRemoteSensing/TransitionMatrices.jl/actions/workflows/CI.yml?query=branch%3Amain)
 [![Coverage](https://codecov.io/gh/JuliaRemoteSensing/TransitionMatrices.jl/branch/main/graph/badge.svg)](https://codecov.io/gh/JuliaRemoteSensing/TransitionMatrices.jl)
 
-The transition matrix method, or T-Matrix method, is one of the most powerful and widely used tools for rigorously computing electromagnetic scattering by single and compounded particles. As a package focusing on this method, `TransitionMatrices.jl` provides the following features:
+The transition matrix method, or **T-Matrix method**, is one of the most powerful and widely used tools for rigorously computing electromagnetic scattering by single and compounded particles. `TransitionMatrices.jl` is a generic, arbitrary-precision Julia implementation focused on this method.
+
+📖 **[Documentation](https://JuliaRemoteSensing.github.io/TransitionMatrices.jl/dev/)**
+
+## Installation
+
+`TransitionMatrices.jl` requires Julia ≥ 1.10. In the Julia REPL's package mode (press `]`):
+
+```julia-repl
+pkg> add TransitionMatrices
+```
+
+To track the development version, or if the package is not yet in your registry, add it by URL:
+
+```julia-repl
+pkg> add https://github.com/JuliaRemoteSensing/TransitionMatrices.jl
+```
+
+## Quick start
+
+```julia
+using TransitionMatrices
+
+# A prolate spheroid: semi-axes a = 1, c = 2, relative refractive index m = 1.5 + 0.01im
+s = Spheroid(1.0, 2.0, 1.5 + 0.01im)
+
+# Auto-converged (classic EBCM) T-matrix at wavelength λ = 2π
+𝐓 = transition_matrix(s, 2π)
+
+# Orientation-averaged far-field observables
+Qsca = scattering_cross_section(𝐓, 2π)
+Qext = extinction_cross_section(𝐓, 2π)
+ω    = albedo(𝐓)
+g    = asymmetry_parameter(𝐓, 2π)
+```
+
+## Features
 
 - Calculate the T-Matrix of various types of scatterers
   - Homogeneous spheres (via `bhmie`)
@@ -41,3 +77,127 @@ are re-exported by `TransitionMatrices.jl` and can be directly used after
 `using TransitionMatrices`.
 
 The `0.5` compatibility line uses `Quadmath.jl` 1.x and `Wigxjpf.jl` 0.3.x.
+
+## How to cite
+
+If you use `TransitionMatrices.jl` in your research, please cite it — a
+[`CITATION.bib`](CITATION.bib) is provided in this repository. Please also cite
+the original publication(s) for the specific method you use; see
+[Methods & references](#methods--references) below.
+
+## Methods & references
+
+Each solver and numerical trick follows the published literature below; the
+in-source docstrings and comments point to the specific equations used.
+
+**T-Matrix framework & far-field conventions**
+
+- P. C. Waterman, *Symmetry, unitarity, and geometry in electromagnetic
+  scattering*, [Phys. Rev. D **3**, 825–839 (1971)](https://doi.org/10.1103/PhysRevD.3.825)
+  — the null-field / EBCM origin of the T-Matrix method.
+- M. I. Mishchenko, L. D. Travis & A. A. Lacis, *Scattering, Absorption, and
+  Emission of Light by Small Particles*, Cambridge University Press (2002) — the
+  amplitude / phase / scattering matrices, cross sections, asymmetry parameter,
+  T-Matrix rotation, analytic random-orientation average, and the Wigner-d
+  recursions (the pervasive `Eq. (x.y)` references throughout the source).
+
+**Mie & coated spheres** (`bhmie`, `bhcoat`)
+
+- C. F. Bohren & D. R. Huffman, *Absorption and Scattering of Light by Small
+  Particles*, Wiley (1983) — the `bhmie` and `bhcoat` algorithms. The Mie
+  T-Matrix uses Mishchenko et al. (2002), Eqs. (5.42)–(5.44).
+
+**EBCM for axisymmetric shapes** (spheroids, cylinders, Chebyshev particles)
+
+- P. C. Waterman (1971), above.
+- M. I. Mishchenko & L. D. Travis, *Capabilities and limitations of a current
+  FORTRAN implementation of the T-matrix method for randomly oriented,
+  rotationally symmetric scatterers*, JQSRT **60**, 309–324 (1998) — the automatic
+  convergence procedure (`routine_mishchenko`: the choice of `nₘₐₓ` and the
+  Gauss division `Ng`), which the axisymmetric assembly is translated from.
+
+**Numerically stable EBCM — the F⁺ formulation** (`stable = true`)
+
+- W. R. C. Somerville, B. Auguié & E. C. Le Ru, *Severe loss of precision in
+  calculations of T-matrix integrals*, JQSRT **113**, 524 (2012) — the
+  catastrophic-cancellation diagnosis.
+- W. R. C. Somerville, B. Auguié & E. C. Le Ru, JQSRT **123**, 153 (2013),
+  [doi:10.1016/j.jqsrt.2012.07.017](https://doi.org/10.1016/j.jqsrt.2012.07.017)
+  — the cancellation-free `F⁺_{nk}(s,x)` projection (their Eq. 45–62, Table 2)
+  used for high-aspect-ratio spheroids.
+
+**Sh-matrix moment separation** (`prepare_sh`, fast `(λ, mᵣ)` sweeps)
+
+- The separation of the size / material parameters from the particle geometry is
+  the *parameter separation* of V. G. Farafonov, V. B. Il'in & M. S. Prokopjeva,
+  *Light scattering by multilayered nonspherical particles: a set of methods*,
+  JQSRT **79–80**, 599–626 (2003),
+  [doi:10.1016/S0022-4073(02)00310-2](https://doi.org/10.1016/S0022-4073(02)00310-2);
+  the "Sh-matrix" name and formalism are due to D. Petrov, Yu. Shkuratov, E. Zubko
+  & G. Videen, *Sh-matrices method as applied to scattering by particles with
+  layered structure*, JQSRT **106**, 437–454 (2007),
+  [doi:10.1016/j.jqsrt.2007.01.027](https://doi.org/10.1016/j.jqsrt.2007.01.027)
+  (extended in Petrov, Shkuratov & Videen, *The Sh-matrices method applied to
+  light scattering by small lenses*, JQSRT **110**, 1448–1459 (2009),
+  [doi:10.1016/j.jqsrt.2009.01.016](https://doi.org/10.1016/j.jqsrt.2009.01.016)).
+- The term-by-term radial integration reuses the `F⁺` projection of Somerville
+  et al. (2013) and the Riccati–Bessel power series of
+  [DLMF §10.53](https://dlmf.nist.gov/10.53).
+
+**Invariant Imbedding T-Matrix (IITM)** (axisymmetric, N-fold, and arbitrary shapes)
+
+- B. R. Johnson, *Invariant imbedding T matrix approach to electromagnetic
+  scattering*, [Appl. Opt. **27**, 4861–4873 (1988)](https://opg.optica.org/ao/abstract.cfm?uri=ao-27-23-4861)
+  — Eq. (97).
+- L. Bi, P. Yang, G. W. Kattawar & M. I. Mishchenko, *Efficient implementation of
+  the invariant imbedding T-matrix method and the separation of variables method
+  applied to large nonspherical inhomogeneous particles*, JQSRT **116**, 169–183
+  (2013) — Eq. (38).
+- A. Doicu & T. Wriedt, *The Invariant Imbedding T Matrix Approach*, in *The
+  Generalized Multipole Technique for Light Scattering* (Springer, 2018),
+  [doi:10.1007/978-3-319-74890-0_2](https://doi.org/10.1007/978-3-319-74890-0_2)
+  — Eq. (2.40).
+- B. Sun, L. Bi, P. Yang, M. Kahnert & G. Kattawar, *Invariant Imbedding T-matrix
+  Method for Light Scattering by Nonspherical and Inhomogeneous Particles*,
+  Elsevier (2019) — Eq. (4.2.36).
+- S. Hu (胡帅), *Research on the Numerical Computational Models and Application of
+  the Scattering Properties of Nonspherical Atmospheric Particles*, PhD
+  dissertation, National University of Defense Technology (2018) — Eq. (5.71).
+
+**Far-field observables & orientation averaging**
+
+- Mishchenko et al. (2002) — cross sections (Eqs. (5.102), (5.107), (5.140), (5.141)),
+  phase matrix (Eqs. (2.106)–(2.121)), asymmetry parameter (Eq. (4.92)), and the
+  analytic random-orientation T-Matrix (Eq. (5.96)).
+- L. Bi & P. Yang, *Accurate simulation of the optical properties of atmospheric
+  ice crystals with the invariant imbedding T-matrix method*, JQSRT **138**, 17–35
+  (2014), [doi:10.1016/j.jqsrt.2014.01.013](https://doi.org/10.1016/j.jqsrt.2014.01.013)
+  — the scattering-matrix expansion coefficients for a general T-Matrix (Eqs. (24)–(74)).
+
+**Shapes**
+
+- A. Mugnai & W. J. Wiscombe, *Scattering from nonspherical Chebyshev particles. 1*,
+  [Appl. Opt. **25**, 1235 (1986)](https://opg.optica.org/ao/abstract.cfm?uri=ao-25-7-1235)
+  — the Chebyshev-particle definition `r(ϑ) = r₀(1 + ε·Tₙ(cos ϑ))`.
+
+**Linearization & Jacobians**
+
+- R. Spurr, J. Wang, J. Zeng & M. I. Mishchenko, *Linearized T-matrix and Mie
+  scattering computations*, JQSRT (2012).
+- F. Xu & A. B. Davis, *Derivatives of light scattering properties of a
+  nonspherical particle computed with the T-matrix method*, Opt. Lett. **36**(22),
+  4464–4466 (2011), [doi:10.1364/OL.36.004464](https://doi.org/10.1364/OL.36.004464).
+- B. Sun, M. Gao, L. Bi & R. Spurr, *Analytical Jacobians of single scattering
+  optical properties using the invariant imbedding T-matrix method*, Opt. Express
+  **29**(6), 9635–9669 (2021), [doi:10.1364/OE.421886](https://doi.org/10.1364/OE.421886).
+- M. Gao & B. Sun, *Improvement and application of linearized invariant imbedding
+  T-matrix scattering method*, JQSRT **290**, 108322 (2022),
+  [doi:10.1016/j.jqsrt.2022.108322](https://doi.org/10.1016/j.jqsrt.2022.108322).
+
+**Numerical kernels (via dependencies)**
+
+- Wigner 3-j symbols via `Wigxjpf.jl`: H. T. Johansson & C. Forssén, *Fast and
+  accurate evaluation of Wigner 3j, 6j, and 9j symbols using prime factorization
+  and multiword integer arithmetic*, SIAM J. Sci. Comput. **38**(1), A376–A384 (2016).
+- Gauss–Legendre quadrature via
+  [`FastGaussQuadrature.jl`](https://github.com/JuliaApproximation/FastGaussQuadrature.jl).
diff --git a/docs/make.jl b/docs/make.jl
index 5a3100d..56833a6 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -74,10 +74,13 @@ makedocs(;
         assets = String[]),
     pages = [
         "Home" => "index.md",
+        "Theory & conventions" => "theory.md",
         "Usage" => "usage.md",
-        "Examples" => NB_PAGES,
+        "Examples" => vcat(["Overview" => "examples/index.md"], NB_PAGES),
         "Linearization" => "linearization.md",
-        "API" => "api.md"
+        "Performance" => "performance.md",
+        "API" => "api.md",
+        "Methods & references" => "references.md"
     ])
 
 if get(ENV, "DOCUMENTER_SKIP_DEPLOY", "false") != "true"
diff --git a/docs/src/api.md b/docs/src/api.md
index ee80457..b29cf5d 100644
--- a/docs/src/api.md
+++ b/docs/src/api.md
@@ -4,10 +4,75 @@ CurrentModule = TransitionMatrices
 
 # API
 
+The full public reference, grouped by topic. Internal (`_`-prefixed) helpers are
+omitted. Use the index below to jump to any symbol.
+
 ```@index
 ```
 
+## Shapes
+
+Scatterer geometries and their common interface.
+
+```@autodocs
+Modules = [TransitionMatrices]
+Pages = ["shapes/index.jl", "shapes/spheroid.jl", "shapes/cylinder.jl",
+    "shapes/chebyshev.jl", "shapes/prism.jl"]
+Filter = t -> !startswith(string(Base.nameof(t)), "_")
+```
+
+## Solvers & T-matrix construction
+
+The two-layer solver API ([`AbstractSolver`](@ref), the fixed solvers, and
+[`Iterative`](@ref)) and the underlying engines (Mie, EBCM with its `stable` and
+Sh-matrix variants, IITM).
+
+```@autodocs
+Modules = [TransitionMatrices]
+Pages = ["solvers.jl",
+    "EBCM/index.jl", "EBCM/axisymmetric.jl", "EBCM/routines.jl",
+    "EBCM/shmatrix.jl", "EBCM/stabilization.jl",
+    "IITM/index.jl", "IITM/axisymmetric.jl", "IITM/arbitrary.jl",
+    "IITM/nfold.jl", "IITM/fourier.jl",
+    "Mie/index.jl", "Mie/bhmie.jl", "Mie/bhcoat.jl", "Mie/MieTransitionMatrix.jl"]
+Filter = t -> !startswith(string(Base.nameof(t)), "_")
+```
+
+## Transition matrices, observables & post-processing
+
+The transition-matrix types and the far-field quantities derived from them
+(cross sections, amplitude/phase/scattering matrices, orientation averaging).
+
+```@autodocs
+Modules = [TransitionMatrices]
+Pages = ["common/index.jl", "common/AbstractTransitionMatrix.jl",
+    "common/AxisymmetricTransitionMatrix.jl",
+    "common/RandomOrientationTransitionMatrix.jl", "common/utils.jl"]
+Filter = t -> !startswith(string(Base.nameof(t)), "_")
+```
+
+## Linearization
+
+The differentiation framework and its backends (see
+[Linearization Framework](@ref) for the narrative).
+
+```@autodocs
+Modules = [TransitionMatrices]
+Pages = ["linearization.jl", "Mie/linearization.jl",
+    "EBCM/linearization.jl", "IITM/linearization.jl"]
+Filter = t -> !startswith(string(Base.nameof(t)), "_")
+```
+
+## Special functions
+
+The Riccati–Bessel, Wigner-d/D, Clebsch–Gordan, factorial, and quadrature
+kernels used internally (a few are exported for convenience).
+
 ```@autodocs
 Modules = [TransitionMatrices]
+Pages = ["special_functions/index.jl", "special_functions/bessel.jl",
+    "special_functions/clebschgordan.jl", "special_functions/factorial.jl",
+    "special_functions/quadrature.jl", "special_functions/wignerd.jl",
+    "compat/index.jl"]
 Filter = t -> !startswith(string(Base.nameof(t)), "_")
 ```
diff --git a/docs/src/examples/index.md b/docs/src/examples/index.md
new file mode 100644
index 0000000..2272f79
--- /dev/null
+++ b/docs/src/examples/index.md
@@ -0,0 +1,43 @@
+# Examples
+
+These pages are [Pluto.jl](https://plutojl.org) notebooks rendered to HTML (with
+plots). Each one is a self-contained walkthrough of one capability; the live,
+editable sources live in the [`examples/`](https://github.com/JuliaRemoteSensing/TransitionMatrices.jl/tree/main/examples)
+directory of the repository.
+
+## The notebooks
+
+- [**Shapes quick-start gallery**](shapes_gallery.md) — define a spheroid,
+  cylinder, Chebyshev particle, and N-fold prism; compute each T-matrix and its
+  `Qsca`/`Qext`/`g`/`ω`, with notes on which solver to use.
+- [**Solver landscape & convergence**](solver_landscape.md) — the two-layer API
+  (`EBCM`, `IITM`, `Iterative`, `stable`) on a high-aspect spheroid, showing
+  classic EBCM diverging with `nₘₐₓ` while the stabilized path holds.
+- [**Angular scattering & polarization**](angular_scattering.md) — the
+  orientation-averaged scattering (Mueller) matrix vs angle (phase function
+  `F₁₁`, linear polarization `-F₁₂/F₁₁`, …), plus fixed-geometry amplitude and
+  phase matrices.
+- [**Orientation averaging**](orientation_averaging.md) — Mishchenko's analytic
+  `RandomOrientationTransitionMatrix` vs the numerical `orientation_average`, and
+  the numerical average converging to it.
+- [**Spectral sensitivity**](spectral_sensitivity.md) — fast wavelength /
+  refractive-index sensitivities via the Sh-matrix moment-separation backend:
+  `prepare_sh` runs the geometry quadrature once, then `ForwardDiff` over the
+  cheap reconstruction gives exact `∂Cₛ𝚌ₐ/∂λ` and `∂Cₛ𝚌ₐ/∂mᵣ` across a spectrum.
+- [**Rain radar observables**](rain_radar.md) — a self-contained
+  radiative-transfer prototype: spheroidal IITM T-matrices, a rain-drop
+  axis-ratio model, water/ice refractive-index fits, and a drop-size distribution
+  feeding dual-pol radar moments and brightness-temperature integrals.
+
+## Running them locally
+
+The notebooks use the `examples/` environment, which develops the parent package:
+
+```julia-repl
+julia> import Pkg; Pkg.activate("examples"); Pkg.instantiate()  # first time only
+julia> using Pluto; Pluto.run()
+```
+
+Then open a notebook from the Pluto start page. Each notebook's first cell
+activates this environment and develops `TransitionMatrices` from the repository
+root, so it works from a fresh checkout.
diff --git a/docs/src/index.md b/docs/src/index.md
index b66fbe7..8c3dc77 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -1,10 +1,58 @@
+```@meta
+CurrentModule = TransitionMatrices
+```
+
 # TransitionMatrices.jl
 
 [![Dev](https://img.shields.io/badge/docs-dev-blue.svg)](https://JuliaRemoteSensing.github.io/TransitionMatrices.jl/dev/)
 [![Build Status](https://github.com/JuliaRemoteSensing/TransitionMatrices.jl/actions/workflows/CI.yml/badge.svg?branch=main)](https://github.com/JuliaRemoteSensing/TransitionMatrices.jl/actions/workflows/CI.yml?query=branch%3Amain)
 [![Coverage](https://codecov.io/gh/JuliaRemoteSensing/TransitionMatrices.jl/branch/main/graph/badge.svg)](https://codecov.io/gh/JuliaRemoteSensing/TransitionMatrices.jl)
 
-The transition matrix method, or T-Matrix method, is one of the most powerful and widely used tools for rigorously computing electromagnetic scattering by single and compounded particles. As a package focusing on this method, `TransitionMatrices.jl` provides the following features:
+The transition matrix method, or **T-Matrix method**, is one of the most powerful
+and widely used tools for rigorously computing electromagnetic scattering by
+single and compounded particles. `TransitionMatrices.jl` is a generic,
+arbitrary-precision Julia implementation focused on this method.
+
+## Installation
+
+`TransitionMatrices.jl` requires Julia ≥ 1.10. From the Julia REPL's package mode
+(press `]`):
+
+```julia-repl
+pkg> add TransitionMatrices
+```
+
+To track the development version, or if the package is not yet in your registry,
+add it by URL — this always works:
+
+```julia-repl
+pkg> add https://github.com/JuliaRemoteSensing/TransitionMatrices.jl
+```
+
+## Quick start
+
+```julia
+using TransitionMatrices
+
+# A prolate spheroid: semi-axes a = 1, c = 2, relative refractive index m = 1.5 + 0.01im
+s = Spheroid(1.0, 2.0, 1.5 + 0.01im)
+
+# Auto-converged (classic EBCM) T-matrix at wavelength λ = 2π
+𝐓 = transition_matrix(s, 2π)
+
+# Orientation-averaged far-field observables
+Qsca = scattering_cross_section(𝐓, 2π)   # scattering cross section
+Qext = extinction_cross_section(𝐓, 2π)   # extinction cross section
+ω    = albedo(𝐓)                         # single-scattering albedo
+g    = asymmetry_parameter(𝐓, 2π)        # asymmetry parameter ⟨cos Θ⟩
+```
+
+That is the whole loop: **define a shape → build its T-matrix → derive far-field
+quantities**. From here, [Usage](@ref) walks through every step, [Theory &
+conventions](@ref) defines the quantities and conventions used, and [Examples](@ref)
+collects runnable notebooks.
+
+## Features
 
 - Calculate the T-Matrix of various types of scatterers
   - Homogeneous spheres (via `bhmie`)
@@ -27,18 +75,40 @@ The transition matrix method, or T-Matrix method, is one of the most powerful an
   - Analytical fixed-geometry IITM material/wavelength slices for axisymmetric,
     n-fold, and arbitrary-shape solvers
 
-Compared to existing packages, `TransitionMatrices.jl` is special in that it is generic and supports various floating-point types, e.g.:
+## Floating-point genericity
+
+Compared to existing packages, `TransitionMatrices.jl` is special in that it is
+generic and supports various floating-point types, e.g.:
 
 - `Float64` and `BigFloat` from [`Base`](https://docs.julialang.org/en/v1/base/)
 - `Double64` from [`DoubleFloats.jl`](https://github.com/JuliaMath/DoubleFloats.jl)
 - `Float128` from [`Quadmath.jl`](https://github.com/JuliaMath/Quadmath.jl)
 - `Arb` and `Acb` from [`Arblib.jl`](https://github.com/kalmarek/Arblib.jl)
 
-For EBCM, by using higher-precision floating-point types, the maximum size
-parameter that can be handled is greatly improved.
+For EBCM, using a higher-precision floating-point type greatly improves the
+maximum size parameter that can be handled (see [Choosing a solver](@ref) for
+when to reach for extended precision versus the `stable` formulation).
 
 The precision types `Double64`, `Float128`, `ComplexF128`, `Arb`, and `Acb`
 are re-exported by `TransitionMatrices.jl` and can be directly used after
 `using TransitionMatrices`.
 
 The `0.5` compatibility line uses `Quadmath.jl` 1.x and `Wigxjpf.jl` 0.3.x.
+
+## Where to go next
+
+- [Theory & conventions](@ref) — the T-matrix definition, the index/sign
+  conventions, the size parameter, and how each observable is defined.
+- [Usage](@ref) — define a shape, choose a solver, build the T-matrix, and
+  compute far-field quantities, orientation averages, and derivatives.
+- [Examples](@ref) — runnable [Pluto.jl](https://plutojl.org) notebooks with plots.
+- [Linearization Framework](@ref) — analytical and automatic differentiation.
+- [API](@ref) — the full reference, grouped by topic.
+- [Methods & references](@ref) — the literature each method is built on.
+
+## How to cite
+
+If you use `TransitionMatrices.jl` in your research, please cite it — a
+[`CITATION.bib`](https://github.com/JuliaRemoteSensing/TransitionMatrices.jl/blob/main/CITATION.bib)
+is provided in the repository. Please also cite the original publication(s) for
+the specific method you use; see [Methods & references](@ref).
diff --git a/docs/src/linearization.md b/docs/src/linearization.md
index cb53054..d38e590 100644
--- a/docs/src/linearization.md
+++ b/docs/src/linearization.md
@@ -60,6 +60,17 @@ throw `UnsupportedLinearization`.
 This explicit failure is part of the design. It prevents users from mistaking a
 numerical finite-difference fallback for an analytical Jacobian.
 
+### Support at a glance
+
+| Backend | Model / shapes | Canonical variables | Notes |
+| --- | --- | --- | --- |
+| `MieLinearization()` | sphere (direct size-parameter model) | `:x`, `:mᵣ`, `:mᵢ`, `:λ` | `linearize_observable` covers `scattering_cross_section`, `extinction_cross_section`, `absorption_cross_section`, `albedo`, and fixed-angle `amplitude_matrix` |
+| `EBCMLinearization()` | spheroid, Chebyshev, cylinder (fixed `nₘₐₓ`, `Ng`) | spheroid `:a`, `:c`, `:mᵣ`, `:mᵢ`, `:λ`; Chebyshev `:r₀`, `:ε`, `:mᵣ`, `:mᵢ`, `:λ`; cylinder `:r`, `:h`, `:mᵣ`, `:mᵢ`, `:λ` | forward tangents propagated through `T = -P * inv(P + im * U)` |
+| `IITMLinearization()` (`:axisymmetric` / `:nfold` / `:arbitrary`) | fixed geometry | `:mᵣ`, `:mᵢ`, `:λ` | geometry variables not yet supported analytically |
+
+Always confirm a specific combination at runtime with
+`supports_linearization(problem, backend; output = :transition_matrix)`.
+
 The Mie backend is the first implemented backend. It supports the direct size
 parameter model with unique canonical variables drawn from `:x`, `:mᵣ`,
 `:mᵢ`, and `:λ`.
@@ -86,6 +97,21 @@ wavelength-dependent wavenumber, and Ricatti-Bessel functions, then propagate
 `ForwardDiff.jl` is still used in tests as the reference, but not in the EBCM
 production linearization path.
 
+A spheroid EBCM Jacobian is requested the same way as the Mie one — the rebuild
+function returns the canonical `shape`, `λ`, `nₘₐₓ`, and `Ng` fields:
+
+```julia
+problem = LinearizationProblem(
+    [2.0, 3.0, 1.311, 0.02, 2π];
+    variables = (:a, :c, :mᵣ, :mᵢ, :λ),
+) do x
+    (; shape = Spheroid(x[1], x[2], complex(x[3], x[4])), λ = x[5], nₘₐₓ = 10, Ng = 100)
+end
+
+result = linearize_transition_matrix(problem, EBCMLinearization())
+∂T_∂a = derivative(result, :a)   # ∂T / ∂a  (semi-major axis)
+```
+
 The IITM backend currently has fixed-geometry analytical slices for
 axisymmetric, n-fold, and arbitrary-shape solvers. `IITMLinearization()`,
 `IITMLinearization(:axisymmetric)`, `IITMLinearization(:nfold)`, and
@@ -136,8 +162,9 @@ dependence.
 
 - Spurr, Wang, Zeng, and Mishchenko, "Linearized T-matrix and Mie scattering
   computations", 2012.
-- Xu and Davis, "Analytical derivatives of T-matrix light scattering
-  quantities", 2011.
+- Feng Xu and Anthony B. Davis, "Derivatives of light scattering properties of a
+  nonspherical particle computed with the T-matrix method", Optics Letters 36(22),
+  4464--4466, 2011, doi:10.1364/OL.36.004464.
 - Sun, Gao, Bi, and Spurr, "Analytical Jacobians of single scattering optical
   properties using the invariant imbedding T-matrix method", Optics Express
   29(6), 9635--9669, 2021, doi:10.1364/OE.421886.
diff --git a/docs/src/performance.md b/docs/src/performance.md
new file mode 100644
index 0000000..2e7e123
--- /dev/null
+++ b/docs/src/performance.md
@@ -0,0 +1,61 @@
+# Performance & benchmarks
+
+## Performance notes
+
+A few rules of thumb when scattering calculations get expensive:
+
+- **Cost grows with the size parameter.** The truncation order needed for
+  convergence scales roughly with ``x = 2\pi r/\lambda``, and the dense linear
+  algebra is super-linear in that order. Large particles are intrinsically
+  costly.
+- **Reach for `stable` before raising precision on high-aspect spheroids.** For
+  an elongated spheroid in `Float64`, `Iterative(EBCM; stable = true)` is both
+  *more accurate and faster* than bumping the element type to `Double64` without
+  it; see [Usage](@ref) → *High-aspect-ratio spheroids*.
+- **Extended precision is several times slower.** `Double64`/`Arb` buy headroom
+  against round-off and larger size parameters, but at a real runtime cost. Use
+  them when accuracy demands it, not by default.
+- **Reuse a preparation for parameter sweeps.** When you need one shape at many
+  wavelengths or refractive indices, [`prepare_sh`](@ref) / [`ShMatrix`](@ref)
+  amortise the geometry quadrature across all points.
+- **Multithreading.** Several kernels use threads; start Julia with
+  `julia --threads=auto` (or set `JULIA_NUM_THREADS`) to use them.
+
+## The benchmark suite
+
+The repository ships a benchmark suite built on
+[BenchmarkTools.jl](https://github.com/JuliaCI/BenchmarkTools.jl), in
+[`benchmark/`](https://github.com/JuliaRemoteSensing/TransitionMatrices.jl/tree/main/benchmark).
+It follows the standard Julia convention — `benchmark/benchmarks.jl` defines a
+single global `const SUITE::BenchmarkGroup` — so it works unchanged with both
+[PkgBenchmark.jl](https://github.com/JuliaCI/PkgBenchmark.jl) and
+[AirspeedVelocity.jl](https://github.com/MilesCranmer/AirspeedVelocity.jl).
+
+Benchmark groups:
+
+- `special_functions` — Riccati–Bessel and Wigner-d recursions (inner kernels).
+- `ebcm` — EBCM blocks, fixed-order solves, and end-to-end auto-convergence.
+- `iitm` — invariant-imbedding T-matrix (axisymmetric and N-fold).
+- `postprocessing` — far-field observables from a precomputed T-matrix.
+- `linearization` — the analytical EBCM Jacobian.
+- `precision` — `Float64` vs `Double64` on the same EBCM block.
+
+### Running
+
+The simplest path uses the repository's `justfile`:
+
+```sh
+just bench            # instantiate the env and run the whole suite
+```
+
+Or directly, in the isolated `benchmark/` environment:
+
+```sh
+julia --project=benchmark -e 'using Pkg; Pkg.develop(path="."); Pkg.instantiate()'
+julia --project=benchmark benchmark/runbenchmarks.jl results-baseline.json
+```
+
+Absolute timings are machine-specific, so no numbers are reproduced here; the
+suite is meant for tracking *relative* changes (regressions and speed-ups)
+between revisions. See [`benchmark/README.md`](https://github.com/JuliaRemoteSensing/TransitionMatrices.jl/blob/main/benchmark/README.md)
+for the full layout.
diff --git a/docs/src/references.md b/docs/src/references.md
new file mode 100644
index 0000000..ff9f28e
--- /dev/null
+++ b/docs/src/references.md
@@ -0,0 +1,118 @@
+# Methods & references
+
+Each solver and numerical trick follows the published literature below; the
+in-source docstrings and comments point to the specific equations used. (This
+list mirrors the `Methods & references` section of the repository `README.md`,
+which is the authoritative copy.)
+
+**T-Matrix framework & far-field conventions**
+
+- P. C. Waterman, *Symmetry, unitarity, and geometry in electromagnetic
+  scattering*, [Phys. Rev. D **3**, 825–839 (1971)](https://doi.org/10.1103/PhysRevD.3.825)
+  — the null-field / EBCM origin of the T-Matrix method.
+- M. I. Mishchenko, L. D. Travis & A. A. Lacis, *Scattering, Absorption, and
+  Emission of Light by Small Particles*, Cambridge University Press (2002) — the
+  amplitude / phase / scattering matrices, cross sections, asymmetry parameter,
+  T-Matrix rotation, analytic random-orientation average, and the Wigner-d
+  recursions (the pervasive `Eq. (x.y)` references throughout the source).
+
+**Mie & coated spheres** (`bhmie`, `bhcoat`)
+
+- C. F. Bohren & D. R. Huffman, *Absorption and Scattering of Light by Small
+  Particles*, Wiley (1983) — the `bhmie` and `bhcoat` algorithms. The Mie
+  T-Matrix uses Mishchenko et al. (2002), Eqs. (5.42)–(5.44).
+
+**EBCM for axisymmetric shapes** (spheroids, cylinders, Chebyshev particles)
+
+- P. C. Waterman (1971), above.
+- M. I. Mishchenko & L. D. Travis, *Capabilities and limitations of a current
+  FORTRAN implementation of the T-matrix method for randomly oriented,
+  rotationally symmetric scatterers*, JQSRT **60**, 309–324 (1998) — the automatic
+  convergence procedure (`routine_mishchenko`: the choice of `nₘₐₓ` and the
+  Gauss division `Ng`), which the axisymmetric assembly is translated from.
+
+**Numerically stable EBCM — the F⁺ formulation** (`stable = true`)
+
+- W. R. C. Somerville, B. Auguié & E. C. Le Ru, *Severe loss of precision in
+  calculations of T-matrix integrals*, JQSRT **113**, 524 (2012) — the
+  catastrophic-cancellation diagnosis.
+- W. R. C. Somerville, B. Auguié & E. C. Le Ru, JQSRT **123**, 153 (2013),
+  [doi:10.1016/j.jqsrt.2012.07.017](https://doi.org/10.1016/j.jqsrt.2012.07.017)
+  — the cancellation-free `F⁺_{nk}(s,x)` projection (their Eq. 45–62, Table 2)
+  used for high-aspect-ratio spheroids.
+
+**Sh-matrix moment separation** (`prepare_sh`, fast `(λ, mᵣ)` sweeps)
+
+- The separation of the size / material parameters from the particle geometry is
+  the *parameter separation* of V. G. Farafonov, V. B. Il'in & M. S. Prokopjeva,
+  *Light scattering by multilayered nonspherical particles: a set of methods*,
+  JQSRT **79–80**, 599–626 (2003),
+  [doi:10.1016/S0022-4073(02)00310-2](https://doi.org/10.1016/S0022-4073(02)00310-2);
+  the "Sh-matrix" name and formalism are due to D. Petrov, Yu. Shkuratov, E. Zubko
+  & G. Videen, *Sh-matrices method as applied to scattering by particles with
+  layered structure*, JQSRT **106**, 437–454 (2007),
+  [doi:10.1016/j.jqsrt.2007.01.027](https://doi.org/10.1016/j.jqsrt.2007.01.027)
+  (extended in Petrov, Shkuratov & Videen, *The Sh-matrices method applied to
+  light scattering by small lenses*, JQSRT **110**, 1448–1459 (2009),
+  [doi:10.1016/j.jqsrt.2009.01.016](https://doi.org/10.1016/j.jqsrt.2009.01.016)).
+- The term-by-term radial integration reuses the `F⁺` projection of Somerville
+  et al. (2013) and the Riccati–Bessel power series of
+  [DLMF §10.53](https://dlmf.nist.gov/10.53).
+
+**Invariant Imbedding T-Matrix (IITM)** (axisymmetric, N-fold, and arbitrary shapes)
+
+- B. R. Johnson, *Invariant imbedding T matrix approach to electromagnetic
+  scattering*, [Appl. Opt. **27**, 4861–4873 (1988)](https://opg.optica.org/ao/abstract.cfm?uri=ao-27-23-4861)
+  — Eq. (97).
+- L. Bi, P. Yang, G. W. Kattawar & M. I. Mishchenko, *Efficient implementation of
+  the invariant imbedding T-matrix method and the separation of variables method
+  applied to large nonspherical inhomogeneous particles*, JQSRT **116**, 169–183
+  (2013) — Eq. (38).
+- A. Doicu & T. Wriedt, *The Invariant Imbedding T Matrix Approach*, in *The
+  Generalized Multipole Technique for Light Scattering* (Springer, 2018),
+  [doi:10.1007/978-3-319-74890-0_2](https://doi.org/10.1007/978-3-319-74890-0_2)
+  — Eq. (2.40).
+- B. Sun, L. Bi, P. Yang, M. Kahnert & G. Kattawar, *Invariant Imbedding T-matrix
+  Method for Light Scattering by Nonspherical and Inhomogeneous Particles*,
+  Elsevier (2019) — Eq. (4.2.36).
+- S. Hu (胡帅), *Research on the Numerical Computational Models and Application of
+  the Scattering Properties of Nonspherical Atmospheric Particles*, PhD
+  dissertation, National University of Defense Technology (2018) — Eq. (5.71).
+
+**Far-field observables & orientation averaging**
+
+- Mishchenko et al. (2002) — cross sections (Eqs. (5.102), (5.107), (5.140), (5.141)),
+  phase matrix (Eqs. (2.106)–(2.121)), asymmetry parameter (Eq. (4.92)), and the
+  analytic random-orientation T-Matrix (Eq. (5.96)).
+- L. Bi & P. Yang, *Accurate simulation of the optical properties of atmospheric
+  ice crystals with the invariant imbedding T-matrix method*, JQSRT **138**, 17–35
+  (2014), [doi:10.1016/j.jqsrt.2014.01.013](https://doi.org/10.1016/j.jqsrt.2014.01.013)
+  — the scattering-matrix expansion coefficients for a general T-Matrix (Eqs. (24)–(74)).
+
+**Shapes**
+
+- A. Mugnai & W. J. Wiscombe, *Scattering from nonspherical Chebyshev particles. 1*,
+  [Appl. Opt. **25**, 1235 (1986)](https://opg.optica.org/ao/abstract.cfm?uri=ao-25-7-1235)
+  — the Chebyshev-particle definition `r(ϑ) = r₀(1 + ε·Tₙ(cos ϑ))`.
+
+**Linearization & Jacobians**
+
+- R. Spurr, J. Wang, J. Zeng & M. I. Mishchenko, *Linearized T-matrix and Mie
+  scattering computations*, JQSRT (2012).
+- F. Xu & A. B. Davis, *Derivatives of light scattering properties of a
+  nonspherical particle computed with the T-matrix method*, Opt. Lett. **36**(22),
+  4464–4466 (2011), [doi:10.1364/OL.36.004464](https://doi.org/10.1364/OL.36.004464).
+- B. Sun, M. Gao, L. Bi & R. Spurr, *Analytical Jacobians of single scattering
+  optical properties using the invariant imbedding T-matrix method*, Opt. Express
+  **29**(6), 9635–9669 (2021), [doi:10.1364/OE.421886](https://doi.org/10.1364/OE.421886).
+- M. Gao & B. Sun, *Improvement and application of linearized invariant imbedding
+  T-matrix scattering method*, JQSRT **290**, 108322 (2022),
+  [doi:10.1016/j.jqsrt.2022.108322](https://doi.org/10.1016/j.jqsrt.2022.108322).
+
+**Numerical kernels (via dependencies)**
+
+- Wigner 3-j symbols via `Wigxjpf.jl`: H. T. Johansson & C. Forssén, *Fast and
+  accurate evaluation of Wigner 3j, 6j, and 9j symbols using prime factorization
+  and multiword integer arithmetic*, SIAM J. Sci. Comput. **38**(1), A376–A384 (2016).
+- Gauss–Legendre quadrature via
+  [`FastGaussQuadrature.jl`](https://github.com/JuliaApproximation/FastGaussQuadrature.jl).
diff --git a/docs/src/theory.md b/docs/src/theory.md
new file mode 100644
index 0000000..6d5d8d3
--- /dev/null
+++ b/docs/src/theory.md
@@ -0,0 +1,145 @@
+```@meta
+CurrentModule = TransitionMatrices
+```
+
+# Theory & conventions
+
+This page summarises the definitions and conventions used throughout
+`TransitionMatrices.jl`. They follow Mishchenko, Travis & Lacis, *Scattering,
+Absorption, and Emission of Light by Small Particles* (Cambridge University
+Press, 2002) — referred to below as **MTL** — which the in-source docstrings
+cite by equation number. Full references are listed in
+[Methods & references](@ref).
+
+## The transition matrix
+
+For a single scatterer embedded in a non-absorbing host medium, both the
+incident and the scattered electric fields are expanded in **vector spherical
+wave functions** (VSWFs). The incident field uses the *regular* VSWFs
+``\mathrm{Rg}\mathbf{M}_{mn}, \mathrm{Rg}\mathbf{N}_{mn}`` with expansion
+coefficients ``(a_{mn}, b_{mn})``; the scattered field uses the *outgoing*
+VSWFs ``\mathbf{M}_{mn}, \mathbf{N}_{mn}`` with coefficients ``(p_{mn}, q_{mn})``.
+Because Maxwell's equations are linear, the two sets of coefficients are linearly
+related, and the matrix of that relation is the **transition matrix** ``T`` (MTL §5):
+
+```math
+\begin{bmatrix} p_{mn} \\ q_{mn} \end{bmatrix}
+  = \sum_{n'=1}^{\infty} \sum_{m'=-n'}^{n'}
+    \begin{bmatrix} T^{11}_{mn m'n'} & T^{12}_{mn m'n'} \\
+                    T^{21}_{mn m'n'} & T^{22}_{mn m'n'} \end{bmatrix}
+    \begin{bmatrix} a_{m'n'} \\ b_{m'n'} \end{bmatrix}.
+```
+
+The transition matrix depends only on the particle (its shape, size, and
+refractive index relative to the host) and on the wavelength — **not** on the
+incident field. Computing it once therefore gives access to scattering for *any*
+incidence direction and polarisation, and to orientation averages in closed form.
+
+### Storage
+
+[`AbstractTransitionMatrix`](@ref) stores ``T`` as a six-dimensional array
+indexed ``T_{mn m'n'}^{kl}``: the unprimed ``(m, n)`` label the scattered
+harmonic, the primed ``(m', n')`` the incident one, and ``k, l \in \{1, 2\}``
+select the magnetic (``\mathbf{M}``) and electric (``\mathbf{N}``) parts. The
+degree runs ``1 \le n \le n_{\max}`` and the order ``-n \le m \le n``, so a
+truncation order ``n_{\max}`` fixes the matrix size. Concrete subtypes
+([`TransitionMatrix`](@ref), [`AxisymmetricTransitionMatrix`](@ref),
+[`MieTransitionMatrix`](@ref), [`RandomOrientationTransitionMatrix`](@ref))
+exploit symmetry to store and apply ``T`` more compactly.
+
+## Size parameter and units
+
+Lengths and the wavelength are expressed in the same (arbitrary) unit; only their
+ratio matters. The wavenumber in the host medium is ``k = 2\pi/\lambda``, and the
+dimensionless **size parameter** of a scatterer of characteristic radius ``r`` is
+
+```math
+x = k r = \frac{2\pi r}{\lambda}.
+```
+
+The refractive index `m` passed to a shape is **relative** to the host medium
+(``m = m_\text{particle}/m_\text{host}``); a complex `m` carries absorption in
+its imaginary part. The default wavelength in the API is ``\lambda = 2\pi``, for
+which ``k = 1`` and the size parameter equals the radius numerically — convenient
+for quick experiments.
+
+The truncation order needed for convergence grows with the size parameter
+(roughly ``n_{\max} \gtrsim x + c\,x^{1/3}``), which is why large particles are
+expensive and eventually require extended precision or the `stable` formulation;
+see [Choosing a solver](@ref).
+
+## Orientation
+
+Because ``T`` is an operator between VSWF bases, a rotation of the particle by
+Euler angles ``(\alpha, \beta, \gamma)`` acts on it through Wigner ``D``-functions
+(MTL Eq. (5.29)); [`rotate`](@ref) returns the rotated transition matrix without
+recomputing it from scratch.
+
+For randomly oriented particles the orientation-averaged transition matrix has a
+closed form (MTL Eq. (5.96)), implemented by
+[`RandomOrientationTransitionMatrix`](@ref); [`orientation_average`](@ref) instead
+integrates numerically against a user-supplied orientation distribution and
+converges to the analytic result. The orientation-averaged **scattering matrix**
+is expanded in generalised spherical functions with coefficients computed by
+[`expansion_coefficients`](@ref).
+
+## Far-field observables
+
+All of the following are methods of an [`AbstractTransitionMatrix`](@ref) (see
+[Post-processing](@ref) for runnable examples).
+
+### Amplitude (Jones) matrix
+
+For a fixed incidence direction ``\hat{\mathbf{n}}^{\text{inc}}`` and scattering
+direction ``\hat{\mathbf{n}}^{\text{sca}}``, the ``2\times 2`` complex amplitude
+matrix ``\mathbf{S}`` maps the incident to the scattered transverse field
+(MTL Eqs. (5.11)–(5.17)):
+
+```math
+\begin{bmatrix} E_\vartheta^{\text{sca}} \\ E_\varphi^{\text{sca}} \end{bmatrix}
+  = \frac{e^{ikr}}{r}\,\mathbf{S}\,
+    \begin{bmatrix} E_\vartheta^{\text{inc}} \\ E_\varphi^{\text{inc}} \end{bmatrix}.
+```
+
+[`amplitude_matrix`](@ref) returns ``\mathbf{S}`` for given incidence/scattering
+angles.
+
+### Phase (Mueller) matrix
+
+The real ``4\times 4`` phase matrix ``\mathbf{Z}`` propagates the Stokes vector
+and follows from ``\mathbf{S}`` (MTL Eqs. (2.106)–(2.121)); [`phase_matrix`](@ref)
+builds it from an amplitude matrix.
+
+### Cross sections, albedo, asymmetry
+
+The orientation-averaged scattering and extinction cross sections (MTL
+Eqs. (5.140)/(5.141) and (5.102)/(5.107)) are returned by
+[`scattering_cross_section`](@ref) and [`extinction_cross_section`](@ref); their
+difference is the [`absorption_cross_section`](@ref). The single-scattering
+[`albedo`](@ref) is ``\omega = C_\text{sca}/C_\text{ext}`` and the
+[`asymmetry_parameter`](@ref) is ``g = \langle\cos\Theta\rangle = \alpha_1^1/3``
+(MTL Eq. (4.92)), where ``\alpha_1^1`` is the first scattering-matrix expansion
+coefficient.
+
+### Scattering matrix
+
+For a macroscopically isotropic and mirror-symmetric ensemble the scattering
+matrix ``\mathbf{F}(\Theta)`` has the block-diagonal structure with up to ten
+independent elements ``(\alpha_1\!-\!\alpha_4, \beta_1\!-\!\beta_6)``;
+[`scattering_matrix`](@ref) evaluates it on a grid of scattering angles ``\Theta``.
+
+## Methods for building ``T``
+
+`TransitionMatrices.jl` provides several engines, all producing the same kind of
+transition matrix:
+
+| Method | Scatterers | Notes |
+| --- | --- | --- |
+| **Mie** (`bhmie`, `bhcoat`) | homogeneous & coated spheres | analytic; the reference for validation |
+| **EBCM** (null-field) | axisymmetric (spheroid, cylinder, Chebyshev) | Waterman's method; fast, but the standard `Float64` integrands lose precision for elongated particles — use `stable` or extended precision |
+| **IITM** (invariant imbedding) | axisymmetric *and* arbitrary (incl. prisms) | numerically robust where EBCM struggles; the only built-in route for non-axisymmetric shapes |
+| **Sh-matrix** (moment separation) | axisymmetric | separates the geometry from ``(\lambda, m)`` so parameter **sweeps** reuse one quadrature |
+
+Which one to pick, and when to enable `stable` or raise the precision, is covered
+in [Choosing a solver](@ref). The literature behind each method is collected in
+[Methods & references](@ref).
diff --git a/docs/src/usage.md b/docs/src/usage.md
index f065ffd..90e90b2 100644
--- a/docs/src/usage.md
+++ b/docs/src/usage.md
@@ -36,6 +36,19 @@ The following defines a Chebyshev particle with radius `r₀=1.0`, deformation c
 chebyshev = Chebyshev{Arb, Acb}(1.0, 0.1, 4, 1.5+0.01im);
 ```
 
+### Prism
+
+The following defines a hexagonal prism (`N = 6` base edges) with base edge
+length `a = 1.0`, height `h = 2.0`, and refractive index `m = 1.5 + 0.01im`. The
+number of base edges is a type parameter, supplied as the first argument:
+
+```julia
+prism = Prism(6, 1.0, 2.0, 1.5 + 0.01im);
+```
+
+A prism is **not** axisymmetric, so its T-matrix is built with the [`IITM`](@ref)
+solver rather than EBCM (see [Choosing a solver](@ref)).
+
 ## Calculate the T-Matrix
 
 After defining a shape, the T-Matrix is computed with `transition_matrix` (alias
@@ -66,10 +79,39 @@ For a fixed build, pass a fixed solver with its discretization:
 ```julia
 𝐓 = transition_matrix(spheroid, 2π, EBCM(10, 100))       # order 10, 100 quadrature points
 𝐓 = transition_matrix(cylinder, 2π, IITM(10, 20, 24))    # IITM: order, Nr, Nϑ
+𝐓 = transition_matrix(prism, 2π, IITM(10, 20, 24, 24))   # non-axisymmetric: add Nφ
 𝐓 = transition_matrix(spheroid, 2π, ShMatrix(10, 200))   # Sh-matrix moment separation
 ```
 
-### High-aspect-ratio spheroids (`stable`)
+## Choosing a solver
+
+All engines produce the same kind of transition matrix; pick by the scatterer
+and the regime:
+
+| Situation | Recommended |
+| --- | --- |
+| Sphere or coated sphere | Mie (`bhmie` / `bhcoat`) — analytic and exact |
+| Axisymmetric, moderate size and aspect ratio | bare `transition_matrix(s, λ)` (auto-converged classic EBCM) |
+| High-aspect-ratio **spheroid** | [`Iterative(EBCM; stable = true)`](@ref Iterative) (see the next section) |
+| Non-axisymmetric shape (e.g. a prism) | [`IITM`](@ref) — the only built-in route |
+| Axisymmetric but EBCM ill-conditioned | [`IITM`](@ref) |
+| Many wavelengths / refractive indices for **one** shape | [`ShMatrix`](@ref) / [`prepare_sh`](@ref) — prepare once, reuse |
+
+Two orthogonal knobs control accuracy in hard cases:
+
+- **`stable`** rebuilds the EBCM `𝐔`-matrix integrals with a cancellation-free
+  formulation. It is the right tool for high-aspect-ratio spheroids in `Float64`,
+  and is **spheroid-only** (see the next section).
+- **Extended precision** (`Double64`, `Float128`, `Arb`) raises the size
+  parameter EBCM can reach and lowers round-off floors in general. It works for
+  any shape but is several times slower. The two stack: `stable` with a
+  `Double64` element type reaches `~1e-25` at high aspect ratio.
+
+When in doubt, start with the bare call; if it fails to converge or loses
+accuracy, switch to `Iterative(EBCM; stable = true)` for spheroids, or to
+[`IITM`](@ref) otherwise.
+
+## High-aspect-ratio spheroids: the `stable` formulation
 
 For spheroids of high aspect ratio the standard EBCM surface integrals lose all
 precision in `Float64`: the irregular Riccati–Bessel products develop large
@@ -101,7 +143,7 @@ combine the two: `stable` with a `Double64` element type reaches `~1e-25` at hig
 aspect ratio, and also lowers the residual `Float64` round-off floor that appears
 as the refractive index approaches `s → 1`.
 
-### Parameter sweeps with the Sh-matrix method (`prepare_sh`)
+## Parameter sweeps with the Sh-matrix method (`prepare_sh`)
 
 When you need the T-matrix of one fixed shape at *many* wavelengths or
 refractive indices — a spectrum, a dispersion curve, a refractive-index scan —
@@ -156,7 +198,8 @@ still applies, but the `𝐔` reconstruction is not stabilized.
 
 ## Post-processing
 
-After getting the T-Matrix, you can calculate the far-field scattering properties using the following functions:
+After getting the T-Matrix, you can calculate the far-field scattering
+properties. Each function has a short `calc_*` alias:
 
 - [`amplitude_matrix`](@ref), a.k.a. `calc_S`
 - [`phase_matrix`](@ref), a.k.a. `calc_Z`
@@ -165,10 +208,64 @@ After getting the T-Matrix, you can calculate the far-field scattering propertie
 - [`absorption_cross_section`](@ref), a.k.a. `calc_Cabs`
 - [`albedo`](@ref), a.k.a. `calc_ω`
 - [`asymmetry_parameter`](@ref), a.k.a. `calc_g`
+- [`scattering_matrix`](@ref), a.k.a. `calc_F` (orientation-averaged)
 
-And the orientation-averaged scattering matrix:
+The orientation-averaged scalar quantities take the T-matrix and the wavelength:
+
+```julia
+𝐓 = transition_matrix(spheroid, 2π)
 
-- [`scattering_matrix`](@ref), a.k.a. `calc_F`
+Csca = scattering_cross_section(𝐓, 2π)
+Cext = extinction_cross_section(𝐓, 2π)
+Cabs = absorption_cross_section(𝐓, 2π)
+ω    = albedo(𝐓)                     # = Csca / Cext
+g    = asymmetry_parameter(𝐓, 2π)    # ⟨cos Θ⟩
+```
+
+For a **fixed scattering geometry**, the ``2\times2`` amplitude (Jones) matrix is
+evaluated at an incidence direction `(ϑᵢ, φᵢ)` and a scattering direction
+`(ϑₛ, φₛ)` (all angles in radians); the ``4\times4`` Mueller / phase matrix
+follows from it:
+
+```julia
+𝐒 = amplitude_matrix(𝐓, 0.0, 0.0, π / 2, 0.0; λ = 2π)   # forward-to-side geometry
+𝐙 = phase_matrix(𝐒)
+```
+
+The orientation-averaged scattering matrix is evaluated on a grid of scattering
+angles ``\Theta``:
+
+```julia
+θs = range(0, π; length = 181)
+𝐅 = scattering_matrix(𝐓, 2π, θs)     # one matrix per angle
+```
+
+## Orientation averaging
+
+For randomly oriented particles you usually want orientation-averaged
+quantities. There are two routes:
+
+```julia
+𝐓 = transition_matrix(spheroid, 2π)
+
+# analytic random-orientation average — fast and exact (Mishchenko et al. (2002), Eq. 5.96)
+𝐓ᵣ = RandomOrientationTransitionMatrix(𝐓)
+
+# numerical average over an orientation distribution p(α, β, γ); uniform here
+𝐓ₙ = orientation_average(𝐓, (α, β, γ) -> 1 / (8π^2); Nα = 40, Nβ = 40, Nγ = 1)
+```
+
+The numerical average converges to the analytic one as the angular grid is
+refined (see the *Orientation averaging* example in [Examples](@ref)). The cross
+sections and `scattering_matrix` above apply to either result.
+
+A single **fixed** orientation is obtained by rotating the T-matrix by Euler
+angles (Z-Y-Z), without recomputing it from the shape (`RotZYZ` is re-exported):
+
+```julia
+𝐓rot = rotate(𝐓, RotZYZ(0.3, 0.5, 0.0))
+𝐒    = amplitude_matrix(𝐓rot, 0.0, 0.0, π / 2, 0.0)
+```
 
 ## Differentiation
 
diff --git a/src/EBCM/axisymmetric.jl b/src/EBCM/axisymmetric.jl
index 103ff42..e861486 100644
--- a/src/EBCM/axisymmetric.jl
+++ b/src/EBCM/axisymmetric.jl
@@ -180,17 +180,17 @@ Parameters:
 - `nₘₐₓ`: the maximum order of the T-Matrix.
 - `Ng`: the number of Gauss-Legendre quadrature points to be used.
 - `stable`: when `true`, assemble the `𝐔`-matrix integrals with the
-  cancellation-free `F⁺` formulation of Somerville, Auguié & Le Ru, JQSRT 123
-  (2013). The standard integrands lose all precision for spheroids of high aspect
-  ratio (the irregular `χ_n·ψ_k` products develop huge Laurent terms that should
+  cancellation-free `F⁺` formulation of Somerville et al. (2013). The standard
+  integrands lose all precision for spheroids of high aspect
+  ratio (the irregular ``\chi_n\psi_k`` products develop huge Laurent terms that should
   cancel on integration but do not numerically); `stable=true` removes that
   cancellation, recovering a relative accuracy of about `1e-9` in `Float64`
   regardless of aspect ratio. It is **only valid for `Spheroid`** (the
   cancellation relies on the spheroid surface) and costs roughly 2–3× the default
   assembly, so it is opt-in. For it to help, `nₘₐₓ` must be large enough that
-  `nₘₐₓ+1 ≳ k·c + 15`, where `c` is the largest semi-axis — comparable to the
+  ``n_{\max}+1 \gtrsim k\cdot c + 15``, where `c` is the largest semi-axis — comparable to the
   order needed for convergence at that size anyway. The remaining `Float64`
-  round-off (in particular a `~1e-9` floor as the refractive index `s → 1`) is
+  round-off (in particular a `~1e-9` floor as the refractive index ``s \to 1``) is
   orthogonal to the cancellation and is lowered by using an extended-precision
   element type: `stable=true` with `Double64` reaches `~1e-25` at high aspect
   ratio.
@@ -204,7 +204,7 @@ function transition_matrix(s::AbstractAxisymmetricShape{T, CT}, λ, nₘₐₓ,
     stable && !(s isa Spheroid) &&
         throw(ArgumentError("stable=true is only valid for spheroids: the \
               cancellation-free F⁺ integrands rely on the spheroid surface making \
-              the divergent Laurent terms integrate to zero (Somerville et al. 2013)."))
+              the divergent Laurent terms integrate to zero (Somerville et al. (2013))."))
     𝐓 = Vector{Matrix{CT}}(undef, nₘₐₓ + 1)
     𝐓[1],
     cache = transition_matrix_m₀(s, λ, nₘₐₓ, Ng; zerofn = zerofn, reuse = true,
@@ -342,7 +342,7 @@ function ebcm_matrices_m₀(s::AbstractAxisymmetricShape{T, CT}, λ, nₘₐₓ,
     # For spheroids the off-diagonal `𝐔` integrands lose all precision at high
     # aspect ratio (the χ_n·ψ_k products have huge cancelling negative powers).
     # When `stable`, replace those products by the cancellation-free `F⁺` matrix
-    # (Somerville et al. 2013); the diagonal (n=n′) has no cancellation and `𝐏`
+    # (Somerville et al. (2013)); the diagonal (n=n′) has no cancellation and `𝐏`
     # is regular, so both keep the standard direct products. The x-independent
     # last-row series coefficients are computed once and shared across all points.
     Fmats = if stable
@@ -577,9 +577,9 @@ function _transition_matrix_m_core(m, s::AbstractAxisymmetricShape{T, CT}, k, n
                 if stable
                     F = Fmats[i]
                     xi = k * r[i]
-                    # Cancellation-free K¹/K² integrands (Somerville, Auguié & Le Ru,
-                    # JQSRT 123 (2013), Eqs. 53–54): K¹ ← [x·χ_n·ψ′_{n′}]⁺/x (their
-                    # Eq. 59), K² ← [x·χ′_n·ψ_{n′}]⁺/x (their Eq. 60).
+                    # Cancellation-free K¹/K² integrands (Somerville et al. (2013),
+                    # Eqs. (53)–(54)): K¹ ← [x·χ_n·ψ′_{n′}]⁺/x (their
+                    # Eq. (59)), K² ← [x·χ′_n·ψ_{n′}]⁺/x (their Eq. (60)).
                     UK₁ += w[i] * k * r′[i] * 𝜋[i, n] * d[i, n′] *
                            (_xχψ′⁺_mat(F, n, n′) / xi)
                     UK₂ += w[i] * k * r′[i] * 𝜋[i, n] * d[i, n′] *
diff --git a/src/EBCM/shmatrix.jl b/src/EBCM/shmatrix.jl
index 146e692..374b4c9 100644
--- a/src/EBCM/shmatrix.jl
+++ b/src/EBCM/shmatrix.jl
@@ -22,9 +22,10 @@
 # analytically, so — computed at high precision — they contribute ≈0 and the
 # catastrophic cancellation in the irregular `χ_n·ψ_{n′}` products never forms.
 #
-# References:
-# - V. G. Farafonov, Sh-matrix family (parameter separation).
-# - Somerville, Auguié & Le Ru, JQSRT 123 (2013) 153 (the `F⁺` projection).
+# References (full entries in the README "Methods & references"):
+# - Farafonov et al. (2003) — the size/material vs geometry parameter separation;
+#   Petrov et al. (2007, 2009) — the "Sh-matrix" name and formalism.
+# - Somerville et al. (2013) — the `F⁺` projection reused for the radial integrals.
 
 # ── Radial power-series term tables (reuse `_β`, `_γ`) ────────────────────────
 #
@@ -33,12 +34,12 @@
 #   ψ′_n(z) : base = n,     c_b = β_{n,b}·(n+1+2b)
 #   χ_n(z)  : base = -n,    c_b = γ_{n,b}
 #   χ′_n(z) : base = -n-1,  c_b = γ_{n,b}·(2b-n)
-"""
+@doc raw"""
     _sh_series(kind, n, B, R) -> (base, coeffs::Vector{R})
 
 Coefficients `c_b` (`b = 0:B-1`) and the lowest exponent `base` of the power
 series of a Riccati–Bessel function `kind ∈ (:ψ, :ψ′, :χ, :χ′)` of order `n`,
-i.e. `f_n(z) = Σ_b coeffs[b+1] · z^{base + 2b}`.
+i.e. ``f_n(z) = \sum_b \text{coeffs}[b+1] \cdot z^{\text{base} + 2b}``.
 """
 function _sh_series(kind::Symbol, n::Integer, B::Integer, ::Type{R}) where {R}
     coeffs = Vector{R}(undef, B)
@@ -83,11 +84,11 @@ function _sh_coeff_tables(nmax::Integer, B::Integer, ::Type{R}) where {R}
 end
 
 # ── Per-evaluation power weighting ────────────────────────────────────────────
-"""
+@doc raw"""
     _sh_weighted(tbl, nmin, nmax, z) -> OffsetVector of (base, w)
 
-Fold the argument power `z^{base+2b}` into the series coefficients, giving
-`w[b+1] = coeffs[b+1]·z^{base+2b}`, with `z = k` for the regular factor `f(kr)`
+Fold the argument power ``z^{\text{base}+2b}`` into the series coefficients, giving
+``w[b+1] = \text{coeffs}[b+1]\cdot z^{\text{base}+2b}``, with `z = k` for the regular factor `f(kr)`
 or `z = s·k` for the internal factor `g(s·kr)`. Built once per `(k, mᵣ)` point
 per order and reused across all integrands. Powers advance incrementally — no
 `^` in the hot loop. The `base` is kept so `_sh_conv` can index the
@@ -112,15 +113,15 @@ function _sh_weighted(tbl, nmin::Integer, nmax::Integer, z::Number)
 end
 
 # ── Folded coefficient×moment reconstruction ──────────────────────────────────
-"""
+@doc raw"""
     _sh_conv(f, g, M, shift) -> Complex
 
-Reconstruct `∫ (geometry) · f_n(kr) · g_{n′}(s·kr) · r^{shift} dϑ` from the
+Reconstruct ``\int (\text{geometry}) \cdot f_n(kr) \cdot g_{n'}(s\cdot kr) \cdot r^{\text{shift}} \, d\vartheta`` from the
 weighted radial terms `f = (basef, u)`, `g = (baseg, v)` (see `_sh_weighted`)
 and the shape-moment lookup `M(q)`. The double sum is folded along the
-anti-diagonals `j = b+c`: the `r`-power is `q₀ + 2j` (`q₀ = basef+baseg+shift`)
+anti-diagonals `j = b+c`: the `r`-power is ``q_0 + 2j`` (``q_0 = \text{basef}+\text{baseg}+\text{shift}``)
 and depends only on `j`, so the moment `M` is fetched `O(B)` times instead of
-`O(B²)`, and the inner convolution `Σ_{b+c=j} u[b]·v[c]` reuses the precomputed
+`O(B²)`, and the inner convolution ``\sum_{b+c=j} u[b]\cdot v[c]`` reuses the precomputed
 weighted terms with no `^` and no per-term coefficient work.
 """
 @inline function _sh_conv(f, g, M, shift::Integer)
@@ -151,23 +152,23 @@ end
 _sh_widen(s, ::Type{R}) where {R} = s
 
 # ── Shape moments for one azimuthal index `m` ────────────────────────────────
-"""
+@doc raw"""
     _sh_geometry(shape, Ng, qlo, qhi; momtype, store) -> NamedTuple
 
 Compute the `m`-independent geometry once: the (high-precision) Gauss–Legendre
 nodes/weights, the polar angles `ϑ`, and the two `r`-power weight tables shared
 by every azimuthal block,
 
-  WR[i,q] = wᵢ r′ᵢ rᵢ^q     (for the `τd`/`πd` families)
-  Wq[i,q] = wᵢ rᵢ^q         (for the `ππττ` family, which carries no `r′`)
+  ``WR[i,q] = w_i r'_i r_i^q``     (for the `τd`/`πd` families)
+  ``Wq[i,q] = w_i r_i^q``          (for the `ππττ` family, which carries no `r′`)
 
-`gaussquad` and the `rᵢ^q` table do not depend on `m`, so hoisting them out of
+`gaussquad` and the ``r_i^q`` table do not depend on `m`, so hoisting them out of
 the per-`m` loop avoids recomputing them `nmax+1` times.
 
 The weights are also returned cast to the fast accumulation type
 `Tf = store<:IEEEFloat ? store : momtype` as `WRf`/`Wqf`: the catastrophic
 cancellation lives entirely in the negative-`r`-power (`q<0`) moments, so the
-well-conditioned `q≥0` bulk can be contracted in hardware precision while only
+well-conditioned ``q \ge 0`` bulk can be contracted in hardware precision while only
 the `q<0` part keeps the slow `momtype` accumulation (see `_sh_moments_m`).
 """
 function _sh_geometry(shape, Ng::Integer, qlo::Integer, qhi::Integer;
@@ -209,29 +210,29 @@ function _sh_geometry(shape, Ng::Integer, qlo::Integer, qhi::Integer;
     return (; sym, ng, ϑ, WR, Wq, WRf, Wqf, qlo, qhi)
 end
 
-"""
+@doc raw"""
     _sh_moments_m(geom, m, nmax; store) -> NamedTuple
 
 Compute the shape-only moment tables for azimuthal index `m` from the shared
 `_sh_geometry` `geom`, stored as `store` (default `Float64`). Families
 (sum over `i in 1:ng`, mirroring `ebcm_matrices_m₀`):
 
-  Mτd[n,n′,q]  = Σ WR[i,q] τ[i,n] d[i,n′]
-  Mπd[n,n′,q]  = Σ WR[i,q] π[i,n] d[i,n′]   (only `m>0`)
-  Mππττ[n,q]   = Σ Wq[i,q] (π[i,n]²+τ[i,n]²)
+  ``M_{\tau d}[n,n',q]  = \sum_i WR[i,q]\,\tau[i,n]\,d[i,n']``
+  ``M_{\pi d}[n,n',q]  = \sum_i WR[i,q]\,\pi[i,n]\,d[i,n']``   (only `m>0`)
+  ``M_{\pi\pi\tau\tau}[n,q]   = \sum_i Wq[i,q]\,(\pi[i,n]^2+\tau[i,n]^2)``
 
-`Mdτ[n,n′,q] = Mτd[n′,n,q]` exactly (swap the two orders), so it is not stored —
+``Md\tau[n,n',q] = M\tau d[n',n,q]`` exactly (swap the two orders), so it is not stored —
 the assembly reads `Mτd` with the indices transposed. The angular products are
 formed once per `(n,n′)` instead of once per `q`.
 
 **Precision split.** The catastrophic cancellation that motivates the
 high-precision accumulation lives entirely in the negative-`r`-power (`q<0`)
 moments — for a spheroid those vanish analytically and must reach ≈0 rather than
-round-off. The `q≥0` bulk is well conditioned, so it is contracted in the fast
+round-off. The ``q \ge 0`` bulk is well conditioned, so it is contracted in the fast
 type `Tf` (hardware float when `store` is one), while only the `q<0` slice keeps
 the slow `momtype` accumulation. This is where most of the precompute time goes,
 so the split is the dominant speedup; it leaves results unchanged because the
-`q≥0` moments never needed the extra precision.
+``q \ge 0`` moments never needed the extra precision.
 """
 function _sh_moments_m(geom, m::Integer, nmax::Integer;
         store::Type{Ts} = Float64, vanish_negative::Bool = false) where {Ts}
@@ -447,13 +448,13 @@ function _sh_matrices_m₀(mom, coeffs, k::Real, s::Number, nmax::Integer,
 end
 
 # ── m>0 assembly from moments (mirrors `_transition_matrix_m_core`) ───────────
-"""
+@doc raw"""
     _sh_matrices_m(mom, coeffs, m, k, s, nmax, CT) -> (𝐏, 𝐔)
 
-Reconstruct the `m`-th `𝐏` and `𝐔` blocks (`2nn × 2nn`, `nn = nmax-m+1` for
-`m≥1`) from the precomputed moment tables `mom` for that `m` and the
+Reconstruct the `m`-th `𝐏` and `𝐔` blocks (``2nn \times 2nn``, ``nn = \text{nmax}-m+1`` for
+``m \ge 1``) from the precomputed moment tables `mom` for that `m` and the
 size/material-independent coefficient tables `coeffs`. The `K`-blocks use the
-`M^{πd}` family; the `L`-blocks reuse the same reconstruction as the `m=0` case
+``M^{\pi d}`` family; the `L`-blocks reuse the same reconstruction as the `m=0` case
 (with the `m`-dependent Wigner functions baked into `mom`). The algebra
 reproduces `_transition_matrix_m_core` term by term.
 """
@@ -587,7 +588,7 @@ struct ShPreparation{ST, MT, CF}
     stable::Bool
 end
 
-"""
+@doc raw"""
     prepare_sh(shape, nmax, Ng; B, momtype, store) -> ShPreparation
 
 Precompute the Sh-matrix shape-only moment tables for `shape` up to order `nmax`
@@ -599,7 +600,7 @@ Keyword arguments:
 
 - `B`: number of radial power-series terms per Riccati–Bessel function. It sets
   the largest size parameter the reconstruction resolves (the series, like the
-  `F⁺` evaluation, needs more terms as `k·rₘₐₓ·|mᵣ|` grows); it is fixed at
+  `F⁺` evaluation, needs more terms as ``k\cdot r_{\max}\cdot|m_r|`` grows); it is fixed at
   preparation time because it determines the moment band. Default `max(30,
   nmax+15)`.
 - `momtype`: precision used to accumulate the moments. The default (`nothing`)
diff --git a/src/EBCM/stabilization.jl b/src/EBCM/stabilization.jl
index bdd5308..eaa7fd1 100644
--- a/src/EBCM/stabilization.jl
+++ b/src/EBCM/stabilization.jl
@@ -5,16 +5,16 @@
 # and a regular one ψ_k. For a spheroid, the negative-power part of their
 # Laurent expansion integrates *exactly* to zero, but evaluated directly it
 # dominates the integrand by many orders of magnitude and destroys precision in
-# `Float64` (Somerville, Auguié & Le Ru, JQSRT 113:524 (2012)).
+# `Float64` (Somerville et al. (2012)).
 #
-# The fix (Somerville et al., JQSRT 123:153 (2013)) replaces χ_n(x)·ψ_k(sx) by
-# F⁺_{nk}(s,x)/x, where (their Eq. 45)
+# The fix (Somerville et al. (2013)) replaces χ_n(x)·ψ_k(sx) by
+# F⁺_{nk}(s,x)/x, where (their Eq. (45))
 #
 #     F_{nk}(s,x)  = x · χ_n(x) · ψ_k(sx),
 #     F⁺_{nk}(s,x) = P⁺[F_{nk}]   (the part with non-negative powers of x),
 #
 # computed directly from the power series so the cancelling negative powers are
-# never formed (their Eq. 46 and §4.1).
+# never formed (their Eq. (46) and §4.1).
 #
 # Power series (DLMF 10.53.1/10.53.2). This package uses ψ_k = z·j_k and the
 # convention χ_n = z·y_n (NO minus sign; see `ricattibessely`):
@@ -42,7 +42,7 @@
 """
     _odd_dfact(m) -> BigInt
 
-Odd double factorial `m!!` for odd `m ≥ -1`, with the convention `(-1)!! = 1`.
+Odd double factorial ``m!!`` for odd `m ≥ -1`, with the convention ``(-1)!! = 1``.
 """
 function _odd_dfact(m::Integer)
     isodd(m) || throw(ArgumentError("_odd_dfact expects an odd integer, got $m"))
@@ -56,10 +56,10 @@ function _odd_dfact(m::Integer)
     return r
 end
 
-"""
+@doc raw"""
     _γ(n, a) -> Rational{BigInt}
 
-`γ_{n,a}`: coefficient of `z^{2a-n}` in the series of `χ_n(z) = -z y_n(z)`
+``\gamma_{n,a}``: coefficient of ``z^{2a-n}`` in the series of ``\chi_n(z) = -z y_n(z)``
 (DLMF 10.53.2).
 """
 function _γ(n::Integer, a::Integer)
@@ -77,10 +77,10 @@ function _γ(n::Integer, a::Integer)
     end
 end
 
-"""
+@doc raw"""
     _β(k, b) -> Rational{BigInt}
 
-`β_{k,b}`: coefficient of `z^{k+1+2b}` in the series of `ψ_k(z) = z j_k(z)`
+``\beta_{k,b}``: coefficient of ``z^{k+1+2b}`` in the series of ``\psi_k(z) = z j_k(z)``
 (DLMF 10.53.1).
 """
 function _β(k::Integer, b::Integer)
@@ -89,15 +89,15 @@ function _β(k::Integer, b::Integer)
     return sgn // den
 end
 
-"""
+@doc raw"""
     _F⁺_coeffs(n, k, s; nterms, prec) -> (qmin, coeffs)
 
-High-precision coefficients of `F⁺_{nk}(s,x) = Σ c_q x^{2q+k-n+2}` for
+High-precision coefficients of ``F^+_{nk}(s,x) = \sum_q c_q x^{2q+k-n+2}`` for
 `q = qmin, …, qmin+nterms-1`, returned as a `Vector{Complex{BigFloat}}`.
 `s` is the (generally complex) relative refractive index. The sum defining each
 `c_q` is accumulated in `BigFloat` at precision `prec` bits to defeat the
-`(2n-1)!!`-scale cancellation. This `BigFloat` accumulation also subsumes the
-additional `s≈1` cancellation that Somerville et al. 2013 handle analytically in
+``(2n-1)!!``-scale cancellation. This `BigFloat` accumulation also subsumes the
+additional ``s \approx 1`` cancellation that Somerville et al. (2013) handle analytically in
 their Appendix B (the leading coefficients vanish at `s=1`): the base precision
 has many digits of margin for it, so no special treatment is needed — the s≈1
 accuracy floor of the stabilized assembly is instead the `Float64` round-off in
@@ -124,17 +124,17 @@ function _F⁺_coeffs(n::Integer, k::Integer, s::Number; nterms::Integer = 48,
     return qmin, coeffs
 end
 
-"""
+@doc raw"""
     _eval_F⁺(qmin, coeffs, n, k, x::T) -> Complex{T}
 
-Evaluate `F⁺_{nk}(s,x) = Σ_q c_q x^{2q+k-n+2}` at `x` in floating type `T`,
+Evaluate ``F^+_{nk}(s,x) = \sum_q c_q x^{2q+k-n+2}`` at `x` in floating type `T`,
 converting the high-precision `coeffs` to `Complex{T}` first.
 
-The terms `c_q x^{2q+…}` first grow then decay; the sum is truncated once it has
+The terms ``c_q x^{2q+\ldots}`` first grow then decay; the sum is truncated once it has
 started to decay and a term is negligible relative to the running total
-(Somerville et al. 2013, §4.1). This both avoids summing past convergence and,
+(Somerville et al. (2013), §4.1). This both avoids summing past convergence and,
 crucially, prevents the monomial `xp = x^{2q+…}` from overflowing to `Inf`/`NaN`
-for large `x` (it would reach `x^{2·length(coeffs)}`). The series is only well
+for large `x` (it would reach ``x^{2\cdot\text{length(coeffs)}}``). The series is only well
 conditioned in `Float64` when the order `n` exceeds the argument `x`; in that
 regime convergence is reached long before `xp` overflows.
 """
@@ -162,12 +162,12 @@ function _eval_F⁺(qmin::Integer, coeffs::AbstractVector{<:Complex}, n::Integer
     return val
 end
 
-"""
+@doc raw"""
     F⁺(n, k, s, x; nterms, prec) -> Complex
 
-Numerically stable `F⁺_{nk}(s,x)` (Somerville et al. 2013, Eq. 45–46),
-evaluated in the precision of `x`. `F⁺/x` is the cancellation-free replacement
-for `χ_n(x)·ψ_k(sx)` in the spheroid EBCM `U`-matrix integrand.
+Numerically stable ``F^+_{nk}(s,x)`` (Somerville et al. (2013), Eqs. (45)–(46)),
+evaluated in the precision of `x`. ``F^+/x`` is the cancellation-free replacement
+for ``\chi_n(x)\psi_k(sx)`` in the spheroid EBCM ``\mathbf{U}``-matrix integrand.
 """
 function F⁺(n::Integer, k::Integer, s::Number, x::T;
         nterms::Integer = max(24, ceil(Int, 2 * abs(s) * x + 16)),
@@ -176,38 +176,38 @@ function F⁺(n::Integer, k::Integer, s::Number, x::T;
     return _eval_F⁺(qmin, coeffs, n, k, x)
 end
 
-# ── Cancellation-free "modified" Bessel products via F⁺ (Somerville 2013) ─────
+# ── Cancellation-free "modified" Bessel products via F⁺ (Somerville et al. (2013)) ─
 # Each returns the P⁺ (non-negative-power) part of a product appearing in the
 # spheroid EBCM U-matrix integrand, assembled from F⁺ at shifted indices through
-# the Riccati–Bessel recurrences (their Eqs. 57–62). All carry the leading `x`
+# the Riccati–Bessel recurrences (their Eqs. (57)–(62)). All carry the leading `x`
 # factor, like `F⁺` itself. For small enough n−k these equal the full products
 # (no negative powers), which the tests exploit for validation.
 
-"""
+@doc raw"""
     _xχψ′⁺(n, k, s, x) -> Complex
 
-`[x·χ_n(x)·ψ′_k(sx)]⁺` (Somerville 2013, Eq. 59), from
-`(2k+1)ψ′_k(sx) = (k+1)ψ_{k-1}(sx) − k ψ_{k+1}(sx)`.
+``[x\chi_n(x)\psi'_k(sx)]^+`` (Somerville et al. (2013), Eq. (59)), from
+``(2k+1)\psi'_k(sx) = (k+1)\psi_{k-1}(sx) - k\psi_{k+1}(sx)``.
 """
 function _xχψ′⁺(n, k, s, x; kw...)
     ((k + 1) * F⁺(n, k - 1, s, x; kw...) - k * F⁺(n, k + 1, s, x; kw...)) / (2k + 1)
 end
 
-"""
+@doc raw"""
     _xχ′ψ⁺(n, k, s, x) -> Complex
 
-`[x·χ′_n(x)·ψ_k(sx)]⁺` (Somerville 2013, Eq. 60), from
-`(2n+1)χ′_n(x) = (n+1)χ_{n-1}(x) − n χ_{n+1}(x)`.
+``[x\chi'_n(x)\psi_k(sx)]^+`` (Somerville et al. (2013), Eq. (60)), from
+``(2n+1)\chi'_n(x) = (n+1)\chi_{n-1}(x) - n\chi_{n+1}(x)``.
 """
 function _xχ′ψ⁺(n, k, s, x; kw...)
     ((n + 1) * F⁺(n - 1, k, s, x; kw...) - n * F⁺(n + 1, k, s, x; kw...)) / (2n + 1)
 end
 
-"""
+@doc raw"""
     _L⁷⁺(n, k, s, x) -> Complex
 
-`[x·(χ′_n ψ′_k + n(n+1) χ_n ψ_k /(s x²))]⁺` (Somerville 2013, Eq. 62) — the
-cancellation-free radial factor of the `L⁷` integrand.
+``[x(\chi'_n \psi'_k + n(n+1)\chi_n\psi_k/(sx^2))]^+`` (Somerville et al. (2013), Eq. (62)) — the
+cancellation-free radial factor of the ``L^7`` integrand.
 """
 function _L⁷⁺(n, k, s, x; kw...)
     Fmm = F⁺(n - 1, k - 1, s, x; kw...)
@@ -218,11 +218,11 @@ function _L⁷⁺(n, k, s, x; kw...)
             (n - k) * ((n + 1) * Fmp + n * Fpm)) / ((2n + 1) * (2k + 1))
 end
 
-"""
+@doc raw"""
     _L⁸⁺(n, k, s, x) -> Complex
 
-`[x·(χ′_n ψ′_k + k(k+1) χ_n ψ_k /(s x²))]⁺` (Somerville 2013, Eq. 61) — the
-cancellation-free radial factor of the `L⁸` integrand.
+``[x(\chi'_n \psi'_k + k(k+1)\chi_n\psi_k/(sx^2))]^+`` (Somerville et al. (2013), Eq. (61)) — the
+cancellation-free radial factor of the ``L^8`` integrand.
 """
 function _L⁸⁺(n, k, s, x; kw...)
     Fmm = F⁺(n - 1, k - 1, s, x; kw...)
@@ -233,11 +233,11 @@ function _L⁸⁺(n, k, s, x; kw...)
             (k - n) * ((k + 1) * Fpm + k * Fmp)) / ((2n + 1) * (2k + 1))
 end
 
-"""
+@doc raw"""
     _F⁺_lastrow(s, N, xmax; prec, nterms) -> Vector
 
 Precompute the *x-independent* power-series coefficients of the F⁺-matrix last row
-`n = N+1` (the entries with `n-k ≥ 4`, `n+k` even, that seed the recursion). These
+``n = N+1`` (the entries with ``n-k \ge 4``, ``n+k`` even, that seed the recursion). These
 are the only expensive (`BigFloat`) part of `_F⁺_matrix`, so computing them once
 and passing them to every per-quadrature-point `_F⁺_matrix` call removes the
 dominant cost of the stabilized assembly. `nterms` is sized for the *largest*
@@ -260,28 +260,28 @@ end
 @doc raw"""
     _F⁺_matrix(s, x, N; prec, nterms, lastrow) -> Matrix{Complex{T}}
 
-Stably evaluate the whole matrix of `F⁺_{nk}(s,x)` for `0 ≤ n,k ≤ N+1`, returned
-with a 1-based offset (`F[n+1, k+1]` holds `F⁺_{nk}`). Only `n+k`-even entries are
+Stably evaluate the whole matrix of ``F^+_{nk}(s,x)`` for ``0 \le n,k \le N+1``, returned
+with a 1-based offset (`F[n+1, k+1]` holds ``F^+_{nk}``). Only `n+k`-even entries are
 filled; the rest stay zero (they are only needed for `m>0`).
 
 The bare power series is well conditioned in `Float64` *only* where the order
 exceeds the argument, so it cannot be used alone for the small-`n`/large-`x`
-corner. This routine combines the three regimes of Somerville, Auguié & Le Ru,
-JQSRT 123 (2013) §4.2, so the result is accurate for all `x`:
+corner. This routine combines the three regimes of Somerville et al. (2013),
+§4.2, so the result is accurate for all `x`:
 
-  * `n ≤ k+2` — `F⁺ = F = x·χ_n(x)·ψ_k(sx)`, the direct Riccati–Bessel product
+  * ``n \le k+2`` — ``F^+ = F = x\chi_n(x)\psi_k(sx)``, the direct Riccati–Bessel product
     (no negative powers ⇒ no cancellation, stable for all `x`);
-  * last row `n = N+1`, deep entries `n-k ≥ 4` — the power series (Eq. 46), which
-    is accurate because `N` is taken large enough that `N+1 ≫ x`;
-  * `n ≥ k+4` — the stable recursion (Eq. 51, scheme (c) of their Fig. 3), filled
-    along diagonals `j = n-k = 4, 6, …` from high `n` down, seeded by the last row
-    and the `n = k+2` sub-diagonal:
+  * last row ``n = N+1``, deep entries ``n-k \ge 4`` — the power series (their Eq. (46)), which
+    is accurate because `N` is taken large enough that ``N+1 \gg x``;
+  * ``n \ge k+4`` — the stable recursion (their Eq. (51), scheme (c) of their Fig. 3), filled
+    along diagonals ``j = n-k = 4, 6, \ldots`` from high `n` down, seeded by the last row
+    and the ``n = k+2`` sub-diagonal:
 
     ```math
     F^+_{n,k} = \frac{2k+3}{s(2n+1)}\,(F^+_{n+1,k+1} + F^+_{n-1,k+1}) - F^+_{n,k+2}.
     ```
 
-`N` must satisfy roughly `N+1 ≳ x + 15` for the last-row seed (and hence the whole
+`N` must satisfy roughly ``N+1 \gtrsim x + 15`` for the last-row seed (and hence the whole
 bottom-left block) to be accurate; this is comparable to the multipole order
 needed for convergence anyway.
 """
@@ -325,7 +325,7 @@ function _F⁺_matrix(s::Number, x::T, N::Integer;
         F[M + 1, k + 1] = _eval_F⁺(qmin, coeffs, M, k, x)
     end
 
-    # (3) stable recursion (Eq. 51, scheme c): diagonals j = 4, 6, … filled from
+    # (3) stable recursion (their Eq. (51), scheme c): diagonals j = 4, 6, … filled from
     #     high n down; each entry uses the same-diagonal entry at n+1 and the j-2
     #     diagonal at (n-1,k+1) and (n,k+2), all already known.
     @inbounds for j in 4:2:M
@@ -346,13 +346,13 @@ end
 # matrix returned by `_F⁺_matrix` (`F[n+1,k+1] = F⁺_{nk}`).
 @inline _Fp(F, n, k) = @inbounds F[n + 1, k + 1]
 
-"`[x·χ_n·ψ′_k]⁺` (Somerville et al., JQSRT 123 (2013), Eq. 59) from a precomputed F⁺ matrix."
+raw"``[x\chi_n\psi'_k]^+`` (Somerville et al. (2013), Eq. (59)) from a precomputed F⁺ matrix."
 _xχψ′⁺_mat(F, n, k) = ((k + 1) * _Fp(F, n, k - 1) - k * _Fp(F, n, k + 1)) / (2k + 1)
 
-"`[x·χ′_n·ψ_k]⁺` (Somerville et al., JQSRT 123 (2013), Eq. 60) from a precomputed F⁺ matrix."
+raw"``[x\chi'_n\psi_k]^+`` (Somerville et al. (2013), Eq. (60)) from a precomputed F⁺ matrix."
 _xχ′ψ⁺_mat(F, n, k) = ((n + 1) * _Fp(F, n - 1, k) - n * _Fp(F, n + 1, k)) / (2n + 1)
 
-"L⁷ radial factor (Somerville et al., JQSRT 123 (2013), Eq. 62) from a precomputed F⁺ matrix."
+raw"``L^7`` radial factor (Somerville et al. (2013), Eq. (62)) from a precomputed F⁺ matrix."
 function _L⁷⁺_mat(F, n, k)
     Fmm = _Fp(F, n - 1, k - 1);
     Fpp = _Fp(F, n + 1, k + 1)
@@ -362,7 +362,7 @@ function _L⁷⁺_mat(F, n, k)
             (n - k) * ((n + 1) * Fmp + n * Fpm)) / ((2n + 1) * (2k + 1))
 end
 
-"L⁸ radial factor (Somerville et al., JQSRT 123 (2013), Eq. 61) from a precomputed F⁺ matrix."
+raw"``L^8`` radial factor (Somerville et al. (2013), Eq. (61)) from a precomputed F⁺ matrix."
 function _L⁸⁺_mat(F, n, k)
     Fmm = _Fp(F, n - 1, k - 1);
     Fpp = _Fp(F, n + 1, k + 1)
@@ -510,7 +510,7 @@ end
     using TransitionMatrices: Spheroid, transition_matrix, scattering_cross_section,
                               extinction_cross_section
 
-    # Somerville, Auguié & Le Ru, JQSRT 123 (2013), Table 2 (Model 1 of their
+    # Somerville et al. (2013), Table 2 (Model 1 of their
     # Ref. [35]): prolate, s = 1.55+0.01i, aspect h = 4, xmax = k·c = 10.079368.
     # Reference (their New DP / AP): Qsca = 3.21290554203156, Qext = 3.36721292620922.
     c = 10.079368
diff --git a/src/IITM/arbitrary.jl b/src/IITM/arbitrary.jl
index 05e29a4..18fe5d3 100644
--- a/src/IITM/arbitrary.jl
+++ b/src/IITM/arbitrary.jl
@@ -172,11 +172,11 @@ function transition_matrix_iitm(s::AbstractShape{T, CT}, λ, nₘₐₓ, Nr, Nϑ
         𝐐ₕⱼ = im * k * transpose(𝐇) * 𝐐 * 𝐉
         𝐐ₕₕ = im * k * transpose(𝐇) * 𝐐 * 𝐇
 
-        # Eq (97) in Johnson (1988) Note: 𝐐ⱼₕ'==𝐐ₕⱼ only holds for spheres
-        # Eq (38) in Bi et al. (2013)
-        # Eq (2.40) in Doicu & Wriedt (2018)
-        # Eq (5.71) in Hu (2018)
-        # Eq (4.2.36) in Sun et al. (2019) Note: incorrect multiplication order
+        # Johnson (1988), Eq. (97). Note: 𝐐ⱼₕ'==𝐐ₕⱼ only holds for spheres
+        # Bi et al. (2013), Eq. (38)
+        # Doicu & Wriedt (2018), Eq. (2.40)
+        # Hu (2018), Eq. (5.71)
+        # Sun et al. (2019), Eq. (4.2.36). Note: incorrect multiplication order
         𝐓 = 𝐐ⱼⱼ + (𝐈 + 𝐐ⱼₕ) * _iitm_ldiv(𝐈 - 𝐓 * 𝐐ₕₕ, 𝐓) * (𝐈 + 𝐐ₕⱼ)
     end
 
diff --git a/src/IITM/axisymmetric.jl b/src/IITM/axisymmetric.jl
index 097af16..d8ab495 100644
--- a/src/IITM/axisymmetric.jl
+++ b/src/IITM/axisymmetric.jl
@@ -157,11 +157,11 @@ function transition_matrix_iitm(s::AbstractAxisymmetricShape{T, CT}, λ, nₘₐ
             𝐐ₕⱼ = im * k * transpose(𝐇ᵥ) * 𝐐 * 𝐉ᵥ
             𝐐ₕₕ = im * k * transpose(𝐇ᵥ) * 𝐐 * 𝐇ᵥ
 
-            # Eq (97) in Johnson (1988) Note: 𝐐ⱼₕ'==𝐐ₕⱼ only holds for spheres
-            # Eq (38) in Bi et al. (2013)
-            # Eq (2.40) in Doicu & Wriedt (2018)
-            # Eq (5.71) in Hu (2018)
-            # Eq (4.2.36) in Sun et al. (2019) Note: incorrect multiplication order
+            # Johnson (1988), Eq. (97). Note: 𝐐ⱼₕ'==𝐐ₕⱼ only holds for spheres
+            # Bi et al. (2013), Eq. (38)
+            # Doicu & Wriedt (2018), Eq. (2.40)
+            # Hu (2018), Eq. (5.71)
+            # Sun et al. (2019), Eq. (4.2.36). Note: incorrect multiplication order
             Ts[m + 1] = 𝐐ⱼⱼ +
                         (𝐈 + 𝐐ⱼₕ) * _iitm_ldiv(𝐈 - Ts[m + 1] * 𝐐ₕₕ, Ts[m + 1]) *
                         (𝐈 + 𝐐ₕⱼ)
diff --git a/src/IITM/nfold.jl b/src/IITM/nfold.jl
index 86283d5..5bbca90 100644
--- a/src/IITM/nfold.jl
+++ b/src/IITM/nfold.jl
@@ -187,11 +187,11 @@ function transition_matrix_iitm(s::AbstractNFoldShape{N, T, CT}, λ, nₘₐₓ,
             𝐐ₕⱼ = im * k * transpose(𝐇) * 𝐐 * 𝐉
             𝐐ₕₕ = im * k * transpose(𝐇) * 𝐐 * 𝐇
 
-            # Eq (97) in Johnson (1988) Note: 𝐐ⱼₕ'==𝐐ₕⱼ only holds for spheres
-            # Eq (38) in Bi et al. (2013)
-            # Eq (2.40) in Doicu & Wriedt (2018)
-            # Eq (5.71) in Hu (2018)
-            # Eq (4.2.36) in Sun et al. (2019) Note: incorrect multiplication order
+            # Johnson (1988), Eq. (97). Note: 𝐐ⱼₕ'==𝐐ₕⱼ only holds for spheres
+            # Bi et al. (2013), Eq. (38)
+            # Doicu & Wriedt (2018), Eq. (2.40)
+            # Hu (2018), Eq. (5.71)
+            # Sun et al. (2019), Eq. (4.2.36). Note: incorrect multiplication order
             𝐓 .= 𝐐ⱼⱼ + (𝐈 + 𝐐ⱼₕ) * _iitm_ldiv(𝐈 - 𝐓 * 𝐐ₕₕ, 𝐓) * (𝐈 + 𝐐ₕⱼ)
         end
     end
diff --git a/src/Mie/MieTransitionMatrix.jl b/src/Mie/MieTransitionMatrix.jl
index d082a9b..9f034f7 100644
--- a/src/Mie/MieTransitionMatrix.jl
+++ b/src/Mie/MieTransitionMatrix.jl
@@ -1,5 +1,5 @@
 @doc raw"""
-According to Eq. (5.42) -- Eq. (5.44) in Mishchenko et al. (2002), the T-Matrix for a Mie particle can be written as:
+According to Mishchenko et al. (2002), Eqs. (5.42)–(5.44), the T-Matrix for a Mie particle can be written as:
 
 ```math
 \begin{array}{l}
diff --git a/src/Mie/bhcoat.jl b/src/Mie/bhcoat.jl
index dc2f50c..c675b48 100644
--- a/src/Mie/bhcoat.jl
+++ b/src/Mie/bhcoat.jl
@@ -22,7 +22,7 @@ Outputs:
 
 References:
 
-- Bohren, C.F., Huffman, D.R., 1983. Absorption and scattering of light by small particles. John Wiley & Sons.
+- C. F. Bohren & D. R. Huffman, *Absorption and Scattering of Light by Small Particles*, Wiley (1983).
 """
 function bhcoat(T, xᵢₙ,
         xₒᵤₜ, mᵢₙ, mₒᵤₜ;
diff --git a/src/Mie/bhmie.jl b/src/Mie/bhmie.jl
index 04a9b95..b1f767c 100644
--- a/src/Mie/bhmie.jl
+++ b/src/Mie/bhmie.jl
@@ -19,7 +19,7 @@ Outputs:
 
 References:
 
-- Bohren, C.F., Huffman, D.R., 1983. Absorption and scattering of light by small particles. John Wiley & Sons.
+- C. F. Bohren & D. R. Huffman, *Absorption and Scattering of Light by Small Particles*, Wiley (1983).
 """
 function bhmie(T, x, m; nₘₐₓ = ceil(Int, Float64(max(x + 4 * ∛x + 2, x * abs(m)))))
     C = complex(T)
diff --git a/src/common/AbstractTransitionMatrix.jl b/src/common/AbstractTransitionMatrix.jl
index d2d735e..4607264 100644
--- a/src/common/AbstractTransitionMatrix.jl
+++ b/src/common/AbstractTransitionMatrix.jl
@@ -32,7 +32,7 @@ rotate(𝐓::AbstractTransitionMatrix{CT, N}, rot::Rotation{3})
 
 Rotate the given T-Matrix `𝐓` by the Euler angle `rot` and generate a new T-Matrix.
 
-For a general T-Matrix, Eq. (5.29) in Mishchenko et al. (2002) is used as a fallback. A `TransitionMatrix` will be returned, which is the most general yet concrete type.
+For a general T-Matrix, Mishchenko et al. (2002), Eq. (5.29), is used as a fallback. A `TransitionMatrix` will be returned, which is the most general yet concrete type.
 
 ```math
 T_{m n m^{\prime} n^{\prime}}^{p p′}(L ; \alpha, \beta, \gamma)=\sum_{m_1=-n}^n \sum_{m_2=-n^{\prime}}^{n^{\prime}} D_{m m_1}^n(\alpha, \beta, \gamma) T_{m_1 n m_2 n^{\prime}}^{p p′}(P) D_{m_2 m^{\prime}}^{n^{\prime}}(-\gamma,-\beta,-\alpha)\quad p,p′=1,2
@@ -95,7 +95,7 @@ Parameters:
 - `φₛ`: the scattering azimuthal angle in radians.
 - `λ`: the wavelength of the incident wave in the host medium. Default to 2π.
 
-For a general T-Matrix, Eq. (5.11) -- Eq. (5.17) in Mishchenko et al. (2002) is used as a fallback.
+For a general T-Matrix, Mishchenko et al. (2002), Eqs. (5.11)–(5.17), is used as a fallback.
 
 ```math
 \begin{array}{l}
@@ -257,7 +257,7 @@ end
 scattering_cross_section(𝐓::AbstractTransitionMatrix{CT, N}, λ=2π) where {CT, N}
 ```
 
-Calculate the scattering cross section per particle averaged over the uniform orientation distribution, according to Eq. (5.140) in Mishchenko et al. (2002).
+Calculate the scattering cross section per particle averaged over the uniform orientation distribution, according to Mishchenko et al. (2002), Eq. (5.140).
 
 ```math
 \left\langle C_{\mathrm{sca}}\right\rangle=\frac{2 \pi}{k_1^2} \sum_{n=1}^{\infty} \sum_{m=-n}^n \sum_{n^{\prime}=1}^{\infty} \sum_{m^{\prime}=-n^{\prime}}^{n^{\prime}} \sum_{k=1}^2 \sum_{l=1}^2\left|T_{m n m^{\prime} n^{\prime}}^{k l}(P)\right|^2
@@ -278,7 +278,7 @@ end
 extinction_cross_section(𝐓::AbstractTransitionMatrix{CT, N}, λ=2π) where {CT, N}
 ```
 
-Calculate the extinction cross section per particle averaged over the uniform orientation distribution, according to Eq. (5.102) in Mishchenko et al. (2002).
+Calculate the extinction cross section per particle averaged over the uniform orientation distribution, according to Mishchenko et al. (2002), Eq. (5.102).
 
 ```math
 \left\langle C_{\mathrm{ext}}\right\rangle=-\frac{2 \pi}{k_1^2} \operatorname{Re} \sum_{n=1}^{\infty} \sum_{m=-n}^n\left[T_{m n n n}^{11}(P)+T_{m n m n}^{22}(P)\right]
@@ -307,7 +307,7 @@ end
 phase_matrix(𝐒::AbstractMatrix)
 ```
 
-Calculate the phase matrix `𝐙` from the amplitude matrix `𝐒`, according to Eq. (2.106) -- Eq. (2.121) in Mishchenko et al. (2002).
+Calculate the phase matrix `𝐙` from the amplitude matrix `𝐒`, according to Mishchenko et al. (2002), Eqs. (2.106)–(2.121).
 """
 function phase_matrix(𝐒::AbstractMatrix)
     𝐙₁₁ = 0.5 * (𝐒[1, 1] * 𝐒[1, 1]' + 𝐒[1, 2] * 𝐒[1, 2]' + 𝐒[2, 1] * 𝐒[2, 1]' +
@@ -406,7 +406,7 @@ end
 expansion_coefficients(𝐓::AbstractTransitionMatrix{CT, N}, λ) where {CT, N}
 ```
 
-Calculate the expansion coefficients from an arbitrary T-Matrix, using Eq. (24) -- (74) in Bi et al. (2014).
+Calculate the expansion coefficients from an arbitrary T-Matrix, using Bi & Yang (2014), Eqs. (24)–(74).
 
 Parameters:
 
@@ -664,7 +664,7 @@ end
 asymmetry_parameter(𝐓, λ)
 ```
 
-Calculate the asymmetry parameter from the given transition matrix, using Eq. (4.92) in Mishchenko et al. (2002):
+Calculate the asymmetry parameter from the given transition matrix, using Mishchenko et al. (2002), Eq. (4.92):
 
 ```math
 \langle\cos\Theta\rangle=\frac{\alpha_1^1}{3}
diff --git a/src/common/AxisymmetricTransitionMatrix.jl b/src/common/AxisymmetricTransitionMatrix.jl
index 4256ee0..4c29a9f 100644
--- a/src/common/AxisymmetricTransitionMatrix.jl
+++ b/src/common/AxisymmetricTransitionMatrix.jl
@@ -1,3 +1,11 @@
+@doc raw"""
+Transition matrix of an **axisymmetric** scatterer (spheroid, cylinder, or
+Chebyshev particle). Axial symmetry makes the matrix block-diagonal in the
+azimuthal order `m`, so only the per-`m` blocks are stored (in the field `𝐓`) —
+a far more compact representation than the general [`TransitionMatrix`](@ref).
+Indexing as `T[m, n, m′, n′, p, p′]` transparently returns the corresponding
+element of the dense layout.
+"""
 struct AxisymmetricTransitionMatrix{CT, N, V <: AbstractVector{<:AbstractMatrix{CT}}, T} <:
        AbstractTransitionMatrix{CT, N}
     𝐓::V
@@ -25,7 +33,7 @@ end
 scattering_cross_section(axi::AxisymmetricTransitionMatrix{CT, N}, λ=2π) where {CT, N}
 ```
 
-Calculate the scattering cross section per particle averaged over the uniform orientation distribution, according to Eq. (5.141) in Mishchenko et al. (2002).
+Calculate the scattering cross section per particle averaged over the uniform orientation distribution, according to Mishchenko et al. (2002), Eq. (5.141).
 
 ```math
 \left\langle C_{\text {sca }}\right\rangle=\frac{2 \pi}{k_1^2} \sum_{n=1}^{\infty} \sum_{n^{\prime}=1}^{\infty} \sum_{m=0}^{\min \left(n, n^{\prime}\right)} \sum_{k=1}^2 \sum_{l=1}^2\left(2-\delta_{m 0}\right)\left|T_{m n m n^{\prime}}^{k l}(P)\right|^2
@@ -72,7 +80,7 @@ end
 extinction_cross_section(axi::AxisymmetricTransitionMatrix{CT, N}, λ=2π) where {CT, N}
 ```
 
-Calculate the extinction cross section per particle averaged over the uniform orientation distribution, according to Eq. (5.107) in Mishchenko et al. (2002).
+Calculate the extinction cross section per particle averaged over the uniform orientation distribution, according to Mishchenko et al. (2002), Eq. (5.107).
 
 ```math
 \left\langle C_{\text {ext }}\right\rangle=-\frac{2 \pi}{k_1^2} \operatorname{Re} \sum_{n=1}^{\infty} \sum_{m=0}^n\left(2-\delta_{m 0}\right)\left[T_{m n m n}^{11}(P)+T_{m n m n}^{22}(P)\right]
diff --git a/src/common/RandomOrientationTransitionMatrix.jl b/src/common/RandomOrientationTransitionMatrix.jl
index 459fa64..4fc58b6 100644
--- a/src/common/RandomOrientationTransitionMatrix.jl
+++ b/src/common/RandomOrientationTransitionMatrix.jl
@@ -8,7 +8,7 @@ end
 RandomOrientationTransitionMatrix(𝐓::AbstractTransitionMatrix{CT, N}) where {CT, N}
 ```
 
-Calculate the random-orientation T-Matrix according to Eq. (5.96) in Mishchenko et al. (2002).
+Calculate the random-orientation T-Matrix according to Mishchenko et al. (2002), Eq. (5.96).
 
 ```math
 \left\langle T_{m n m^{\prime} n^{\prime}}^{k l}(L)\right\rangle=\frac{\delta_{n n^{\prime}} \delta_{m m^{\prime}}}{2 n+1} \sum_{m_1=-n}^n T_{m_1 n m_1 n}^{k l}(P), \quad k, l=1,2
diff --git a/src/solvers.jl b/src/solvers.jl
index 388fed6..2bf0579 100644
--- a/src/solvers.jl
+++ b/src/solvers.jl
@@ -120,12 +120,12 @@ function transition_matrix(s::AbstractAxisymmetricShape{T}, λs::AbstractVector,
 end
 
 # ── Iterative solver ──────────────────────────────────────────────────────────
-"""
+@doc raw"""
     ConvergencePolicy(; threshold, ndgs, maxiter, nₘₐₓ_only, nₛₜₐᵣₜ, Ngₛₜₐᵣₜ, routine_generator)
 
 Convergence settings for an [`Iterative`](@ref) solver. `threshold` is the
 relative `Qsca`/`Qext` tolerance, `ndgs` the quadrature points added per order,
-`nₛₜₐᵣₜ`/`Ngₛₜₐᵣₜ` the starting resolution (`0` ⇒ auto from `k·rₘₐₓ`),
+`nₛₜₐᵣₜ`/`Ngₛₜₐᵣₜ` the starting resolution (`0` ⇒ auto from ``k \cdot r_{\max}``),
 `nₘₐₓ_only` stops once `nₘₐₓ` converges, and `routine_generator` builds the
 stepping routine (default `routine_mishchenko`).
 """
@@ -190,7 +190,7 @@ function _solver_initial(it::Iterative{EBCM}, s, λ)
     it.opts.stable && !(s isa Spheroid) &&
         throw(ArgumentError("Iterative(EBCM; stable=true) is only valid for spheroids: the \
               cancellation-free F⁺ integrands rely on the spheroid surface making \
-              the divergent Laurent terms integrate to zero (Somerville et al. 2013)."))
+              the divergent Laurent terms integrate to zero (Somerville et al. (2013))."))
     pol = it.policy
     nₛₜₐᵣₜ = pol.nₛₜₐᵣₜ
     if nₛₜₐᵣₜ == 0
diff --git a/src/special_functions/wignerd.jl b/src/special_functions/wignerd.jl
index 68b150b..efd9606 100644
--- a/src/special_functions/wignerd.jl
+++ b/src/special_functions/wignerd.jl
@@ -10,7 +10,7 @@ end
 wigner_d([T=Float64,], m::Integer, n::Integer, s::Integer, ϑ::Number) where {T}
 ```
 
-Calculate Wigner (small) d-function ``d_{mn}^s(\theta)`` for a single ``(m, n, s)`` combination, using Eq. (B.1) of Mishchenko et al. (2002).
+Calculate Wigner (small) d-function ``d_{mn}^s(\theta)`` for a single ``(m, n, s)`` combination, using Mishchenko et al. (2002), Eq. (B.1).
 
 ```math
 \begin{aligned}
@@ -61,7 +61,7 @@ end
 wigner_d_recursion([T=Float64,], m::Integer, n::Integer, sₘₐₓ::Integer, ϑ::Number; deriv::Bool = false) where {T}
 ```
 
-Calculate Wigner (small) d-function ``d_{mn}^s(\theta)`` for ``s\in[s_{\min}=\max(|m|, |n|),s_{\max}]`` (alternatively, its derivative as well) via upward recursion, using Eq. (B.22) of Mishchenko et al. (2002).
+Calculate Wigner (small) d-function ``d_{mn}^s(\theta)`` for ``s\in[s_{\min}=\max(|m|, |n|),s_{\max}]`` (alternatively, its derivative as well) via upward recursion, using Mishchenko et al. (2002), Eq. (B.22).
 
 ```math
 \begin{aligned}
@@ -70,7 +70,7 @@ d_{m n}^{s+1}(\vartheta)=& \frac{1}{s \sqrt{(s+1)^{2}-m^{2}} \sqrt{(s+1)^{2}-n^{
 \end{aligned}
 ```
 
-The initial terms are given by Eq. (B.23) and Eq. (B.24).
+The initial terms are given by Mishchenko et al. (2002), Eqs. (B.23) and (B.24).
 
 ```math
 \begin{array}{l}