diff --git a/benchmark/benchmarks.jl b/benchmark/benchmarks.jl index c15b428..7f4820a 100644 --- 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 index 56833a6..343b765 100644 --- 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 index b29cf5d..296464e 100644 --- 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 index 2272f79..f014ca8 100644 --- 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 new file mode 100644 index 0000000..5b43312 --- /dev/null +++ b/docs/src/examples/near_field.md @@ -0,0 +1,124 @@ +```@raw html + + + + + +

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

$$\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.

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.

+ +
begin
+    import Pkg
+    Pkg.activate(@__DIR__)                 # the examples/ environment
+    Pkg.develop(; path = dirname(@__DIR__)) # point TransitionMatrices at this repo
+    Pkg.instantiate()
+    using TransitionMatrices, Plots
+end
+ + + +

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}\).

+ +
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
+ + + +

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).

+ +
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 = incident_field(λ, ϑ_inc, φ_inc, Eθ, Eφ, pos) +
+            scattered_field(p, q, λ, pos)
+        enhancement[i, j] = sqrt(abs2(E[1]) + abs2(E[2]) + abs2(E[3]))
+    end
+    nothing
+end
+ + + +

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.

+ +
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
+ + + +

4. A line cut along the optical axis

Along \(x = 0\) the field shows the shadow/illumination asymmetry and the forward enhancement directly.

+ +
let
+    zline = range(-3a, 3a; length = 400)
+    amp = map(zline) do zz
+        abs(zz) ≤ a && return NaN
+        pos = [0.0, 0.0, zz]
+        E = incident_field(λ, ϑ_inc, φ_inc, Eθ, Eφ, pos) + scattered_field(p, q, λ, pos)
+        sqrt(abs2(E[1]) + abs2(E[2]) + abs2(E[3]))
+    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
+ + + +

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, see the internal-field reconstruction (TODO, Tier 2).

+ + +``` + diff --git a/docs/src/usage.md b/docs/src/usage.md index 90e90b2..5a6f6f7 100644 --- 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 new file mode 100644 index 0000000..f52ad3d --- /dev/null +++ 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 index 2261ae2..c59b86e 100644 --- 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 index 357cb65..dd08cdd 100644 --- 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 diff --git a/examples/radiative_transfer/wrf_dualpol_radar.jl b/examples/radiative_transfer/wrf_dualpol_radar.jl index 83729ba..6635ed7 100644 --- a/examples/radiative_transfer/wrf_dualpol_radar.jl +++ b/examples/radiative_transfer/wrf_dualpol_radar.jl @@ -95,8 +95,8 @@ function water_refractive_index(freq_hz, T_celsius) λ_m = C0_M_PER_S / freq_hz ε0 = 8.854e-12 εs = 78.54 * (1 - 4.579e-3 * (T_celsius - 25) + - 1.19e-5 * (T_celsius - 25)^2 - - 2.8e-8 * (T_celsius - 25)^3) + 1.19e-5 * (T_celsius - 25)^2 - + 2.8e-8 * (T_celsius - 25)^3) ε∞ = 5.27137 + 2.16474e-2 * T_celsius - 1.31198e-3 * T_celsius^2 α = -16.8129 / (T_celsius + 273) + 6.09265e-2 λs = 3.3836e-6 * exp(2513.98 / (T_celsius + 273)) @@ -114,9 +114,9 @@ end function rain_axis_ratio(D_mm) D_mm < 0.7 && return 1.0 D_mm < 1.5 && return 1.173 - 0.5265D_mm + 0.4698D_mm^2 - - 0.1317D_mm^3 - 8.5e-3D_mm^4 + 0.1317D_mm^3 - 8.5e-3D_mm^4 1.065 - 6.25e-2D_mm - 3.99e-3D_mm^2 + - 7.66e-4D_mm^3 - 4.095e-5D_mm^4 + 7.66e-4D_mm^3 - 4.095e-5D_mm^4 end function scattering_table(Ds, freq_hz; temperature_celsius = 0.0, @@ -269,11 +269,11 @@ function radar_moments(table, freq_hz, N0, Λ) Zhh = hh * λ^4 / (π^5 * Kw2) Zvv = vv * λ^4 / (π^5 * Kw2) cov = sum(conj(row.Shh_back) * row.Svv_back * weights[i] - for (i, row) in pairs(table)) + for (i, row) in pairs(table)) ρhv = abs(cov) / sqrt(hh * vv) kdp = 180 / π * 1e-3 * λ * sum(real(row.Shh_forward - row.Svv_forward) * weights[i] - for (i, row) in pairs(table)) + for (i, row) in pairs(table)) (; ZH = 10log10(Zhh), ZV = 10log10(Zvv), diff --git a/examples/rain_radar.jl b/examples/rain_radar.jl index e145990..cbdc1ae 100644 --- a/examples/rain_radar.jl +++ b/examples/rain_radar.jl @@ -48,8 +48,8 @@ begin λ_m = c0_m_per_s() / freq_hz ε0 = 8.854e-12 εs = 78.54 * (1 - 4.579e-3 * (T_celsius - 25) + - 1.19e-5 * (T_celsius - 25)^2 - - 2.8e-8 * (T_celsius - 25)^3) + 1.19e-5 * (T_celsius - 25)^2 - + 2.8e-8 * (T_celsius - 25)^3) ε∞ = 5.27137 + 2.16474e-2 * T_celsius - 1.31198e-3 * T_celsius^2 α = -16.8129 / (T_celsius + 273) + 6.09265e-2 λs = 3.3836e-6 * exp(2513.98 / (T_celsius + 273)) @@ -85,16 +85,15 @@ begin function rain_axis_ratio(D_mm) D_mm < 0.7 && return 1.0 D_mm < 1.5 && return 1.173 - 0.5265D_mm + 0.4698D_mm^2 - - 0.1317D_mm^3 - 8.5e-3D_mm^4 + 0.1317D_mm^3 - 8.5e-3D_mm^4 1.065 - 6.25e-2D_mm - 3.99e-3D_mm^2 + - 7.66e-4D_mm^3 - 4.095e-5D_mm^4 + 7.66e-4D_mm^3 - 4.095e-5D_mm^4 end end # ╔═╡ 07e1ca0a-3bf8-45a9-ae4d-8c51971c6f57 begin - round_complex(z; digits = 4) = - complex(round(real(z); digits), round(imag(z); digits)) + round_complex(z; digits = 4) = complex(round(real(z); digits), round(imag(z); digits)) material_summary = [ (; material = "water", freq_GHz = 9.375, T = 0.0, @@ -149,12 +148,10 @@ begin solver = IITM(6, 10, 16) table = scattering_table(Ds, freq_hz; temperature_celsius, solver) - table_summary = [ - (; D = row.D_mm, - axis_ratio = round(row.axis_ratio; digits = 3), - Cext = round(row.Cext_mm2; digits = 4)) - for row in table - ] + table_summary = [(; D = row.D_mm, + axis_ratio = round(row.axis_ratio; digits = 3), + Cext = round(row.Cext_mm2; digits = 4)) + for row in table] end # ╔═╡ aaef0a83-7657-4c40-a875-bc7f18a7b42e @@ -172,8 +169,9 @@ algorithm. # ╔═╡ d78d1d83-508b-4521-ae44-9068d15f4d14 begin - drop_distribution(D_mm, rain_rate_mm_h) = + function drop_distribution(D_mm, rain_rate_mm_h) rain_rate_mm_h == 0 ? 0.0 : 8.0e3 * exp(-4.1 * rain_rate_mm_h^(-0.21) * D_mm) + end function dualpol_moments(table, freq_hz, rain_rate_mm_h; temperature_celsius = 0.0) @@ -181,21 +179,19 @@ begin ΔD = spacing(table) m = water_refractive_index(freq_hz, temperature_celsius) Kw2 = abs2((m^2 - 1) / (m^2 + 2)) - weights = [ - drop_distribution(row.D_mm, rain_rate_mm_h) * ΔD - for row in table - ] + weights = [drop_distribution(row.D_mm, rain_rate_mm_h) * ΔD + for row in table] hh = sum(abs2(row.Shh_back) * weights[i] for (i, row) in pairs(table)) vv = sum(abs2(row.Svv_back) * weights[i] for (i, row) in pairs(table)) Zhh = hh * λ^4 / (π^5 * Kw2) Zvv = vv * λ^4 / (π^5 * Kw2) cov = sum(conj(row.Shh_back) * row.Svv_back * weights[i] - for (i, row) in pairs(table)) + for (i, row) in pairs(table)) ρhv = abs(cov) / sqrt(hh * vv) kdp = 180 / π * 1e-3 * λ * sum(real(row.Shh_forward - row.Svv_forward) * weights[i] - for (i, row) in pairs(table)) + for (i, row) in pairs(table)) (; R = rain_rate_mm_h, ZH = round(10log10(Zhh); digits = 2), @@ -217,14 +213,10 @@ end # ╔═╡ 383e1078-e216-4810-b83d-cd7a6590749b begin rain_rates = [1.0, 10.0, 50.0] - radar = [ - dualpol_moments(table, freq_hz, R; temperature_celsius) - for R in rain_rates - ] - tb = [ - (; R, TB = round(brightness_temperature(table, R); digits = 4)) - for R in rain_rates - ] + radar = [dualpol_moments(table, freq_hz, R; temperature_celsius) + for R in rain_rates] + tb = [(; R, TB = round(brightness_temperature(table, R); digits = 4)) + for R in rain_rates] (; radar, tb) end @@ -246,9 +238,11 @@ begin for f in frequencies_hz local_table = scattering_table(Ds, f; temperature_celsius, solver) for R in rain_rates - push!(rows, (; freq_GHz = round(f / 1e9; digits = 3), R, - TB = round(brightness_temperature(local_table, R; - path_km, Tatm_K); digits = 4))) + push!(rows, + (; freq_GHz = round(f / 1e9; digits = 3), R, + TB = round( + brightness_temperature(local_table, R; + path_km, Tatm_K); digits = 4))) end end rows diff --git a/examples/shapes_gallery.jl b/examples/shapes_gallery.jl index f0924ce..9455054 100644 --- a/examples/shapes_gallery.jl +++ b/examples/shapes_gallery.jl @@ -40,10 +40,10 @@ begin λ = 2π m = 1.5 + 0.02im shapes = ( - spheroid = TransitionMatrices.Spheroid{Float64, ComplexF64}(2.0, 1.0, m), - cylinder = TransitionMatrices.Cylinder{Float64, ComplexF64}(1.0, 2.0, m), + spheroid = TransitionMatrices.Spheroid{Float64, ComplexF64}(2.0, 1.0, m), + cylinder = TransitionMatrices.Cylinder{Float64, ComplexF64}(1.0, 2.0, m), chebyshev = TransitionMatrices.Chebyshev{Float64, ComplexF64}(1.0, 0.1, 4, m), - prism6 = TransitionMatrices.Prism(6, 1.0, 2.0, m), + prism6 = TransitionMatrices.Prism(6, 1.0, 2.0, m) ) end @@ -57,8 +57,10 @@ IITM solve. Pluto renders the resulting vector of named tuples as a table. # ╔═╡ 2173874a-6082-11f1-8f4b-f3cee77b5fae begin - Tmatrix(s) = s isa TransitionMatrices.Prism ? + function Tmatrix(s) + s isa TransitionMatrices.Prism ? transition_matrix(s, λ, IITM(8, 12, 24, 24)) : calc_T(s, λ) + end summary = map(collect(pairs(shapes))) do (name, s) T = Tmatrix(s) diff --git a/examples/solver_landscape.jl b/examples/solver_landscape.jl index cfc8266..cb2933d 100644 --- a/examples/solver_landscape.jl +++ b/examples/solver_landscape.jl @@ -39,12 +39,16 @@ Three routes to the same cross section: # ╔═╡ 246053c2-6082-11f1-9425-cd50e8e63bdd begin λ = 2π - prolate = TransitionMatrices.Spheroid{Float64, ComplexF64}(2.5198421, 10.079368, 1.55 + 0.01im) + prolate = TransitionMatrices.Spheroid{Float64, ComplexF64}(2.5198421, 10.079368, 1.55 + + 0.01im) routes = [ - (; route = "EBCM fixed (n=24)", Csca = calc_Csca(transition_matrix(prolate, λ, EBCM(24, 384)), λ)), - (; route = "Iterative(EBCM; stable)", Csca = calc_Csca(transition_matrix(prolate, λ, Iterative(EBCM; stable = true)), λ)), - (; route = "IITM", Csca = calc_Csca(transition_matrix(prolate, λ, IITM(24, 30, 60)), λ)), + (; route = "EBCM fixed (n=24)", + Csca = calc_Csca(transition_matrix(prolate, λ, EBCM(24, 384)), λ)), + (; route = "Iterative(EBCM; stable)", + Csca = calc_Csca(transition_matrix(prolate, λ, Iterative(EBCM; stable = true)), λ)), + (; route = "IITM", + Csca = calc_Csca(transition_matrix(prolate, λ, IITM(24, 30, 60)), λ)) ] end @@ -60,8 +64,10 @@ classic EBCM vs the `stable = true` EBCM. Watch the classic curve blow up. begin ref = calc_Csca(transition_matrix(prolate, λ, EBCM(40, 600; stable = true)), λ) ns = 16:2:40 - classic_err = [abs(calc_Csca(transition_matrix(prolate, λ, EBCM(n, 16n)), λ) - ref) / ref for n in ns] - stable_err = [abs(calc_Csca(transition_matrix(prolate, λ, EBCM(n, 16n; stable = true)), λ) - ref) / ref for n in ns] + classic_err = [abs(calc_Csca(transition_matrix(prolate, λ, EBCM(n, 16n)), λ) - ref) / + ref for n in ns] + stable_err = [abs(calc_Csca(transition_matrix(prolate, λ, EBCM(n, 16n; stable = true)), λ) - + ref) / ref for n in ns] nothing end diff --git a/examples/spectral_sensitivity.jl b/examples/spectral_sensitivity.jl index 2bd8d28..ab494ad 100644 --- a/examples/spectral_sensitivity.jl +++ b/examples/spectral_sensitivity.jl @@ -116,7 +116,8 @@ begin f_sh(λ) = ForwardDiff.derivative(l -> Csca(l, mᵣ), λ) f_classic(λ) = ForwardDiff.derivative(l -> Csca_classic(l, mᵣ), λ) - f_sh(λs[1]); f_classic(λs[1]) # warm up + f_sh(λs[1]); + f_classic(λs[1]) # warm up t_prepare = @belapsed prepare_sh($spheroid, $nmax, $Ng) samples=1 evals=1 t_sh = @belapsed [f_sh(λ) for λ in $λs] samples=1 evals=1 diff --git a/src/EBCM/index.jl b/src/EBCM/index.jl index 2fc6449..ed3e1f3 100644 --- a/src/EBCM/index.jl +++ b/src/EBCM/index.jl @@ -3,3 +3,4 @@ include("stabilization.jl") include("shmatrix.jl") include("routines.jl") include("linearization.jl") +include("near_field.jl") diff --git a/src/EBCM/near_field.jl b/src/EBCM/near_field.jl new file mode 100644 index 0000000..579e811 --- /dev/null +++ b/src/EBCM/near_field.jl @@ -0,0 +1,147 @@ +# ── Internal field (Tier 2): general axisymmetric particle (EBCM) ───────────── +# +# Extends `internal_coefficients` / `internal_field` (defined for the homogeneous +# sphere in `common/near_field.jl`) to a general axisymmetric particle. The +# interior coefficients come from the EBCM matrices rather than an analytic +# formula. With the package's 𝐐 = 𝐏 + i𝐔 (so 𝐓 = -𝐏 𝐐⁻¹) the Waterman relation +# gives the internal coefficients as 𝐐⁻¹ 𝐚 per azimuthal block; an overall factor +# ½ — universal across size and refractive index, fixed by matching the analytic +# Mie internal field on a sphere to machine precision — converts to this package's +# VSWF normalization: +# +# 𝐜⁽ᵐ⁾ = ½ 𝐐⁽ᵐ⁾⁻¹ 𝐚⁽ᵐ⁾ (degree +m). +# +# The −m block follows the same ±m symmetry the T-matrix uses (see +# `AxisymmetricTransitionMatrix`): the M↔N cross blocks flip sign, so +# 𝐜⁽⁻ᵐ⁾ = 𝐃 (½ 𝐐⁽ᵐ⁾⁻¹) 𝐃 𝐚⁽⁻ᵐ⁾ with 𝐃 = diag(+𝐈_M, -𝐈_N). Each block stacks the +# M (p=1) coefficients above the N (p=2) coefficients over n = max(1,m):nmax. +# +# VALIDITY. The regular-VSWF interior expansion is mathematically guaranteed to +# converge to the true field at least within the inscribed sphere (radius +# `rmin(shape)`). Empirically (tested) it does NOT diverge outside it: the +# reconstruction stays numerically stable and physical throughout the interior — +# tips included, well beyond the inscribed sphere — as long as the underlying EBCM +# build is well conditioned. The real failure mode is EBCM 𝐐-matrix +# ill-conditioning at high aspect ratio with high nmax (the same Float64 +# cancellation the F⁺ `stable` path fixes for the T-matrix): since 𝐜 = ½ 𝐐⁻¹ 𝐚 +# inherits 𝐐⁻¹, a poorly-conditioned 𝐐 corrupts the coefficients EVERYWHERE +# (including the core), not just near the surface. Keep aspect×nmax in the regime +# where the T-matrix itself is reliable. +# +# Absolute accuracy beyond the inscribed sphere is not independently validated for +# non-spherical shapes (there is no self-contained surface check — tangential +# boundary continuity would need both interior and exterior expansions to converge +# at the surface, which the Rayleigh hypothesis forbids — and no built-in reference +# method). The convention (the ½ factor and the ±m symmetry) is pinned to machine +# precision by the degenerate-sphere (a = c) comparison against the analytic Mie +# internal field. + +@doc raw""" +``` +internal_coefficients(shape::AbstractAxisymmetricShape, λ, nmax, Ng, ϑ_inc, φ_inc, Eθ, Eφ) -> (c, d) +``` + +Internal-field expansion coefficients `(cₘₙ, dₘₙ)` for a general axisymmetric +particle, obtained from the EBCM matrices: per azimuthal block ``\mathbf{c} = +\tfrac{1}{2}\mathbf{Q}^{-1}\mathbf{a}`` with ``\mathbf{Q} = \mathbf{P} + +\mathrm{i}\mathbf{U}``, under an incident plane wave (`(ϑ_inc, φ_inc)` propagation, +Jones polarization `(Eθ, Eφ)`). Returns `OffsetArray`s indexed `[m, n]`; reuse them +across field points with [`internal_field`](@ref). + +The interior expansion is guaranteed within the inscribed sphere (radius +`rmin(shape)`) and is empirically stable throughout the interior; the practical +limit is EBCM conditioning at high aspect ratio with high `nmax` (keep +`nmax`/aspect where the T-matrix itself is reliable). See [`internal_field`](@ref). +""" +function internal_coefficients(shape::AbstractAxisymmetricShape{T, CT}, λ, nmax::Integer, + Ng::Integer, ϑ_inc, φ_inc, Eθ, Eφ) where {T, CT} + Tr = real(CT) + a, b = _plane_wave_coefficients(nmax, Tr(ϑ_inc), Tr(φ_inc), CT(Eθ), CT(Eφ)) + c = OffsetArray(zeros(CT, 2nmax + 1, nmax), (-nmax):nmax, 1:nmax) + d = OffsetArray(zeros(CT, 2nmax + 1, nmax), (-nmax):nmax, 1:nmax) + for mₐ in 0:nmax + 𝐏, + 𝐔 = mₐ == 0 ? ebcm_matrices_m₀(shape, λ, nmax, Ng) : + ebcm_matrices_m(mₐ, shape, λ, nmax, Ng) + 𝐐 = 𝐏 .+ im .* 𝐔 + ns = max(1, mₐ):nmax + nn = length(ns) + avec = CT[[a[mₐ, n] for n in ns]; [b[mₐ, n] for n in ns]] + cvec = (𝐐 \ avec) ./ 2 + for (i, n) in enumerate(ns) + c[mₐ, n] = cvec[i] + d[mₐ, n] = cvec[nn + i] + end + if mₐ > 0 + 𝐃 = CT[ones(CT, nn); -ones(CT, nn)] + amvec = CT[[a[-mₐ, n] for n in ns]; [b[-mₐ, n] for n in ns]] + cmvec = (𝐃 .* (𝐐 \ (𝐃 .* amvec))) ./ 2 + for (i, n) in enumerate(ns) + c[-mₐ, n] = cmvec[i] + d[-mₐ, n] = cmvec[nn + i] + end + end + end + return c, d +end + +@doc raw""" +``` +internal_field(shape::AbstractAxisymmetricShape, λ, nmax, Ng, ϑ_inc, φ_inc, Eθ, Eφ, r⃗) -> SVector{3} +``` + +Internal electric field inside a general axisymmetric particle at the Cartesian +point `r⃗`, via the EBCM internal coefficients (see [`internal_coefficients`](@ref)). +The relative refractive index is taken from `shape`, and `r⃗` should lie inside the +particle. Convergence is guaranteed within the inscribed sphere (radius +`rmin(shape)`) and is empirically stable throughout the interior; the binding +constraint is EBCM conditioning at high aspect × `nmax`, not the inscribed-sphere +radius. +""" +function internal_field(shape::AbstractAxisymmetricShape, λ, nmax::Integer, Ng::Integer, + ϑ_inc, φ_inc, Eθ, Eφ, r⃗) + c, d = internal_coefficients(shape, λ, nmax, Ng, ϑ_inc, φ_inc, Eθ, Eφ) + return internal_field(c, d, shape.m, λ, r⃗) +end + +@testitem "EBCM internal field matches Mie on a sphere (degenerate spheroid)" begin + using TransitionMatrices: Spheroid, internal_field + + # A sphere is a spheroid with a = c. Its EBCM internal field must reproduce the + # analytic Mie internal field — this locks the EBCM internal-coefficient + # convention (the ½ factor and the ±m symmetry) to machine precision. + λ = 2π + k0 = 2π / λ + @testset "R=$R, mᵣ=$mᵣ" for (R, mᵣ) in ((2.0, 1.5 + 0.05im), (1.5, 1.33 + 0.0im)) + x = k0 * R + nmax = ceil(Int, x + 4cbrt(x) + 2) + Ng = 200 + sph = Spheroid{Float64, ComplexF64}(R, R, mᵣ) + worst = 0.0 + for ϑ in (0.4, 1.1, 2.0), φ in (0.0, 1.3, 4.0), rr in (0.3R, 0.6R, 0.85R) + sϑ, cϑ, sφ, cφ = sin(ϑ), cos(ϑ), sin(φ), cos(φ) + pos = [rr * sϑ * cφ, rr * sϑ * sφ, rr * cϑ] + Eebcm = internal_field(sph, λ, nmax, Ng, 0.0, 0.0, 1.0 + 0im, 0.0im, pos) + Emie = internal_field(x, mᵣ, λ, 0.0, 0.0, 1.0 + 0im, 0.0im, pos; nmax) + worst = max(worst, maximum(abs, Eebcm - Emie) / maximum(abs, Emie)) + end + @test worst < 1e-8 + end +end + +@testitem "EBCM internal field reconstructs inside a spheroid (runs, finite)" begin + using TransitionMatrices: Spheroid, internal_coefficients, internal_field + + # Non-spherical: no self-contained reference, but the coefficients/field must be + # finite and the two-step and one-shot forms must agree, inside the inscribed + # sphere (radius = min semi-axis). + λ = 2π + sph = Spheroid{Float64, ComplexF64}(1.0, 2.0, 1.4 + 0.02im) # prolate, rmin = 1 + nmax, Ng = 10, 200 + c, d = internal_coefficients(sph, λ, nmax, Ng, 0.3, 0.5, 1.0 + 0im, 0.0im) + pos = [0.4, 0.2, 0.5] # |pos| ≈ 0.67 < inscribed radius 1 + E1 = internal_field(c, d, sph.m, λ, pos) + E2 = internal_field(sph, λ, nmax, Ng, 0.3, 0.5, 1.0 + 0im, 0.0im, pos) + @test all(isfinite, abs.(E1)) + @test maximum(abs, E1 - E2) < 1e-12 +end diff --git a/src/TransitionMatrices.jl b/src/TransitionMatrices.jl index 6af0e58..11fc0c7 100644 --- a/src/TransitionMatrices.jl +++ b/src/TransitionMatrices.jl @@ -40,6 +40,10 @@ export AbstractLinearizationBackend, MieLinearization, EBCMLinearization, IITMLi UnsupportedLinearization, variables, rebuild, derivative, supports_linearization, linearize_transition_matrix, linearize_observable +# Near-field reconstruction +export vswf, vswf_cartesian, scattering_coefficients, incident_field, scattered_field, + total_field, internal_coefficients, internal_field + # Utility functions export OrderDegreeIterator, rotate, amplitude_matrix, phase_matrix, scattering_matrix, expansion_coefficients, diff --git a/src/common/index.jl b/src/common/index.jl index ba5adc2..33e9529 100644 --- a/src/common/index.jl +++ b/src/common/index.jl @@ -1,4 +1,6 @@ include("AbstractTransitionMatrix.jl") include("AxisymmetricTransitionMatrix.jl") include("RandomOrientationTransitionMatrix.jl") +include("vswf.jl") +include("near_field.jl") include("utils.jl") diff --git a/src/common/near_field.jl b/src/common/near_field.jl new file mode 100644 index 0000000..683450d --- /dev/null +++ b/src/common/near_field.jl @@ -0,0 +1,477 @@ +# ── Near-field reconstruction (Tier 1: external field) ──────────────────────── +# +# Given a T-matrix `𝐓` and an incident plane wave, reconstruct the electric field +# at points OUTSIDE the circumscribing sphere of the particle: the incident field +# (analytic plane wave), the scattered field, and their sum (the total field). +# +# The incident plane wave propagates along the direction `n̂ = (ϑ_inc, φ_inc)` and +# carries a polarization given as a Jones vector `(Eθ, Eφ)` in the spherical basis +# `(𝛝̂, 𝛗̂)` evaluated at `n̂`: +# +# 𝐄_inc(𝐫) = (Eθ 𝛝̂ + Eφ 𝛗̂) exp(i k n̂·𝐫). +# +# It is expanded in regular VSWFs (`vswf(:regular, …)`) +# 𝐄_inc(𝐫) = Σ_{n,m} [ aₘₙ RgMₘₙ(k𝐫) + bₘₙ RgNₘₙ(k𝐫) ], +# with coefficients (locked numerically against the analytic plane wave to machine +# precision, Mishchenko convention; γₘₙ as in `vswf`): +# aₘₙ = 4π iⁿ γₘₙ e^{-imφ_inc} ( −i πₘₙ(ϑ_inc) Eθ − τₘₙ(ϑ_inc) Eφ ), +# bₘₙ = 4π iⁿ⁻¹ γₘₙ e^{-imφ_inc} ( τₘₙ(ϑ_inc) Eθ − i πₘₙ(ϑ_inc) Eφ ). +# +# The scattered field uses the radiating VSWFs and the scattered coefficients +# (p,q) = 𝐓 (a,b): +# 𝐄_sca(𝐫) = Σ_{n,m} [ pₘₙ Mₘₙ(k𝐫) + qₘₙ Nₘₙ(k𝐫) ], +# pₘₙ = Σ_{n′m′} 𝐓₁₁ aₘ′ₙ′ + 𝐓₁₂ bₘ′ₙ′, qₘₙ = Σ_{n′m′} 𝐓₂₁ aₘ′ₙ′ + 𝐓₂₂ bₘ′ₙ′. +# +# Validity: the radiating expansion converges only OUTSIDE the smallest sphere +# circumscribing the particle (the Rayleigh hypothesis). Inside that sphere the +# scattered series may diverge; callers are responsible for evaluating at r larger +# than the circumscribing radius. + +@inline _dot3(a, b) = a[1] * b[1] + a[2] * b[2] + a[3] * b[3] + +@inline function _cart_to_sph(r⃗) + x, y, z = r⃗[1], r⃗[2], r⃗[3] + r = sqrt(x^2 + y^2 + z^2) + ϑ = acos(clamp(z / r, -one(r), one(r))) + φ = atan(y, x) + return r, ϑ, φ +end + +# Plane-wave expansion coefficients aₘₙ, bₘₙ (degree m, order n), as OffsetArrays +# indexed `[m, n]`, m ∈ -N:N, n ∈ 1:N. +function _plane_wave_coefficients(N::Integer, ϑ::T, φ::T, Eθ::Complex{T}, + Eφ::Complex{T}) where {T} + a = OffsetArray(zeros(Complex{T}, 2N + 1, N), (-N):N, 1:N) + b = OffsetArray(zeros(Complex{T}, 2N + 1, N), (-N):N, 1:N) + for m in (-N):N + dvec, τvec = wigner_d_recursion(T, 0, m, N, ϑ; deriv = true) + for n in max(1, abs(m)):N + P = dvec[n] + τ = τvec[n] + π_ = pi_func(T, m, n, ϑ; d = P) + γ = (-1)^(m & 1) * sqrt(T(2n + 1) / (4 * T(π) * n * (n + 1))) + pref = 4 * T(π) * cis(-m * φ) * γ + a[m, n] = pref * im^(n & 3) * (-im * π_ * Eθ - τ * Eφ) + b[m, n] = pref * im^((n - 1) & 3) * (τ * Eθ - im * π_ * Eφ) + end + end + return a, b +end + +@doc raw""" +``` +scattering_coefficients(𝐓, ϑ_inc, φ_inc, Eθ, Eφ) -> (p, q) +``` + +Scattered-field expansion coefficients `(pₘₙ, qₘₙ)` for the T-matrix `𝐓` under an +incident plane wave propagating along `(ϑ_inc, φ_inc)` with Jones polarization +`(Eθ, Eφ)` in the spherical basis at the incidence direction. Returns two +`OffsetArray`s indexed `[m, n]` (m ∈ -N:N, n ∈ 1:N), where `(p, q) = 𝐓 (a, b)` and +`(a, b)` are the incident plane-wave coefficients. Reuse these across many field +points (see [`scattered_field`](@ref)). +""" +function scattering_coefficients(𝐓::AbstractTransitionMatrix{CT, N}, ϑ_inc, φ_inc, + Eθ, Eφ) where {CT, N} + T = real(CT) + a, b = _plane_wave_coefficients(N, T(ϑ_inc), T(φ_inc), CT(Eθ), CT(Eφ)) + p = OffsetArray(zeros(CT, 2N + 1, N), (-N):N, 1:N) + q = OffsetArray(zeros(CT, 2N + 1, N), (-N):N, 1:N) + for n in 1:N, m in (-n):n + + s1 = zero(CT) + s2 = zero(CT) + for n′ in 1:N, m′ in (-n′):n′ + + aₘ′ₙ′ = a[m′, n′] + bₘ′ₙ′ = b[m′, n′] + s1 += 𝐓[m, n, m′, n′, 1, 1] * aₘ′ₙ′ + 𝐓[m, n, m′, n′, 1, 2] * bₘ′ₙ′ + s2 += 𝐓[m, n, m′, n′, 2, 1] * aₘ′ₙ′ + 𝐓[m, n, m′, n′, 2, 2] * bₘ′ₙ′ + end + p[m, n] = s1 + q[m, n] = s2 + end + return p, q +end + +# Sum Σ c₁ₘₙ Fₘₙ + c₂ₘₙ Gₘₙ at one point, where (F,G)=(RgM,RgN) for kind=:regular +# or (M,N) for kind=:outgoing. Radial factors computed once; angular per m. The +# wavenumber `k` may be COMPLEX (the internal field uses k_int = mᵣ·k); the +# geometry (r, ϑ, φ) and the angular functions stay real, only the radial +# argument x = k·r and the resulting field become complex. +function _field_from_coeffs(c₁, c₂, kind::Symbol, k::Number, r⃗, N::Integer) + T = float(real(promote_type(typeof(k), eltype(r⃗)))) + r⃗T = SVector{3, T}(r⃗[1], r⃗[2], r⃗[3]) + r = sqrt(r⃗T[1]^2 + r⃗T[2]^2 + r⃗T[3]^2) + x = k * r + if _is_regular_origin(kind, x, T) + CT = complex(T) + E = zero(SVector{3, CT}) + N ≥ 1 || return E + for m in -1:1 + _, 𝐆 = _regular_vswf_origin_limit(T, m, 1) + E += c₂[m, 1] * 𝐆 + end + return E + end + _, ϑ, φ = _cart_to_sph(r⃗T) + zₙ, xzₙ′ = _vswf_radial(kind, N, x) + r̂, ϑ̂, φ̂ = _sph_unit_vectors(ϑ, φ) + + CT = complex(T) + E = zero(SVector{3, CT}) + for m in (-N):N + dvec, τvec = wigner_d_recursion(T, 0, m, N, ϑ; deriv = true) + eⁱᵐᵠ = cis(m * φ) + for n in max(1, abs(m)):N + P = dvec[n] + τ = τvec[n] + π_ = pi_func(T, m, n, ϑ; d = P) + γ = (-1)^(m & 1) * sqrt(T(2n + 1) / (4 * T(π) * n * (n + 1))) + pref = γ * eⁱᵐᵠ + zn = zₙ[n] + dzn = xzₙ′[n] + znx = zn / x + 𝐅 = pref * zn * (im * π_ .* ϑ̂ .- τ .* φ̂) + 𝐆 = pref * + (n * (n + 1) * znx * P .* r̂ .+ dzn .* (τ .* ϑ̂ .+ im * π_ .* φ̂)) + E += c₁[m, n] * 𝐅 + c₂[m, n] * 𝐆 + end + end + return E +end + +@doc raw""" +``` +incident_field(λ, ϑ_inc, φ_inc, Eθ, Eφ, r⃗) -> SVector{3} +``` + +The incident plane wave ``(E_\theta\hat{\boldsymbol\vartheta} + E_\varphi\hat{\boldsymbol\varphi})\exp(\mathrm{i}k\,\hat{\mathbf n}\cdot\mathbf r)`` +evaluated at the Cartesian point `r⃗`, propagating along `n̂ = (ϑ_inc, φ_inc)` with +`k = 2π/λ`. Polarization `(Eθ, Eφ)` is given in the spherical basis at `n̂`. +""" +function incident_field(λ, ϑ_inc, φ_inc, Eθ, Eφ, r⃗) + T = float(promote_type(typeof(λ), typeof(ϑ_inc), typeof(φ_inc), eltype(r⃗))) + k = 2 * T(π) / T(λ) + sϑ, cϑ = sincos(T(ϑ_inc)) + sφ, cφ = sincos(T(φ_inc)) + n̂ = SVector(sϑ * cφ, sϑ * sφ, cϑ) + ϑ̂ = SVector(cϑ * cφ, cϑ * sφ, -sϑ) + φ̂ = SVector(-sφ, cφ, zero(T)) + return (Complex{T}(Eθ) .* ϑ̂ .+ Complex{T}(Eφ) .* φ̂) .* cis(k * _dot3(n̂, r⃗)) +end + +@doc raw""" +``` +scattered_field(p, q, λ, r⃗) -> SVector{3} +scattered_field(𝐓, λ, ϑ_inc, φ_inc, Eθ, Eφ, r⃗) -> SVector{3} +``` + +Scattered electric field at the Cartesian point `r⃗`, reconstructed from the +radiating VSWFs. The first form reuses precomputed coefficients `(p, q)` from +[`scattering_coefficients`](@ref) (efficient for evaluating many points); the +second computes them on the fly for the T-matrix `𝐓` and incidence. + +The point `r⃗` must lie OUTSIDE the sphere circumscribing the particle (the +radiating series diverges inside it). +""" +function scattered_field(p::OffsetArray, q::OffsetArray, λ, r⃗) + N = size(p, 2) + Tr = real(eltype(p)) + k = 2 * Tr(π) / λ + return _field_from_coeffs(p, q, :outgoing, k, r⃗, N) +end + +function scattered_field(𝐓::AbstractTransitionMatrix, λ, ϑ_inc, φ_inc, Eθ, Eφ, r⃗) + p, q = scattering_coefficients(𝐓, ϑ_inc, φ_inc, Eθ, Eφ) + return scattered_field(p, q, λ, r⃗) +end + +@doc raw""" +``` +total_field(𝐓, λ, ϑ_inc, φ_inc, Eθ, Eφ, r⃗) -> SVector{3} +``` + +Total external field `𝐄_inc(r⃗) + 𝐄_sca(r⃗)` at the Cartesian point `r⃗` for the +T-matrix `𝐓` under the incident plane wave (see [`incident_field`](@ref) and +[`scattered_field`](@ref)). +""" +function total_field(𝐓::AbstractTransitionMatrix, λ, ϑ_inc, φ_inc, Eθ, Eφ, r⃗) + return incident_field(λ, ϑ_inc, φ_inc, Eθ, Eφ, r⃗) + + scattered_field(𝐓, λ, ϑ_inc, φ_inc, Eθ, Eφ, r⃗) +end + +# ── Internal field (Tier 2): homogeneous sphere (Mie) ───────────────────────── +# +# The field INSIDE the particle is expanded in regular VSWFs at the internal +# wavenumber k_int = mᵣ·k: +# 𝐄_int(𝐫) = Σ_{n,m} c_{mn} RgM_{mn}(k_int 𝐫) + d_{mn} RgN_{mn}(k_int 𝐫). +# Unlike the scattered field, the internal coefficients are NOT determined by the +# T-matrix alone (which maps incident → scattered); they need the method's +# interior solution. For a homogeneous sphere they are the analytic Mie internal +# coefficients (Bohren & Huffman 1983, Eqs. 4.52–4.53), per order n: +# cₙ = mᵣ(ψₙ(x)ξₙ′(x) − ξₙ(x)ψₙ′(x)) / (ψₙ(mᵣx)ξₙ′(x) − mᵣ ξₙ(x)ψₙ′(mᵣx)) (M / magnetic) +# dₙ = mᵣ(ψₙ(x)ξₙ′(x) − ξₙ(x)ψₙ′(x)) / (mᵣ ψₙ(mᵣx)ξₙ′(x) − ξₙ(x)ψₙ′(mᵣx)) (N / electric) +# with ξₙ = ψₙ + iχₙ the package's outgoing Riccati–Bessel. The internal degree-m +# coefficients then ride on the incident ones, cₘₙ = cₙ aₘₙ, dₘₙ = dₙ bₘₙ. +# Validated by Maxwell boundary continuity (tangential 𝐄 and 𝐇 continuous across +# the surface) against the already-validated external field. EBCM/IITM interior +# fields (needing the retained 𝐐 matrix) are a separate, deferred extension. +function _mie_internal_coefficients(x::Real, mᵣ::Number, nmax::Integer) + T = float(real(typeof(float(x)))) + CT = complex(T) + xT = T(x) + m = CT(mᵣ) + ψx, ψx′ = ricattibesselj(nmax, estimate_ricattibesselj_extra_terms(nmax, xT), xT) + χx, χx′ = ricattibessely(nmax, xT) + mx = m * xT + ψm, ψm′ = ricattibesselj(nmax, estimate_ricattibesselj_extra_terms(nmax, mx), mx) + c = Vector{CT}(undef, nmax) + d = Vector{CT}(undef, nmax) + for n in 1:nmax + ξ = ψx[n] + im * χx[n] + ξ′ = ψx′[n] + im * χx′[n] + Dₐ = m * ψm[n] * ξ′ - ξ * ψm′[n] + D_b = ψm[n] * ξ′ - m * ξ * ψm′[n] + num = m * (ψx[n] * ξ′ - ξ * ψx′[n]) + c[n] = num / D_b + d[n] = num / Dₐ + end + return c, d +end + +@doc raw""" +``` +internal_coefficients(x, mᵣ, ϑ_inc, φ_inc, Eθ, Eφ; nmax) -> (c, d) +``` + +Internal-field expansion coefficients `(cₘₙ, dₘₙ)` for a **homogeneous sphere** of +size parameter `x = k a` and relative refractive index `mᵣ`, under an incident +plane wave (`(ϑ_inc, φ_inc)` propagation, Jones polarization `(Eθ, Eφ)` in the +spherical basis at the incidence direction). Returns two `OffsetArray`s indexed +`[m, n]`; reuse them across many field points with [`internal_field`](@ref). +""" +function internal_coefficients(x::Real, mᵣ::Number, ϑ_inc, φ_inc, Eθ, Eφ; + nmax::Integer = ceil(Int, x + 4cbrt(x) + 2)) + T = float(real(promote_type(typeof(float(x)), typeof(float(ϑ_inc))))) + CT = complex(T) + cₙ, dₙ = _mie_internal_coefficients(x, mᵣ, nmax) + a, b = _plane_wave_coefficients(nmax, T(ϑ_inc), T(φ_inc), CT(Eθ), CT(Eφ)) + c = OffsetArray(zeros(CT, 2nmax + 1, nmax), (-nmax):nmax, 1:nmax) + d = OffsetArray(zeros(CT, 2nmax + 1, nmax), (-nmax):nmax, 1:nmax) + for n in 1:nmax, m in (-n):n + + c[m, n] = cₙ[n] * a[m, n] + d[m, n] = dₙ[n] * b[m, n] + end + return c, d +end + +@doc raw""" +``` +internal_field(c, d, mᵣ, λ, r⃗) -> SVector{3} +internal_field(x, mᵣ, λ, ϑ_inc, φ_inc, Eθ, Eφ, r⃗; nmax) -> SVector{3} +``` + +Internal electric field (inside a homogeneous sphere) at the Cartesian point `r⃗`, +reconstructed from regular VSWFs at the internal wavenumber ``k_\text{int} = m_r k``. +The first form reuses precomputed coefficients `(c, d)` from +[`internal_coefficients`](@ref) (efficient for a dense grid); the second computes +them on the fly for size parameter `x = k a` and refractive index `mᵣ`. + +`r⃗` must lie INSIDE the sphere (the regular expansion represents the interior +field). The external field is [`total_field`](@ref); the two are continuous +(tangentially) across the surface. +""" +function internal_field(c::OffsetArray, d::OffsetArray, mᵣ, λ, r⃗) + N = size(c, 2) + Tr = real(eltype(c)) + k = 2 * Tr(π) / λ + return _field_from_coeffs(c, d, :regular, mᵣ * k, r⃗, N) +end + +function internal_field(x::Real, mᵣ, λ, ϑ_inc, φ_inc, Eθ, Eφ, r⃗; + nmax::Integer = ceil(Int, x + 4cbrt(x) + 2)) + c, d = internal_coefficients(x, mᵣ, ϑ_inc, φ_inc, Eθ, Eφ; nmax) + return internal_field(c, d, mᵣ, λ, r⃗) +end + +# The general axisymmetric (EBCM) internal field — extending `internal_coefficients` +# / `internal_field` to non-spherical shapes — lives in `src/EBCM/near_field.jl`, +# included after the EBCM matrices and the shape types it depends on. + +@testitem "Internal field is tangentially continuous with the external field" begin + using TransitionMatrices: MieTransitionMatrix, _plane_wave_coefficients, + _field_from_coeffs, _mie_internal_coefficients, + internal_field, total_field + + # Maxwell boundary conditions on a sphere: tangential 𝐄 and 𝐇 are continuous + # across the surface. Comparing the (validated) external field to the internal + # reconstruction at r = a is a complete physical check of the internal coeffs. + tangential(E, + ϑ̂, + φ̂) = (ϑ̂[1] * E[1] + ϑ̂[2] * E[2] + ϑ̂[3] * E[3], + φ̂[1] * E[1] + φ̂[2] * E[2] + φ̂[3] * E[3]) + + λ = 2π + k0 = 2π / λ + x, mᵣ = 3.0, 1.5 + 0.05im + a = x / k0 + N = ceil(Int, x + 4cbrt(x) + 2) + mie = MieTransitionMatrix{ComplexF64, N}(x, mᵣ) + Eθ, Eφ = 1.0 + 0im, 0.0im + + a_inc, b_inc = _plane_wave_coefficients(N, 0.0, 0.0, ComplexF64(Eθ), ComplexF64(Eφ)) + p = scattering_coefficients(mie, 0.0, 0.0, Eθ, Eφ)[1] + q = scattering_coefficients(mie, 0.0, 0.0, Eθ, Eφ)[2] + cₙ, dₙ = _mie_internal_coefficients(x, mᵣ, N) + c = similar(a_inc) + d = similar(b_inc) + for n in 1:N, m in (-n):n + + c[m, n] = cₙ[n] * a_inc[m, n] + d[m, n] = dₙ[n] * b_inc[m, n] + end + k_int = mᵣ * k0 + + angles = [(ϑ, φ) for ϑ in (0.3, 0.9, 1.5, 2.2, 2.8) for φ in (0.0, 1.9, 3.5)] + discont = map(angles) do (ϑ, φ) + sϑ, cϑ, sφ, cφ = sin(ϑ), cos(ϑ), sin(φ), cos(φ) + r̂ = [sϑ * cφ, sϑ * sφ, cϑ] + ϑ̂ = [cϑ * cφ, cϑ * sφ, -sϑ] + φ̂ = [-sφ, cφ, 0.0] + pos = a .* r̂ + Eext = _field_from_coeffs(a_inc, b_inc, :regular, k0, pos, N) + + _field_from_coeffs(p, q, :outgoing, k0, pos, N) + Eint = _field_from_coeffs(c, d, :regular, k_int, pos, N) + # H̃ = k·(swap M↔N): _field_from_coeffs(c₂, c₁) gives Σ c₁ N + c₂ M ∝ ∇×E. + Hext = k0 * (_field_from_coeffs(b_inc, a_inc, :regular, k0, pos, N) + + _field_from_coeffs(q, p, :outgoing, k0, pos, N)) + Hint = k_int * _field_from_coeffs(d, c, :regular, k_int, pos, N) + Eϑe, Eφe = tangential(Eext, ϑ̂, φ̂) + Eϑi, Eφi = tangential(Eint, ϑ̂, φ̂) + Hϑe, Hφe = tangential(Hext, ϑ̂, φ̂) + Hϑi, Hφi = tangential(Hint, ϑ̂, φ̂) + relE = (abs(Eϑe - Eϑi) + abs(Eφe - Eφi)) / max(abs(Eϑe), abs(Eφe)) + relH = (abs(Hϑe - Hϑi) + abs(Hφe - Hφi)) / max(abs(Hϑe), abs(Hφe)) + (relE, relH) + end + @test maximum(first, discont) < 1e-8 + @test maximum(last, discont) < 1e-8 + + # The one-shot convenience form agrees with the manual reconstruction. + pos = a .* [sin(0.9) * cos(1.9), sin(0.9) * sin(1.9), cos(0.9)] ./ 2 # inside + E1 = internal_field(x, mᵣ, λ, 0.0, 0.0, Eθ, Eφ, pos; nmax = N) + E2 = _field_from_coeffs(c, d, :regular, k_int, pos, N) + @test maximum(abs, E1 - E2) < 1e-12 +end + +@testitem "Internal field converges to the Rayleigh (small-sphere) limit" begin + using TransitionMatrices: internal_field + + # As x → 0 the internal field becomes uniform and equal to the electrostatic + # result 3/(mᵣ²+2)·E₀ (E₀ = x̂ for +z incidence, x-polarization). This is an + # independent absolute-magnitude anchor (boundary continuity alone fixes the + # internal field only relative to the external one). The leading dynamic + # correction is O(x), so halving x roughly halves the deviation. + λ = 2π + E0 = [1.0 + 0im, 0.0im, 0.0im] + devs = Float64[] + for x in (0.02, 0.005), mᵣ in (1.5 + 0.0im, 2.0 + 0.1im) + + a = x / (2π / λ) + pred = 3 / (mᵣ^2 + 2) + pts = (0.3a .* [1.0, 0, 0], 0.3a .* [0, 1.0, 0], 0.2a .* [0, 0, 1.0]) + w = maximum(pts) do p + E = internal_field(x, mᵣ, λ, 0.0, 0.0, 1.0 + 0im, 0.0im, Vector{Float64}(p)) + maximum(abs, E .- pred .* E0) / abs(pred) + end + push!(devs, w) + @test w < 0.01 # within 1% of the static limit already at x ≤ 0.02 + end + # And it converges (smaller x ⇒ smaller deviation) for each (mᵣ). + @test devs[3] < devs[1] # mᵣ=1.5: x=0.005 closer than x=0.02 + @test devs[4] < devs[2] # mᵣ=2.0+0.1im likewise +end + +@testitem "Incident plane-wave expansion reconstructs the analytic plane wave" begin + using TransitionMatrices: _plane_wave_coefficients, _field_from_coeffs, incident_field + + # Σ aₘₙ RgMₘₙ + bₘₙ RgNₘₙ must equal the analytic plane wave inside the + # circumscribing sphere — this locks the incident coefficients and the regular + # VSWF normalization simultaneously. + λ = 2π + k = 2π / λ + N = 20 + pts = ([0.2, 0.1, 0.15], [-0.18, 0.25, -0.2], [0.05, -0.3, 0.22]) + for (Eθ, Eφ) in ((1.0 + 0im, 0.0im), (0.0im, 1.0 + 0im), (0.6 + 0.2im, 0.3 - 0.4im)), + (ϑi, φi) in ((0.7, 0.4), (2.1, -1.3)) + + a, b = _plane_wave_coefficients(N, ϑi, φi, ComplexF64(Eθ), ComplexF64(Eφ)) + for p in pts + Erec = _field_from_coeffs(a, b, :regular, k, p, N) + Epw = incident_field(λ, ϑi, φi, Eθ, Eφ, p) + @test maximum(abs, Erec - Epw) < 1e-10 + end + end +end + +@testitem "Regular field reconstruction is finite at the origin" begin + using TransitionMatrices: _plane_wave_coefficients, _field_from_coeffs, incident_field, + internal_field + + λ = 2π + k = 2π / λ + N = 20 + pos = [0.0, 0.0, 0.0] + for (Eθ, Eφ) in ((1.0 + 0im, 0.0im), (0.0im, 1.0 + 0im), + (0.6 + 0.2im, 0.3 - 0.4im)), + (ϑi, φi) in ((0.7, 0.4), (2.1, -1.3), (0.0, 0.0)) + + a, b = _plane_wave_coefficients(N, ϑi, φi, ComplexF64(Eθ), ComplexF64(Eφ)) + Erec = _field_from_coeffs(a, b, :regular, k, pos, N) + Epw = incident_field(λ, ϑi, φi, Eθ, Eφ, pos) + @test all(isfinite, abs.(Erec)) + @test maximum(abs, Erec - Epw) < 1e-10 + end + + Eint = internal_field(2.5, 1.5 + 0.05im, λ, 0.0, 0.0, 1.0 + 0im, 0.0im, pos) + @test all(isfinite, abs.(Eint)) +end + +@testitem "Scattered field reproduces amplitude_matrix in the far zone" begin + using TransitionMatrices: MieTransitionMatrix, scattering_coefficients, + scattered_field, amplitude_matrix + + dot3(a, b) = a[1] * b[1] + a[2] * b[2] + a[3] * b[3] + function bvecs(ϑ, φ) + sϑ, cϑ = sincos(ϑ) + sφ, cφ = sincos(φ) + [sϑ * cφ, sϑ * sφ, cϑ], [cϑ * cφ, cϑ * sφ, -sϑ], [-sφ, cφ, 0.0] + end + + λ = 2π + k = 2π / λ + x = 3.0 + N = ceil(Int, x + 4cbrt(x) + 2) + 𝐓 = MieTransitionMatrix{ComplexF64, N}(x, 1.5 + 0.01im) + Eθ, Eφ = 0.5 + 0.3im, -0.2 + 0.4im + ϑi, φi = 0.6, 0.5 + p, q = scattering_coefficients(𝐓, ϑi, φi, Eθ, Eφ) + + # The reconstructed far field converges to (e^{ikR}/R)·S·E₀ as O(1/R); checking + # the ratio at two radii a decade apart confirms both the value and the rate, + # which locks the outgoing-VSWF normalization and the (p,q)=𝐓(a,b) application. + for (ϑs, φs) in ((0.4, 0.3), (1.5, 2.0), (2.6, -0.7)) + n̂s, θ̂s, φ̂s = bvecs(ϑs, φs) + S = amplitude_matrix(𝐓, ϑi, φi, ϑs, φs; λ) + predθ(R) = cis(k * R) / R * (S[1, 1] * Eθ + S[1, 2] * Eφ) + predφ(R) = cis(k * R) / R * (S[2, 1] * Eθ + S[2, 2] * Eφ) + e = Float64[] + for R in (2.0e4, 2.0e5) + E = scattered_field(p, q, λ, R .* n̂s) + push!(e, abs(dot3(θ̂s, E) / predθ(R) - 1) + abs(dot3(φ̂s, E) / predφ(R) - 1)) + end + @test e[1] < 2e-3 # close already at R = 2e4 + @test e[2] < e[1] / 5 # and converging as O(1/R) + end +end diff --git a/src/common/vswf.jl b/src/common/vswf.jl new file mode 100644 index 0000000..5541137 --- /dev/null +++ b/src/common/vswf.jl @@ -0,0 +1,235 @@ +# ── Vector spherical wave functions (VSWFs) for near-field reconstruction ───── +# +# The incident, scattered and internal fields of a T-matrix problem are expanded +# in vector spherical wave functions. We use the convention of Mishchenko, Travis +# & Lacis (2002), Eqs. (5.16)–(5.17), which is the convention the rest of this +# package (`amplitude_matrix`, the T-matrix itself) already follows: +# +# 𝐌ₘₙ(k𝐫) = γₘₙ e^{imφ} zₙ(kr) [ i πₘₙ(ϑ) 𝛝̂ − τₘₙ(ϑ) 𝛗̂ ] +# 𝐍ₘₙ(k𝐫) = γₘₙ e^{imφ} { n(n+1)/(kr) zₙ(kr) Pₘₙ(ϑ) 𝐫̂ +# + [kr·zₙ(kr)]′/(kr) [ τₘₙ(ϑ) 𝛝̂ + i πₘₙ(ϑ) 𝛗̂ ] } +# +# with +# γₘₙ = (−1)^m √[(2n+1) / (4π n(n+1))], +# Pₘₙ(ϑ) = d^n_{0m}(ϑ) (Wigner d, = normalized associated Legendre), +# πₘₙ(ϑ) = m d^n_{0m}(ϑ)/sinϑ, +# τₘₙ(ϑ) = d d^n_{0m}(ϑ)/dϑ, +# and zₙ the spherical Bessel function jₙ (`kind = :regular`, used for the regular +# functions RgM/RgN that expand the incident and internal fields) or the spherical +# Hankel function hₙ⁽¹⁾ (`kind = :outgoing`, the radiating M/N that expand the +# scattered field, valid outside the circumscribing sphere). +# +# The radial factors are obtained from the package's Riccati–Bessel routines +# (`ricattibesselj`, `ricattibessely`): +# jₙ(x) = ψₙ(x)/x, yₙ(x) = χₙ(x)/x, hₙ⁽¹⁾(x) = jₙ(x) + i yₙ(x), +# [x·zₙ(x)]′/x = ψ′ₙ(x)/x (regular) resp. (ψ′ₙ + i χ′ₙ)/x (outgoing). +# +# 𝐌 and 𝐍 share the same γₘₙ, so they satisfy ∇×𝐌ₘₙ = k 𝐍ₘₙ and ∇×𝐍ₘₙ = k 𝐌ₘₙ +# regardless of the absolute value of γₘₙ; this curl identity is what the unit +# tests check. The absolute normalization γₘₙ is fixed to the Mishchenko value so +# that, combined with the plane-wave incident coefficients, the reconstructed +# field reproduces `amplitude_matrix` in the far zone (validated against Mie). + +@inline function _sph_unit_vectors(ϑ::T, φ::T) where {T} + sϑ, cϑ = sincos(ϑ) + sφ, cφ = sincos(φ) + r̂ = SVector(sϑ * cφ, sϑ * sφ, cϑ) + ϑ̂ = SVector(cϑ * cφ, cϑ * sφ, -sϑ) + φ̂ = SVector(-sφ, cφ, zero(T)) + return r̂, ϑ̂, φ̂ +end + +@inline _is_regular_origin(kind::Symbol, x, + ::Type{T}) where {T} = kind === :regular && (iszero(x) || abs(x) < eps(T)) + +function _regular_vswf_origin_limit(::Type{T}, m::Integer, n::Integer) where {T} + CT = Complex{T} + zero3 = zero(SVector{3, CT}) + n == 1 || return zero3, zero3 + + α = inv(sqrt(T(12) * T(π))) + β = inv(sqrt(T(6) * T(π))) + 𝐍 = if m == -1 + SVector{3, CT}(α, -im * α, zero(CT)) + elseif m == 0 + SVector{3, CT}(zero(CT), zero(CT), β) + elseif m == 1 + SVector{3, CT}(-α, -im * α, zero(CT)) + else + zero3 + end + return zero3, 𝐍 +end + +function _regular_vswf_radial_origin(::Type{T}, nₘₐₓ::Integer, x::Number) where {T} + XT = x isa Real ? T : Complex{T} + zₙ = zeros(XT, nₘₐₓ) + xzₙ′ = zeros(XT, nₘₐₓ) + nₘₐₓ ≥ 1 && (xzₙ′[1] = XT(T(2) / T(3))) + return zₙ, xzₙ′ +end + +# Radial factors zₙ(x) and [x·zₙ(x)]′/x for n = 1:nₘₐₓ, plus zₙ(x)/x (the factor +# multiplying the 𝐫̂ component of 𝐍). +function _vswf_radial(kind::Symbol, nₘₐₓ::Integer, x::Number) + T = float(real(typeof(x))) + _is_regular_origin(kind, x, T) && return _regular_vswf_radial_origin(T, nₘₐₓ, x) + + if kind === :regular + ψ, ψ′ = ricattibesselj(nₘₐₓ, estimate_ricattibesselj_extra_terms(nₘₐₓ, x), x) + zₙ = ψ ./ x + xzₙ′ = ψ′ ./ x + elseif kind === :outgoing + ψ, ψ′ = ricattibesselj(nₘₐₓ, estimate_ricattibesselj_extra_terms(nₘₐₓ, x), x) + χ, χ′ = ricattibessely(nₘₐₓ, x) + zₙ = (ψ .+ im .* χ) ./ x + xzₙ′ = (ψ′ .+ im .* χ′) ./ x + else + throw(ArgumentError("`kind` must be :regular or :outgoing, got $kind")) + end + return zₙ, xzₙ′ +end + +@doc raw""" +``` +vswf(kind::Symbol, m::Integer, n::Integer, k, r, ϑ, φ) +``` + +Evaluate the vector spherical wave functions ``\mathbf{M}_{mn}`` and +``\mathbf{N}_{mn}`` of degree `m` and order `n` at the point `(r, ϑ, φ)` for size +parameter `k` (so the radial argument is `x = k·r`), returning the pair `(𝐌, 𝐍)` +as `SVector{3}` of **Cartesian** field components. + +`kind` selects the radial dependence: + +- `:regular` — spherical Bessel ``j_n`` (the regular `RgM`/`RgN`, finite at the + origin; used for incident and internal fields). +- `:outgoing` — spherical Hankel ``h_n^{(1)}`` (the radiating `M`/`N` used for the + scattered field; valid outside the circumscribing sphere). + +The convention is Mishchenko, Travis & Lacis (2002), Eqs. (5.16)–(5.17), +consistent with [`amplitude_matrix`](@ref). `𝐌` and `𝐍` satisfy +``\nabla\times\mathbf{M}_{mn} = k\,\mathbf{N}_{mn}`` and +``\nabla\times\mathbf{N}_{mn} = k\,\mathbf{M}_{mn}``. +""" +function vswf(kind::Symbol, m::Integer, n::Integer, k::Real, r::Real, ϑ::Real, φ::Real) + abs(m) ≤ n || throw(ArgumentError("require |m| ≤ n, got m=$m, n=$n")) + T = float(promote_type(typeof(k), typeof(r), typeof(ϑ), typeof(φ))) + ϑ, φ = T(ϑ), T(φ) + x = T(k) * T(r) + _is_regular_origin(kind, x, T) && return _regular_vswf_origin_limit(T, m, n) + + zₙ, xzₙ′ = _vswf_radial(kind, n, x) + zn = zₙ[n] + dzn = xzₙ′[n] + zn_over_x = zn / x + + # Angular factors Pₘₙ = d^n_{0m}, τₘₙ = d(d^n_{0m})/dϑ, πₘₙ = m d^n_{0m}/sinϑ. + dvec, τvec = wigner_d_recursion(T, 0, m, n, ϑ; deriv = true) + P = dvec[n] + τ = τvec[n] + π_ = pi_func(T, m, n, ϑ; d = P) + + γ = (-1)^(m & 1) * sqrt(T(2n + 1) / (4 * T(π) * n * (n + 1))) + eⁱᵐᵠ = cis(m * φ) + r̂, ϑ̂, φ̂ = _sph_unit_vectors(ϑ, φ) + + pref = γ * eⁱᵐᵠ + 𝐌 = pref * zn * (im * π_ .* ϑ̂ .- τ .* φ̂) + 𝐍 = pref * + (n * (n + 1) * zn_over_x * P .* r̂ .+ dzn .* (τ .* ϑ̂ .+ im * π_ .* φ̂)) + + return 𝐌, 𝐍 +end + +""" +``` +vswf_cartesian(kind::Symbol, m, n, k, xyz::AbstractVector) +``` + +Convenience wrapper of [`vswf`](@ref) taking a Cartesian position `xyz = (x, y, z)` +and returning the `(𝐌, 𝐍)` Cartesian field vectors. The polar axis is `+z`. +""" +function vswf_cartesian(kind::Symbol, m::Integer, n::Integer, k::Real, + xyz::AbstractVector) + x, y, z = xyz[1], xyz[2], xyz[3] + T = float(promote_type(typeof(k), typeof(x), typeof(y), typeof(z))) + xT, yT, zT = T(x), T(y), T(z) + r = sqrt(xT^2 + yT^2 + zT^2) + _is_regular_origin(kind, T(k) * r, T) && return _regular_vswf_origin_limit(T, m, n) + ϑ = acos(clamp(zT / r, -one(T), one(T))) + φ = atan(yT, xT) + return vswf(kind, m, n, k, r, ϑ, φ) +end + +@testitem "VSWF satisfy the curl identities ∇×M = kN and ∇×N = kM" begin + using TransitionMatrices: vswf_cartesian + + # Finite-difference Cartesian curl of a vector field F(𝐫). + function fd_curl(F, p, h) + ex = [h, 0.0, 0.0] + ey = [0.0, h, 0.0] + ez = [0.0, 0.0, h] + ∂x = (F(p + ex) - F(p - ex)) / 2h + ∂y = (F(p + ey) - F(p - ey)) / 2h + ∂z = (F(p + ez) - F(p - ez)) / 2h + return [∂y[3] - ∂z[2], ∂z[1] - ∂x[3], ∂x[2] - ∂y[1]] + end + + k = 1.3 + h = 1e-5 + # A handful of generic off-axis points (avoid the poles / origin). + points = ([0.7, 0.5, 0.9], [-0.6, 0.8, 0.4], [1.1, -0.3, -0.7]) + + @testset "kind = $kind" for kind in (:regular, :outgoing) + for p in points + for n in 1:4, m in (-n):n + + M(q) = vswf_cartesian(kind, m, n, k, q)[1] + N(q) = vswf_cartesian(kind, m, n, k, q)[2] + curlM = fd_curl(M, p, h) + curlN = fd_curl(N, p, h) + kN = k * vswf_cartesian(kind, m, n, k, p)[2] + kM = k * vswf_cartesian(kind, m, n, k, p)[1] + @test maximum(abs, curlM - kN) / maximum(abs, kN) < 1e-6 + @test maximum(abs, curlN - kM) / maximum(abs, kM) < 1e-6 + end + end + end +end + +@testitem "Regular VSWFs have finite analytic origin limits" begin + using TransitionMatrices: _vswf_radial, vswf, vswf_cartesian + + α = 1 / sqrt(12π) + β = 1 / sqrt(6π) + zero3 = [0.0, 0.0, 0.0] + + zₙ, xzₙ′ = _vswf_radial(:regular, 4, 0.0) + @test all(isfinite, zₙ) + @test all(isfinite, xzₙ′) + @test maximum(abs, zₙ) == 0 + @test xzₙ′[1] == 2 / 3 + @test maximum(abs, xzₙ′[2:end]) == 0 + + expected_N(m, + n) = n == 1 ? + ( + m == -1 ? α .* [1.0, -im, 0.0] : + m == 0 ? β .* [0.0, 0.0, 1.0] : + m == 1 ? -α .* [1.0, im, 0.0] : + zeros(ComplexF64, 3)) : zeros(ComplexF64, 3) + + for n in 1:4, m in (-n):n + + M, N = vswf_cartesian(:regular, m, n, 1.0, zero3) + @test all(isfinite, abs.(M)) + @test all(isfinite, abs.(N)) + @test maximum(abs, M) < 1e-14 + @test maximum(abs, N - expected_N(m, n)) < 1e-14 + + M_sph, N_sph = vswf(:regular, m, n, 1.0, 0.0, 0.4, 0.2) + @test maximum(abs, M_sph - M) < 1e-14 + @test maximum(abs, N_sph - N) < 1e-14 + end +end