JuliaRemoteSensing · lucifer1004 · Jun 6, 2026 · Jun 6, 2026 · Jun 6, 2026 · Jun 6, 2026
diff --git a/benchmark/benchmarks.jl b/benchmark/benchmarks.jl
@@ -13,6 +13,7 @@
 #   ebcm               – EBCM blocks, fixed-order, and auto-converging solves
 #   iitm               – Invariant Imbedding T-Matrix (axisymmetric + N-fold)
 #   postprocessing     – far-field observables from a precomputed T-matrix
+#   near_field         – external field reconstruction (VSWF) from a T-matrix
 #   linearization      – analytical Jacobian (baseline for the inv→lu work)
 #   precision          – Float64 vs Double64 on the same EBCM block
 
@@ -127,6 +128,22 @@ let g = SUITE["postprocessing"] = BenchmarkGroup(["farfield"])
         Nγ = 1) seconds=30
 end
 
+# --------------------------------------------------------------------------
+# Near-field — external field reconstruction from a precomputed T-matrix.
+# `scattering_coefficients` is the per-incidence cost (the (p,q)=𝐓(a,b) apply);
+# `scattered_field` is the per-point cost that dominates a dense field grid, so
+# it is benchmarked with (p,q) precomputed (reused across points in practice).
+# --------------------------------------------------------------------------
+let g = SUITE["near_field"] = BenchmarkGroup(["nearfield"])
+    pt = [4.0, 0.0, 3.0]   # outside the spheroid's circumscribing sphere (r > 2)
+    p, q = scattering_coefficients(T_REF, 0.0, 0.0, 1.0 + 0im, 0.0im)
+    g["scattering_coefficients"] = @benchmarkable scattering_coefficients($T_REF, 0.0, 0.0,
+        1.0 + 0im, 0.0im)
+    g["scattered_field_point"] = @benchmarkable scattered_field($p, $q, $λ, $pt)
+    g["total_field_point"] = @benchmarkable total_field($T_REF, $λ, 0.0, 0.0, 1.0 + 0im,
+        0.0im, $pt)
+end
+
 # --------------------------------------------------------------------------
 # Linearization — analytical EBCM Jacobian. This is the direct baseline for
 # the planned `inv(𝐐)` → `lu`/factorization-reuse optimization.

diff --git a/docs/make.jl b/docs/make.jl
@@ -13,7 +13,7 @@ const NB_DIR = normpath(joinpath(@__DIR__, "..", "examples"))
 const NB_OUT = joinpath(@__DIR__, "src", "examples")
 const NOTEBOOKS = ["shapes_gallery.jl", "solver_landscape.jl",
     "angular_scattering.jl", "orientation_averaging.jl", "spectral_sensitivity.jl",
-    "rain_radar.jl"]
+    "rain_radar.jl", "near_field.jl"]
 
 function notebook_title(nb)
     in_markdown = false
@@ -33,10 +33,8 @@ function notebook_title(nb)
 end
 
 const NB_PAGE_PATHS = ["examples/" * replace(nb, ".jl" => ".md") for nb in NOTEBOOKS]
-const NB_PAGES = [
-    notebook_title(nb) => path
-    for (nb, path) in zip(NOTEBOOKS, NB_PAGE_PATHS)
-]
+const NB_PAGES = [notebook_title(nb) => path
+                  for (nb, path) in zip(NOTEBOOKS, NB_PAGE_PATHS)]
 
 function build_examples()
     mkpath(NB_OUT)

diff --git a/docs/src/api.md b/docs/src/api.md
@@ -51,6 +51,18 @@ Pages = ["common/index.jl", "common/AbstractTransitionMatrix.jl",
 Filter = t -> !startswith(string(Base.nameof(t)), "_")
 ```
 
+## Near-field reconstruction
+
+Vector spherical wave functions and the external electromagnetic field
+(incident, scattered, total) reconstructed from any transition matrix (see the
+[Near-field maps from a T-matrix](examples/near_field.md) example).
+
+```@autodocs
+Modules = [TransitionMatrices]
+Pages = ["common/vswf.jl", "common/near_field.jl", "EBCM/near_field.jl"]
+Filter = t -> !startswith(string(Base.nameof(t)), "_")
+```
+
 ## Linearization
 
 The differentiation framework and its backends (see

diff --git a/docs/src/examples/index.md b/docs/src/examples/index.md
@@ -28,6 +28,10 @@ directory of the repository.
   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.
+- [**Near-field maps from a T-matrix**](near_field.md) — reconstruct the external
+  electromagnetic field (incident + scattered) from any T-matrix; a
+  field-enhancement ``|E|/|E_0|`` map and axial cut around a Mie sphere, with the
+  Rayleigh-hypothesis region-of-validity note.
 
 ## Running them locally
 

diff --git a/docs/src/examples/near_field.md b/docs/src/examples/near_field.md
diff --git a/docs/src/usage.md b/docs/src/usage.md
@@ -240,6 +240,91 @@ angles ``\Theta``:
 𝐅 = scattering_matrix(𝐓, 2π, θs)     # one matrix per angle
 ```
 
+## Near-field reconstruction
+
+Beyond the far field, the **electromagnetic field around the particle** can be
+reconstructed from any T-matrix. For an incident plane wave the scattered field
+is expanded in radiating vector spherical wave functions with coefficients
+``(p, q) = \mathbf{T}\,(a, b)``, where ``(a, b)`` are the plane-wave coefficients;
+the total external field is ``\mathbf{E}_\text{inc} + \mathbf{E}_\text{sca}``.
+
+The incident plane wave propagates along ``\hat{\mathbf n} = (\vartheta_\text{inc},
+\varphi_\text{inc})`` with a Jones polarization `(Eθ, Eφ)` in the spherical basis
+``(\hat{\boldsymbol\vartheta}, \hat{\boldsymbol\varphi})`` at that direction. Field
+points are Cartesian (any `AbstractVector` of length 3); the polar axis is `+z`.
+
+```julia
+𝐓 = transition_matrix(spheroid, 2π)
+λ = 2π
+
+# +z propagation, x-polarized (θ̂ at ϑ=0 is x̂):
+ϑ_inc, φ_inc, Eθ, Eφ = 0.0, 0.0, 1.0 + 0im, 0.0im
+r⃗ = [4.0, 0.0, 3.0]                                   # a point outside the particle
+
+E_inc = incident_field(λ, ϑ_inc, φ_inc, Eθ, Eφ, r⃗)   # analytic plane wave
+E_sca = scattered_field(𝐓, λ, ϑ_inc, φ_inc, Eθ, Eφ, r⃗)
+E_tot = total_field(𝐓, λ, ϑ_inc, φ_inc, Eθ, Eφ, r⃗)   # = E_inc + E_sca
+```
+
+When evaluating a **dense grid** of field points, compute the scattered
+coefficients once and reuse them:
+
+```julia
+p, q = scattering_coefficients(𝐓, ϑ_inc, φ_inc, Eθ, Eφ)
+field = [scattered_field(p, q, λ, [x, 0.0, z]) for z in zs, x in xs]
+```
+
+The underlying VSWFs are available directly as [`vswf`](@ref) / [`vswf_cartesian`](@ref).
+
+!!! warning "Region of validity"
+    The radiating (scattered) expansion converges only **outside the smallest
+    sphere circumscribing the particle** (the Rayleigh hypothesis). For a sphere
+    that boundary is the surface; for a non-spherical particle, evaluate only
+    outside its circumscribing sphere.
+
+For the field **inside** a homogeneous sphere, [`internal_field`](@ref)
+reconstructs it from the analytic Mie internal coefficients (at the internal
+wavenumber ``k_\text{int} = m_r k``); it is continuous, tangentially, with the
+external [`total_field`](@ref) across the surface:
+
+```julia
+x, mᵣ = 3.0, 1.5 + 0.05im                        # size parameter and index
+r⃗ = [0.3, 0.2, 0.5]                              # a point inside the sphere
+
+E_in = internal_field(x, mᵣ, λ, ϑ_inc, φ_inc, Eθ, Eφ, r⃗)
+
+# reuse the coefficients across a dense interior grid:
+c, d = internal_coefficients(x, mᵣ, ϑ_inc, φ_inc, Eθ, Eφ)
+E_in = internal_field(c, d, mᵣ, λ, r⃗)
+```
+
+For a general **axisymmetric** particle the internal field is reconstructed from
+the EBCM matrices (``\mathbf{c} = \tfrac{1}{2}\mathbf{Q}^{-1}\mathbf{a}`` per
+azimuthal block), via the shape-based methods:
+
+```julia
+spheroid = Spheroid(1.0, 2.0, 1.4 + 0.02im)
+E_in = internal_field(spheroid, λ, nmax, Ng, ϑ_inc, φ_inc, Eθ, Eφ, r⃗)
+# or, reused across a grid:
+c, d = internal_coefficients(spheroid, λ, nmax, Ng, ϑ_inc, φ_inc, Eθ, Eφ)
+E_in = internal_field(c, d, spheroid.m, λ, r⃗)
+```
+
+!!! note "Region of validity and conditioning"
+    The interior expansion is mathematically guaranteed within the inscribed sphere
+    (radius `rmin(shape)`); empirically it stays stable and physical throughout the
+    interior (it does **not** diverge past the inscribed sphere). The binding limit
+    is instead EBCM ``\mathbf{Q}``-matrix conditioning at high aspect ratio with
+    high `nmax` — since ``\mathbf{c} = \tfrac12\mathbf{Q}^{-1}\mathbf{a}`` inherits
+    ``\mathbf{Q}^{-1}``, the same Float64 cancellation that limits the T-matrix
+    corrupts the coefficients everywhere — so keep `nmax`/aspect where the T-matrix
+    is reliable. The convention is validated against Mie on the degenerate sphere
+    (`Spheroid(R, R, mᵣ)`); absolute accuracy beyond the inscribed sphere for
+    non-spherical shapes is not independently verified.
+
+See the [Near-field maps from a T-matrix](examples/near_field.md) example for a
+field-enhancement map.
+
 ## Orientation averaging
 
 For randomly oriented particles you usually want orientation-averaged

diff --git a/examples/near_field.jl b/examples/near_field.jl
@@ -0,0 +1,181 @@
+### A Pluto.jl notebook ###
+# v1.0.1
+
+using Markdown
+using InteractiveUtils
+
+# ╔═╡ 7a1f0001-0000-4000-8000-000000000001
+begin
+    import Pkg
+    Pkg.activate(@__DIR__)                 # the examples/ environment
+    Pkg.develop(; path = dirname(@__DIR__)) # point TransitionMatrices at this repo
+    Pkg.instantiate()
+    using TransitionMatrices, Plots, LinearAlgebra
+end
+
+# ╔═╡ 7a1f0002-0000-4000-8000-000000000002
+md"""
+# Near-field maps from a T-matrix
+
+`TransitionMatrices.jl` reconstructs the electromagnetic field *around* a particle
+— not just the far-field cross sections — from any computed T-matrix. Given an
+incident plane wave, the scattered field is expanded in radiating vector spherical
+wave functions (VSWFs) with coefficients ``(p, q) = \mathbf{T}\,(a, b)``, where
+``(a, b)`` are the plane-wave expansion coefficients. The total external field is
+
+```math
+\mathbf{E}_\text{tot}(\mathbf r) = \mathbf{E}_\text{inc}(\mathbf r) + \mathbf{E}_\text{sca}(\mathbf r),
+\qquad
+\mathbf{E}_\text{sca}(\mathbf r) = \sum_{n,m} p_{mn}\mathbf{M}_{mn}(k\mathbf r) + q_{mn}\mathbf{N}_{mn}(k\mathbf r).
+```
+
+This works for **any** T-matrix (Mie, EBCM, IITM, Sh-matrix). Here we map the field
+enhancement ``|\mathbf{E}_\text{tot}|/|\mathbf{E}_\text{inc}|`` around a dielectric
+sphere.
+
+!!! note "Region of validity"
+    The radiating expansion converges only **outside the sphere circumscribing the
+    particle** (the Rayleigh hypothesis). For a sphere that boundary *is* the
+    surface, so the map is valid right down to it; for a non-spherical particle,
+    evaluate only outside the circumscribing sphere.
+"""
+
+# ╔═╡ 7a1f0003-0000-4000-8000-000000000003
+md"""
+## 1. A scatterer and its incidence
+
+A dielectric sphere of size parameter ``x = k a = 2.5`` and refractive index
+``m_r = 1.5``, illuminated by an ``x``-polarized plane wave propagating along
+``+z``. With ``\lambda = 2\pi`` the size parameter equals the radius, so ``a = 2.5``.
+
+The incidence direction is ``(\vartheta_\text{inc}, \varphi_\text{inc}) = (0, 0)``;
+the polarization ``(E_\vartheta, E_\varphi) = (1, 0)`` is ``\hat{\mathbf x}`` because
+``\hat{\boldsymbol\vartheta}(0,0) = \hat{\mathbf x}``.
+"""
+
+# ╔═╡ 7a1f0004-0000-4000-8000-000000000004
+begin
+    λ = 2π
+    k = 2π / λ
+    a = 2.5
+    mᵣ = 1.5 + 0.0im
+    x = k * a
+    N = ceil(Int, x + 4 * cbrt(x) + 2)
+    𝐓 = TransitionMatrices.MieTransitionMatrix{ComplexF64, N}(x, mᵣ)
+
+    # Incident plane wave: +z propagation, x-polarized.
+    ϑ_inc, φ_inc = 0.0, 0.0
+    Eθ, Eφ = 1.0 + 0im, 0.0im
+
+    # The scattered coefficients (p, q) = 𝐓 (a, b) depend only on the incidence,
+    # so compute them once and reuse them at every field point.
+    p, q = scattering_coefficients(𝐓, ϑ_inc, φ_inc, Eθ, Eφ)
+    nothing
+end
+
+# ╔═╡ 7a1f0005-0000-4000-8000-000000000005
+md"""
+## 2. Sample the total field on a grid
+
+We evaluate ``|\mathbf{E}_\text{tot}|`` in the ``x``–``z`` plane (the plane of
+incidence and polarization). Points inside the sphere are masked — the external
+expansion is not meant to represent the interior field (that is the *internal*
+field, a separate reconstruction).
+"""
+
+# ╔═╡ 7a1f0006-0000-4000-8000-000000000006
+begin
+    xs = range(-3a, 3a; length = 161)
+    zs = range(-3a, 3a; length = 161)
+    enhancement = fill(NaN, length(zs), length(xs))
+    for (i, zz) in enumerate(zs), (j, xx) in enumerate(xs)
+
+        r = hypot(xx, zz)
+        r ≤ a && continue                      # inside the sphere: masked
+        pos = [xx, 0.0, zz]
+        E_inc = incident_field(λ, ϑ_inc, φ_inc, Eθ, Eφ, pos)
+        E = E_inc + scattered_field(p, q, λ, pos)
+        enhancement[i, j] = norm(E) / norm(E_inc)
+    end
+    nothing
+end
+
+# ╔═╡ 7a1f0007-0000-4000-8000-000000000007
+md"""
+## 3. The field-enhancement map
+
+The incident wave travels upward (``+z``). The map shows the standing-wave
+interference of the incident and scattered fields, the forward-scattering
+concentration on the far (``+z``) side, and the shadow on the near side.
+"""
+
+# ╔═╡ 7a1f0008-0000-4000-8000-000000000008
+let
+    hm = heatmap(xs, zs, enhancement;
+        aspect_ratio = 1, c = :inferno, clims = (0, maximum(filter(isfinite, enhancement))),
+        xlabel = "x", ylabel = "z", colorbar_title = "|E| / |E₀|",
+        title = "dielectric sphere, x = $(round(x; digits = 2)), mᵣ = 1.5",
+        size = (640, 560))
+    # outline the sphere
+    θ = range(0, 2π; length = 200)
+    plot!(hm, a .* cos.(θ), a .* sin.(θ); label = "", lw = 1.5, lc = :cyan)
+end
+
+# ╔═╡ 7a1f0009-0000-4000-8000-000000000009
+md"""
+## 4. A line cut along the optical axis
+
+Along ``x = 0`` the field shows the shadow/illumination asymmetry and the
+forward enhancement directly.
+"""
+
+# ╔═╡ 7a1f000a-0000-4000-8000-00000000000a
+let
+    zline = range(-3a, 3a; length = 400)
+    amp = map(zline) do zz
+        abs(zz) ≤ a && return NaN
+        pos = [0.0, 0.0, zz]
+        E_inc = incident_field(λ, ϑ_inc, φ_inc, Eθ, Eφ, pos)
+        E = E_inc + scattered_field(p, q, λ, pos)
+        norm(E) / norm(E_inc)
+    end
+    plot(zline, amp; lw = 2, xlabel = "z (optical axis)", ylabel = "|E| / |E₀|",
+        label = "", title = "axial cut (x = y = 0)", size = (720, 320))
+    vspan!([-a, a]; alpha = 0.15, c = :gray, label = "sphere")
+end
+
+# ╔═╡ 7a1f000b-0000-4000-8000-00000000000b
+md"""
+## Takeaway
+
+`scattering_coefficients` once, then `scattered_field` / `total_field` at any point
+gives the full external field for any T-matrix — sphere, spheroid, cylinder,
+Chebyshev, or a general particle. Reusing the precomputed `(p, q)` keeps a dense
+field grid cheap.
+
+For the field **inside** the particle, use `internal_field`, which reconstructs
+the regular internal solution from `internal_coefficients`. Those coefficients
+are the VSWF expansion amplitudes for the field inside the particle; compute them
+once and reuse them just like `(p, q)`:
+
+```julia
+c, d = internal_coefficients(x, mᵣ, ϑ_inc, φ_inc, Eθ, Eφ; nmax = N)
+E_inside = internal_field(c, d, mᵣ, λ, [0.2a, 0.0, 0.0])
+```
+
+The Tier-2 implementations live in `src/common/near_field.jl` for homogeneous
+spheres and `src/EBCM/near_field.jl` for axisymmetric particles.
+"""
+
+# ╔═╡ Cell order:
+# ╟─7a1f0002-0000-4000-8000-000000000002
+# ╠═7a1f0001-0000-4000-8000-000000000001
+# ╟─7a1f0003-0000-4000-8000-000000000003
+# ╠═7a1f0004-0000-4000-8000-000000000004
+# ╟─7a1f0005-0000-4000-8000-000000000005
+# ╠═7a1f0006-0000-4000-8000-000000000006
+# ╟─7a1f0007-0000-4000-8000-000000000007
+# ╠═7a1f0008-0000-4000-8000-000000000008
+# ╟─7a1f0009-0000-4000-8000-000000000009
+# ╠═7a1f000a-0000-4000-8000-00000000000a
+# ╟─7a1f000b-0000-4000-8000-00000000000b
diff --git a/examples/orientation_averaging.jl b/examples/orientation_averaging.jl
@@ -58,7 +58,7 @@ begin
         (; method = "analytic (Mishchenko)", Qsca = calc_Csca(Tanalytic, λ),
             Qext = calc_Cext(Tanalytic, λ), g = asymmetry_parameter(Tanalytic, λ)),
         (; method = "numerical (20³ grid)", Qsca = calc_Csca(Tnumeric, λ),
-            Qext = calc_Cext(Tnumeric, λ), g = asymmetry_parameter(Tnumeric, λ)),
+            Qext = calc_Cext(Tnumeric, λ), g = asymmetry_parameter(Tnumeric, λ))
     ]
 end
 
@@ -74,7 +74,8 @@ the analytic closed form gives it directly, at a fraction of the cost.
 begin
     Qref = calc_Csca(Tanalytic, λ)
     Ns = 4:2:18
-    errs = [abs(calc_Csca(orientation_average(T, uniform; Nα = 2, Nβ = N, Nγ = 2), λ) - Qref) / Qref
+    errs = [abs(calc_Csca(orientation_average(T, uniform; Nα = 2, Nβ = N, Nγ = 2), λ) -
+                Qref) / Qref
             for N in Ns]
     nothing
 end

diff --git a/examples/radiative_transfer/generate_fake_wrf.jl b/examples/radiative_transfer/generate_fake_wrf.jl
@@ -77,6 +77,7 @@ function synthetic_fields(nx, ny, nz, nt)
 
     for t in 1:nt
         for j in 1:ny, i in 1:nx
+
             xlat[i, j, t] = 34.0f0 + 0.04f0 * Float32(j - 1)
             xlon[i, j, t] = -98.0f0 + 0.04f0 * Float32(i - 1)
         end