vSmartMOM

Appearance

Mie Scattering Deep Dive ​

This tutorial explains Mie theory in vSmartMOM, the NAI2 vs PCW methods, parameter selection, numerical precision, phase reconstruction, delta-M truncation, and how aerosol optics feed into RT.

1. Introduction: What Mie Theory Computes ​

Mie theory solves Maxwell's equations for electromagnetic scattering by a homogeneous sphere. For radiative transfer (RT), we need:

  • Extinction cross-section : total attenuation (scattering + absorption)

  • Scattering cross-section : energy scattered

  • Phase function : angular distribution of scattered light

The phase matrix (for polarized RT) is expanded in generalized spherical functions; the expansion coefficients are the "Greek" coefficients . RT kernels work directly in this Fourier/Legendre space—angle-space values are reconstructed only when needed (e.g., for plotting or diagnostics).

julia
using vSmartMOM.Scattering
using Distributions

2. NAI2 vs PCW Methods ​

NAI2 (Numerical Angular Integration) — Siewert's method: directly integrates the phase matrix over scattering angles to obtain Greek coefficients. No precomputation; suitable when you compute a few aerosols or vary parameters often.

PCW (Precomputed Wigner) — Domke's method: precomputes Wigner d-functions for all indices. The Greek coefficients are then evaluated via fast summations over Mie indices. Use PCW when you have many quadrature angles (e.g., high-resolution RT with many streams)—the one-time Wigner cost is amortized.

Both methods produce identical Greek coefficients; choose based on workflow:

ScenarioPreferred method
Single aerosol, few runsNAI2
Many angles, repeated runsPCW
Parameter sweepsNAI2 (no Wigner)

Example: NAI2 (no Wigner tables needed)

julia
rₘ, σ, nᵣ, nᵢ = 0.3, 2.0, 1.3, 0.001
aero = Aerosol(LogNormal(log(rₘ), log(σ)), nᵣ, nᵢ)
λ = 0.55
model_NAI2 = make_mie_model(NAI2(), aero, λ, Stokes_IQUV(), δBGE(20), 30.0, 500)

aerosol_optics = compute_aerosol_optical_properties(model_NAI2)

PCW requires precomputed Wigner tables (expensive one-time cost): N_max = 120 wigner_A, wigner_B = compute_wigner_values(N_max) model_PCW = make_mie_model(PCW(), aero, λ, Stokes_IQUV(), δBGE(20), 30.0, 500, wigner_A, wigner_B)

3. Parameter Selection Guide ​

r_max — Upper bound for the size-distribution integral. Must cover the bulk of the distribution; truncating too low biases cross-sections. For a lognormal distribution, select a high distribution quantile and check convergence of extinction, SSA, and phase coefficients as the bound increases; a small omitted number fraction can still carry substantial cross section.

nquad_radius — Gauss–Legendre points in over the fixed interval (r_min keyword of make_mie_model, default 1e-3 µm). The nodes do not depend on the size distribution, so they must resolve the narrowest mode and the Mie oscillations of the largest particles. Typical: 500–2000.

l_max / l_trunc — Maximum Legendre order in the Greek expansion. Must exceed the effective size parameter;  . Too low truncates the phase function; too high adds cost with little gain.

Trade-off: Higher nquad_radius and l_max improve accuracy but increase runtime. Start with moderate values (e.g., 500, 20) and increase if needed.

julia
r_max = 30.0
nquad_radius = 500
l_max = 20
truncation_type = δBGE(l_max)

4. Float32 vs Float64 Behavior ​

Mie computations use downward recurrence for the logarithmic derivative and upward recurrence for Bessel functions. Float32 can suffer from:

  • Loss of precision in the recurrence for large size parameters

  • Instability in , for  

Recommendation: Use Float64 for Mie. Float32 is acceptable only for small size parameters (e.g.,   ) and when you have verified stability. The RT model may use Float32 elsewhere; convert aerosol optics to Float64 for Mie, then cast back if needed.

5. Phase Function Reconstruction ​

Greek coefficients encode the phase matrix in Legendre space. To get angle-dependent values:

Use reconstruct_phase(greek_coefs, μ) where μ = cos(Θ). The scalar phase function is f₁₁; it is normalized so   .

Visualization: Plot f₁₁ vs μ (or scattering angle Θ); use log scale for f₁₁ to see the forward peak. The ratio f₁₂/f₁₁ gives the degree of linear polarization.

Minimal reconstruction example (no heavy compute): μ, _ = gausslegendre(100) smat = reconstruct_phase(aerosol_optics.greek_coefs, μ)

plot(μ, smat.f₁₁, yscale=:log10); plot(μ, smat.f₁₂ ./ smat.f₁₁) ​

6. Delta-M Truncation ​

Strong forward peaks (large particles) require many Legendre terms. Delta-M (δ-BGE) approximates the forward peak and renormalizes, reducing the effective expansion length.

Truncation factor : fraction of scattered energy moved to the "delta" term:   , where is the retained scattering fraction. The RT solver uses to adjust the single-scattering albedo and phase function so that the truncated expansion approximates the original phase function on the fitted domain.

is retired and forced to zero. Do not tune an exclusion cone; choose the stream/truncation order and validate the resulting phase/radiance.

7. Connection to RT ​

AerosolOptics provides exactly what the RT model needs:

  • greek_coefs — Fourier/Legendre coefficients for the phase matrix

  • ω̃ — single-scattering albedo ()

  • k — extinction cross-section (per particle or per mass, depending on convention)

  • fᵗ — delta-M truncation factor (set by truncate_phase; the untruncated output of compute_aerosol_optical_properties carries the placeholder 1)

These feed into the doubling-adding kernels. Aerosol layers are combined with molecular (Rayleigh) scattering and absorption; the RT solver never needs angle-space phase functions—only the Greek expansion.

julia
println("Mie deep-dive tutorial complete. Key types: Aerosol, MieModel, AerosolOptics, GreekCoefs.")

This page was generated using Literate.jl.