from __future__ import annotations
import json
from pathlib import Path
def md(*lines: str) -> dict:
return {"cell_type": "markdown", "metadata": {}, "source": list(lines)}
def code(*lines: str) -> dict:
return {
"cell_type": "code",
"metadata": {},
"execution_count": None,
"outputs": [],
"source": list(lines),
}
def notebook(cells: list[dict]) -> dict:
return {
"cells": cells,
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3",
},
"language_info": {
"name": "python",
"version": "3.11",
},
},
"nbformat": 4,
"nbformat_minor": 5,
}
NB01 = [
md(
"# Tutorial 01 — A dielectric sphere at X-band, checked against Mie\n",
"\n",
"A T-matrix solver for nonspherical particles must reduce to classical\n",
"Mie theory in the `axis_ratio = 1` limit. That reduction is the\n",
"simplest end-to-end check we can do. This tutorial builds a 1 mm\n",
"water sphere at X-band, computes its scattering and extinction\n",
"cross-sections via the T-matrix path, and compares against the\n",
"closed-form Mie expressions shipped with rustmatrix.\n",
),
code(
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"\n",
"from rustmatrix import Scatterer, mie_qsca, mie_qext, scatter\n",
"from rustmatrix.tmatrix_aux import geom_horiz_back, wl_X\n",
"from rustmatrix.refractive import m_w_10C\n",
),
md("## Build a sphere scatterer at X-band\n"),
code(
"radius_mm = 1.0\n",
"wavelength_mm = wl_X\n",
"m = m_w_10C[wl_X]\n",
"\n",
"s = Scatterer(radius=radius_mm, wavelength=wavelength_mm, m=m,\n",
" axis_ratio=1.0, ddelt=1e-4, ndgs=2)\n",
"s.set_geometry(geom_horiz_back)\n",
"S, Z = s.get_SZ()\n",
"S\n",
),
md("## Cross-sections — T-matrix vs Mie\n"),
code(
"size_param = 2.0 * np.pi * radius_mm / wavelength_mm\n",
"sigma_sca_tm = scatter.sca_xsect(s, h_pol=True)\n",
"sigma_ext_tm = scatter.ext_xsect(s, h_pol=True)\n",
"\n",
"geom = np.pi * radius_mm ** 2\n",
"sigma_sca_mie = mie_qsca(size_param, m.real, m.imag) * geom\n",
"sigma_ext_mie = mie_qext(size_param, m.real, m.imag) * geom\n",
"\n",
"print(f'rel err σ_sca = {abs(sigma_sca_tm-sigma_sca_mie)/sigma_sca_mie:.2e}')\n",
"print(f'rel err σ_ext = {abs(sigma_ext_tm-sigma_ext_mie)/sigma_ext_mie:.2e}')\n",
),
md("## Q_sca and Q_ext across a size-parameter sweep\n",
"\n",
"Plotting both curves side by side makes the Mie ripple pattern\n",
"visible and confirms that the T-matrix follows it point-for-point.\n"),
code(
"xs = np.linspace(0.1, 10.0, 40)\n",
"q_sca_mie = np.array([mie_qsca(x, m.real, m.imag) for x in xs])\n",
"q_ext_mie = np.array([mie_qext(x, m.real, m.imag) for x in xs])\n",
"q_sca_tm = np.empty_like(xs)\n",
"q_ext_tm = np.empty_like(xs)\n",
"for i, x in enumerate(xs):\n",
" r = x * wavelength_mm / (2.0 * np.pi)\n",
" si = Scatterer(radius=r, wavelength=wavelength_mm, m=m,\n",
" axis_ratio=1.0, ddelt=1e-4, ndgs=2)\n",
" si.set_geometry(geom_horiz_back)\n",
" g = np.pi * r ** 2\n",
" q_sca_tm[i] = scatter.sca_xsect(si) / g\n",
" q_ext_tm[i] = scatter.ext_xsect(si) / g\n",
"\n",
"fig, ax = plt.subplots(figsize=(7, 4))\n",
"ax.plot(xs, q_sca_mie, 'k-', label='Q_sca (Mie)')\n",
"ax.plot(xs, q_ext_mie, 'k--', label='Q_ext (Mie)')\n",
"ax.plot(xs, q_sca_tm, 'C1o', markersize=4, label='Q_sca (T-matrix)')\n",
"ax.plot(xs, q_ext_tm, 'C2s', markersize=4, label='Q_ext (T-matrix)')\n",
"ax.set_xlabel('size parameter x = 2πr/λ')\n",
"ax.set_ylabel('efficiency')\n",
"ax.legend()\n",
"ax.grid(True, alpha=0.3)\n",
"fig.tight_layout();\n",
),
]
NB02 = [
md(
"# Tutorial 02 — Single-drop polarimetric response at S/C/X bands\n",
"\n",
"Falling raindrops are flattened by drag; larger drops are more\n",
"oblate (Thurai et al. 2007). That shape anisotropy imprints four\n",
"distinct signatures on dual-polarisation radar data:\n",
"\n",
"* **Z_h** grows as D⁶ in the Rayleigh regime and then walks out of\n",
" it — earliest at X-band, latest at S-band.\n",
"* **Z_dr** rises with D because oblateness grows with D; the\n",
" wavelength dependence exposes C-band's resonance bump near\n",
" D ≈ 5 mm.\n",
"* **K_dp** scales with Re(f_h(0) − f_v(0)) — strictly stronger at\n",
" shorter wavelengths; X-band K_dp per drop is ≈ 2× C-band and\n",
" ≈ 4× S-band.\n",
"* **LDR** — linear depolarisation ratio — is set by the canting\n",
" distribution. Here we model a σ = 5° Gaussian wobble around the\n",
" flat-lying orientation to produce realistic rain LDR in the\n",
" −30 to −25 dB range for 5+ mm drops.\n",
"\n",
"This notebook sweeps drop equivalent diameter D = 0.1–8 mm at S,\n",
"C, and X bands and plots all four observables. We report Z_h and\n",
"K_dp *per drop/m³* so multiplying by the drop concentration\n",
"N [m⁻³] gives the usual bulk observables.\n",
),
code(
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"\n",
"from rustmatrix import Scatterer, orientation, radar, scatter\n",
"from rustmatrix.tmatrix_aux import (K_w_sqr, dsr_thurai_2007,\n",
" geom_horiz_back, geom_horiz_forw,\n",
" wl_C, wl_S, wl_X)\n",
"from rustmatrix.refractive import m_w_10C\n",
"\n",
"BANDS = [('S', wl_S, 'C0'), ('C', wl_C, 'C1'), ('X', wl_X, 'C2')]\n",
"D_GRID = np.linspace(0.1, 8.0, 40)\n",
"CANTING_STD_DEG = 5.0\n",
),
md("## Build the drop and run the sweep\n",
"\n",
"One `Scatterer` per (D, λ) point, with a σ = 5° Gaussian canting\n",
"PDF around β = 0° (flat-lying oblate drop with modest turbulent\n",
"wobble). Backscatter geometry gives Z_h, Z_dr, and LDR; forward\n",
"geometry gives K_dp.\n"),
code(
"def build_drop(D_mm, wl):\n",
" s = Scatterer(radius=D_mm/2, wavelength=wl, m=m_w_10C[wl],\n",
" axis_ratio=1.0/dsr_thurai_2007(D_mm),\n",
" Kw_sqr=K_w_sqr[wl], ddelt=1e-4, ndgs=2)\n",
" s.orient = orientation.orient_averaged_fixed\n",
" s.or_pdf = orientation.gaussian_pdf(std=CANTING_STD_DEG, mean=0.0)\n",
" s.n_alpha = 6; s.n_beta = 12\n",
" return s\n",
"\n",
"def sweep(wl):\n",
" Zh = np.empty_like(D_GRID); Zdr = np.empty_like(D_GRID)\n",
" Kdp = np.empty_like(D_GRID); LDR = np.empty_like(D_GRID)\n",
" for i, D in enumerate(D_GRID):\n",
" s = build_drop(D, wl)\n",
" s.set_geometry(geom_horiz_back)\n",
" Zh[i] = 10*np.log10(max(radar.refl(s, h_pol=True), 1e-30))\n",
" Zdr[i] = 10*np.log10(max(radar.Zdr(s), 1e-30))\n",
" LDR[i] = 10*np.log10(max(scatter.ldr(s, h_pol=True), 1e-30))\n",
" s.set_geometry(geom_horiz_forw)\n",
" Kdp[i] = radar.Kdp(s)\n",
" return dict(Zh=Zh, Zdr=Zdr, Kdp=Kdp, LDR=LDR)\n",
"\n",
"data = {name: sweep(wl) for name, wl, _ in BANDS}\n",
),
md("## Plot all four observables vs. D\n",
"\n",
"Top-left Z_h tracks the D⁶ line until it bends: X-band breaks\n",
"earliest (shortest λ, smallest χ = πD/λ needed), S-band last.\n",
"Top-right Z_dr rises monotonically; the C-band curve bumps above\n",
"X and S around D ≈ 5 mm — the well-known C-band raindrop\n",
"resonance. Bottom-left K_dp ordering is X > C > S at every D.\n",
"Bottom-right LDR is a clean fingerprint of the canting distribution\n",
"and rises smoothly with D once the drop is oblate enough to leak\n",
"cross-pol power.\n"),
code(
"fig, axes = plt.subplots(2, 2, figsize=(11, 7), sharex=True)\n",
"\n",
"ref_D = D_GRID[(D_GRID > 0.5) & (D_GRID < 2.5)]\n",
"rayleigh = 10*np.log10(ref_D**6) + (data['S']['Zh'][10] - 10*np.log10(D_GRID[10]**6))\n",
"axes[0, 0].plot(ref_D, rayleigh, 'k:', lw=1, label='D⁶ (Rayleigh)')\n",
"for name, _, c in BANDS:\n",
" axes[0, 0].plot(D_GRID, data[name]['Zh'], color=c, lw=1.8, label=f'{name}-band')\n",
"axes[0, 0].set_ylabel('Z_h [dBZ per drop/m³]')\n",
"axes[0, 0].legend(fontsize=9)\n",
"\n",
"for name, _, c in BANDS:\n",
" axes[0, 1].plot(D_GRID, data[name]['Zdr'], color=c, lw=1.8, label=f'{name}-band')\n",
"axes[0, 1].set_ylabel('Z_dr [dB]')\n",
"axes[0, 1].legend(fontsize=9)\n",
"\n",
"for name, _, c in BANDS:\n",
" axes[1, 0].semilogy(D_GRID, np.abs(data[name]['Kdp']),\n",
" color=c, lw=1.8, label=f'{name}-band')\n",
"axes[1, 0].set_ylabel('|K_dp| [°/km per drop/m³]')\n",
"axes[1, 0].legend(fontsize=9)\n",
"\n",
"for name, _, c in BANDS:\n",
" axes[1, 1].plot(D_GRID, data[name]['LDR'], color=c, lw=1.8, label=f'{name}-band')\n",
"axes[1, 1].set_ylim(-60, -20)\n",
"axes[1, 1].set_ylabel('LDR [dB]')\n",
"axes[1, 1].legend(fontsize=9)\n",
"\n",
"for ax in axes.flat:\n",
" ax.set_xlim(0, 8)\n",
" ax.grid(True, alpha=0.3)\n",
"axes[1, 0].set_xlabel('equivalent diameter D [mm]')\n",
"axes[1, 1].set_xlabel('equivalent diameter D [mm]')\n",
"fig.suptitle(f'Single-drop response at S/C/X bands '\n",
" f'(Thurai 2007 shape, 10 °C water, σ_canting = {CANTING_STD_DEG:.0f}°)')\n",
"fig.tight_layout();\n",
),
md("## Spot values at canonical diameters\n",
"\n",
"Six diameters span the regime map: 0.5 mm (nearly spherical,\n",
"pure Rayleigh), 1–3 mm (moderate oblateness, still Rayleigh at\n",
"S/C), 5 mm (C-band resonance territory), and 7 mm (well into\n",
"non-Rayleigh at all three bands).\n"),
code(
"rows = (0.5, 1.0, 2.0, 3.0, 5.0, 7.0)\n",
"idx = [int(np.argmin(np.abs(D_GRID - D))) for D in rows]\n",
"for obs, fmt in (('Zh', '{:+7.2f}'), ('Zdr', '{:+7.3f}'),\n",
" ('Kdp', '{:+7.2e}'), ('LDR', '{:+7.2f}')):\n",
" label = {'Zh': 'Z_h [dBZ]', 'Zdr': 'Z_dr [dB]',\n",
" 'Kdp': 'K_dp [°/km]', 'LDR': 'LDR [dB]'}[obs]\n",
" header = ' '.join(f'D={D:>3.1f}' for D in rows)\n",
" print(f'{label:<14} {header}')\n",
" for name, _, _ in BANDS:\n",
" row = data[name][obs][idx]\n",
" cells = ' '.join(fmt.format(v) for v in row)\n",
" print(f' {name}-band {cells}')\n",
" print()\n",
),
]
NB03 = [
md(
"# Tutorial 03 — Gamma-PSD rain at C-band\n",
"\n",
"Radar volumes sample thousands of drops. The observed Z_h, Z_dr,\n",
"K_dp, and specific attenuation A_i are PSD-weighted integrals of\n",
"the single-drop quantities. This notebook tabulates S(D) and Z(D)\n",
"once, then sweeps the integrated observables across a range of\n",
"normalised gamma PSDs parameterised by the median volume diameter\n",
"D0.\n",
),
code(
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"\n",
"from rustmatrix import Scatterer, radar, psd as rs_psd\n",
"from rustmatrix.tmatrix_aux import (dsr_thurai_2007, geom_horiz_back,\n",
" geom_horiz_forw, K_w_sqr, wl_C)\n",
"from rustmatrix.refractive import m_w_10C\n",
),
md("## Tabulate once, evaluate many\n"),
code(
"s = Scatterer(wavelength=wl_C, m=m_w_10C[wl_C],\n",
" Kw_sqr=K_w_sqr[wl_C], ddelt=1e-4, ndgs=2)\n",
"integ = rs_psd.PSDIntegrator()\n",
"integ.D_max = 8.0\n",
"integ.num_points = 64\n",
"integ.axis_ratio_func = lambda D: 1.0 / dsr_thurai_2007(D)\n",
"integ.geometries = (geom_horiz_back, geom_horiz_forw)\n",
"s.psd_integrator = integ\n",
"s.psd_integrator.init_scatter_table(s)\n",
),
md("## Sweep median diameter D0\n"),
code(
"D0s = np.linspace(0.5, 3.0, 12)\n",
"Zh = np.empty_like(D0s)\n",
"Zdr = np.empty_like(D0s)\n",
"Kdp = np.empty_like(D0s)\n",
"Ai = np.empty_like(D0s)\n",
"\n",
"for i, D0 in enumerate(D0s):\n",
" s.psd = rs_psd.GammaPSD(D0=D0, Nw=8e3, mu=4)\n",
" s.set_geometry(geom_horiz_back)\n",
" Zh[i] = 10 * np.log10(radar.refl(s))\n",
" Zdr[i] = 10 * np.log10(radar.Zdr(s))\n",
" s.set_geometry(geom_horiz_forw)\n",
" Kdp[i] = radar.Kdp(s)\n",
" Ai[i] = radar.Ai(s)\n",
"\n",
"fig, axes = plt.subplots(2, 2, figsize=(9, 6), sharex=True)\n",
"axes[0, 0].plot(D0s, Zh, 'C0-o'); axes[0, 0].set_ylabel('Z_h [dBZ]')\n",
"axes[0, 1].plot(D0s, Zdr, 'C1-o'); axes[0, 1].set_ylabel('Z_dr [dB]')\n",
"axes[1, 0].plot(D0s, Kdp, 'C2-o'); axes[1, 0].set_ylabel('K_dp [°/km]')\n",
"axes[1, 1].plot(D0s, Ai, 'C3-o'); axes[1, 1].set_ylabel('A_i [dB/km]')\n",
"for ax in axes.flat:\n",
" ax.set_xlabel('D0 [mm]')\n",
" ax.grid(True, alpha=0.3)\n",
"fig.suptitle('C-band gamma-PSD rain observables (Nw=8e3, mu=4)')\n",
"fig.tight_layout();\n",
),
]
NB04 = [
md(
"# Tutorial 04 — Oriented ice crystals at W-band\n",
"\n",
"Ice crystals fall with a preferred orientation and a spread of\n",
"canting angles around it. The orientation PDF captures the spread\n",
"and feeds directly into the dual-pol observables. This notebook\n",
"compares the three orientation-averaging strategies rustmatrix\n",
"ships with: none, fixed-quadrature, and adaptive.\n",
),
code(
"import time\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"\n",
"from rustmatrix import Scatterer, orientation as rs_orient, radar\n",
"from rustmatrix.tmatrix_aux import geom_horiz_back, wl_W\n",
"from rustmatrix.refractive import mi\n",
),
md("## One crystal, three averaging schemes\n"),
code(
"D_eq_mm = 0.5\n",
"ice_m = mi(wl_W, 0.9)\n",
"pdf = rs_orient.gaussian_pdf(std=20.0, mean=90.0)\n",
"base = dict(radius=D_eq_mm/2, wavelength=wl_W, m=ice_m,\n",
" axis_ratio=0.5, ddelt=1e-4, ndgs=2)\n",
"\n",
"schemes = [\n",
" ('single', rs_orient.orient_single, {'or_pdf': pdf}),\n",
" ('fixed 6×12', rs_orient.orient_averaged_fixed,\n",
" {'or_pdf': pdf, 'n_alpha': 6, 'n_beta': 12}),\n",
" ('adaptive', rs_orient.orient_averaged_adaptive, {'or_pdf': pdf}),\n",
"]\n",
"rows = []\n",
"for name, orient, extra in schemes:\n",
" s = Scatterer(**base, **extra)\n",
" s.orient = orient\n",
" s.set_geometry(geom_horiz_back)\n",
" t0 = time.perf_counter()\n",
" zdr = radar.Zdr(s)\n",
" elapsed = time.perf_counter() - t0\n",
" rows.append((name, 10*np.log10(zdr), elapsed))\n",
" print(f'{name:<12} Z_dr = {10*np.log10(zdr):+.4f} dB ({elapsed*1000:.1f} ms)')\n",
),
md("## Canting-angle spread sensitivity\n",
"\n",
"Sweep the Gaussian PDF's `std` parameter to see how rapidly Z_dr\n",
"collapses as the orientation spread widens.\n"),
code(
"stds = np.array([1, 5, 10, 15, 20, 30, 45, 60])\n",
"zdrs = np.empty_like(stds, dtype=float)\n",
"for i, sd in enumerate(stds):\n",
" s = Scatterer(**base)\n",
" s.orient = rs_orient.orient_averaged_fixed\n",
" s.or_pdf = rs_orient.gaussian_pdf(std=float(sd), mean=90.0)\n",
" s.n_alpha, s.n_beta = 6, 12\n",
" s.set_geometry(geom_horiz_back)\n",
" zdrs[i] = 10 * np.log10(radar.Zdr(s))\n",
"\n",
"fig, ax = plt.subplots(figsize=(7, 4))\n",
"ax.plot(stds, zdrs, 'C0-o')\n",
"ax.set_xlabel('canting-angle σ [deg]')\n",
"ax.set_ylabel('Z_dr [dB]')\n",
"ax.set_title('W-band prolate ice column, oriented')\n",
"ax.grid(True, alpha=0.3)\n",
"fig.tight_layout();\n",
),
]
NB05 = [
md(
"# Tutorial 05 — One PSD, six radar bands\n",
"\n",
"Retrieval algorithms that use more than one frequency need a\n",
"forward model that works across the whole instrument suite. This\n",
"notebook runs the same moderate-convective gamma PSD through\n",
"rustmatrix at S, C, X, Ku, Ka, and W band and plots the frequency\n",
"response of the dual-pol observables.\n",
),
code(
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"\n",
"from rustmatrix import Scatterer, radar, psd as rs_psd\n",
"from rustmatrix.tmatrix_aux import (dsr_thurai_2007, geom_horiz_back,\n",
" geom_horiz_forw, K_w_sqr,\n",
" wl_S, wl_C, wl_X, wl_Ku,\n",
" wl_Ka, wl_W)\n",
"from rustmatrix.refractive import m_w_20C\n",
"\n",
"BANDS = [('S', wl_S), ('C', wl_C), ('X', wl_X),\n",
" ('Ku', wl_Ku), ('Ka', wl_Ka), ('W', wl_W)]\n",
),
md("## Run the six bands\n"),
code(
"def run(wl):\n",
" s = Scatterer(wavelength=wl, m=m_w_20C[wl],\n",
" Kw_sqr=K_w_sqr[wl], ddelt=1e-4, ndgs=2)\n",
" integ = rs_psd.PSDIntegrator()\n",
" integ.D_max = 8.0\n",
" integ.num_points = 64\n",
" integ.axis_ratio_func = lambda D: 1.0 / dsr_thurai_2007(D)\n",
" integ.geometries = (geom_horiz_back, geom_horiz_forw)\n",
" s.psd_integrator = integ\n",
" s.psd_integrator.init_scatter_table(s)\n",
" s.psd = rs_psd.GammaPSD(D0=1.5, Nw=8e3, mu=4)\n",
" s.set_geometry(geom_horiz_back)\n",
" Zh = 10 * np.log10(radar.refl(s))\n",
" Zdr = 10 * np.log10(radar.Zdr(s))\n",
" s.set_geometry(geom_horiz_forw)\n",
" return Zh, Zdr, radar.Kdp(s), radar.Ai(s)\n",
"\n",
"Zh = {}; Zdr = {}; Kdp = {}; Ai = {}\n",
"for name, wl in BANDS:\n",
" Zh[name], Zdr[name], Kdp[name], Ai[name] = run(wl)\n",
" print(f'{name:<3} Z_h={Zh[name]:6.2f} Z_dr={Zdr[name]:+.3f} '\n",
" f'K_dp={Kdp[name]:+.3f} A_i={Ai[name]:.4f}')\n",
),
md("## Plot the frequency response\n"),
code(
"names = [b[0] for b in BANDS]\n",
"wls = np.array([b[1] for b in BANDS])\n",
"fig, axes = plt.subplots(2, 2, figsize=(9, 6))\n",
"axes[0, 0].semilogx(wls, [Zh[n] for n in names], 'o-'); axes[0, 0].set_ylabel('Z_h [dBZ]')\n",
"axes[0, 1].semilogx(wls, [Zdr[n] for n in names], 'o-'); axes[0, 1].set_ylabel('Z_dr [dB]')\n",
"axes[1, 0].semilogx(wls, [Kdp[n] for n in names], 'o-'); axes[1, 0].set_ylabel('K_dp [°/km]')\n",
"axes[1, 1].semilogx(wls, [Ai[n] for n in names], 'o-'); axes[1, 1].set_ylabel('A_i [dB/km]')\n",
"for ax, label in zip(axes.flat, names*4):\n",
" ax.set_xlabel('wavelength [mm]')\n",
" ax.grid(True, alpha=0.3, which='both')\n",
" for name, wl in BANDS:\n",
" ax.axvline(wl, color='k', alpha=0.1)\n",
"fig.suptitle('Moderate convective rain (D0=1.5 mm): frequency response')\n",
"fig.tight_layout();\n",
),
]
NB06 = [
md(
"# Tutorial 06 — Hydrometeor mixtures at C-band\n",
"\n",
"Real radar volumes often contain more than one species — rain with\n",
"melting aggregates, graupel in convective cores, pristine crystals\n",
"with aggregates in stratiform ice. The combined polarimetric\n",
"signature is the incoherent sum of the per-species amplitude (S) and\n",
"phase (Z) matrices. The *non*-linear observables (Z_dr, ρ_hv, δ_hv)\n",
"cannot be averaged from per-species values; they must be recomputed\n",
"from the summed matrices.\n",
"\n",
"This tutorial demonstrates how the **mixture ρ_hv drops below\n",
"unity** whenever two populations contribute comparable power and\n",
"carry different polarimetric fingerprints — even when each\n",
"individual species has ρ_hv ≈ 1. Three pairings illustrate the\n",
"point:\n",
"\n",
"1. **rain + light snow** (baseline) — rain dominates Z_h, ρ_hv barely\n",
" moves.\n",
"2. **rain + heavy snow** — snow concentration bumped so its Z_h\n",
" contribution matches rain's. The Z_dr mismatch (rain is positive,\n",
" snow is near zero with wide canting) drags ρ_hv down.\n",
"3. **rain + graupel** — near-spherical graupel with ρ = 0.4 g/cm³,\n",
" D_max = 8 mm and wide Gaussian canting (σ = 40°). C-band\n",
" non-Rayleigh resonance at the large-D tail adds a small δ_hv\n",
" signature and further degrades ρ_hv.\n",
"\n",
"Graupel parameters follow Ryzhkov, Zrnić, Burgess (2005, *JAMC*\n",
"44:557) and Kumjian (2013, *J. Operational Meteor.*): density\n",
"0.4 g/cm³, oblateness 0.8, tumbling canting σ ≈ 40°.\n",
),
code(
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"\n",
"from rustmatrix import (HydroMix, MixtureComponent, Scatterer,\n",
" orientation, radar, psd as rs_psd)\n",
"from rustmatrix.tmatrix_aux import (dsr_thurai_2007, geom_horiz_back,\n",
" geom_horiz_forw, K_w_sqr, wl_C)\n",
"from rustmatrix.refractive import m_w_10C, mi\n",
),
md("## Build one scatterer per species\n",
"\n",
"Each species gets its own `PSDIntegrator` and (for snow/graupel) an\n",
"orientation PDF that models tumbling. All integrators register the\n",
"same two geometries the mixture will query (backscatter for Z_h,\n",
"Z_dr, ρ_hv; forward for K_dp, A_i) and share the wavelength.\n"),
code(
"def build_rain():\n",
" s = Scatterer(wavelength=wl_C, m=m_w_10C[wl_C],\n",
" Kw_sqr=K_w_sqr[wl_C], ddelt=1e-4, ndgs=2)\n",
" integ = rs_psd.PSDIntegrator()\n",
" integ.D_max = 6.0; integ.num_points = 64\n",
" integ.axis_ratio_func = lambda D: 1.0 / dsr_thurai_2007(D)\n",
" integ.geometries = (geom_horiz_back, geom_horiz_forw)\n",
" s.psd_integrator = integ\n",
" s.psd_integrator.init_scatter_table(s)\n",
" return s\n",
"\n",
"def build_tumbling(rho, axis_ratio, D_max):\n",
" s = Scatterer(wavelength=wl_C, m=mi(wl_C, rho),\n",
" Kw_sqr=K_w_sqr[wl_C], axis_ratio=axis_ratio,\n",
" ddelt=1e-4, ndgs=2)\n",
" s.orient = orientation.orient_averaged_fixed\n",
" s.or_pdf = orientation.gaussian_pdf(std=40.0, mean=90.0)\n",
" s.n_alpha = 6; s.n_beta = 12\n",
" integ = rs_psd.PSDIntegrator()\n",
" integ.D_max = D_max; integ.num_points = 64\n",
" integ.geometries = (geom_horiz_back, geom_horiz_forw)\n",
" s.psd_integrator = integ\n",
" s.psd_integrator.init_scatter_table(s)\n",
" return s\n",
"\n",
"rain = build_rain()\n",
"snow = build_tumbling(rho=0.2, axis_ratio=0.6, D_max=8.0)\n",
"graupel = build_tumbling(rho=0.4, axis_ratio=0.8, D_max=8.0)\n",
"\n",
"rain_psd = rs_psd.GammaPSD(D0=1.5, Nw=8e3, mu=4, D_max=6.0)\n",
"snow_psd_light = rs_psd.ExponentialPSD(N0=5e3, Lambda=2.0, D_max=8.0)\n",
"snow_psd_heavy = rs_psd.ExponentialPSD(N0=1.5e5, Lambda=2.0, D_max=8.0)\n",
"graupel_psd = rs_psd.ExponentialPSD(N0=4e3, Lambda=1.4, D_max=8.0)\n",
"\n",
"rain.psd = rain_psd\n",
),
md("## Assemble the three mixtures and read the observables\n",
"\n",
"Each mixture is a list of `MixtureComponent(Scatterer, PSD)` pairs.\n",
"`HydroMix` exposes a Scatterer-shaped API so the usual `radar.*`\n",
"helpers work on it directly.\n"),
code(
"def obs(x):\n",
" x.set_geometry(geom_horiz_back)\n",
" out = dict(Zh=10*np.log10(radar.refl(x)),\n",
" Zdr=10*np.log10(radar.Zdr(x)),\n",
" rho=radar.rho_hv(x),\n",
" delta=np.degrees(radar.delta_hv(x)))\n",
" x.set_geometry(geom_horiz_forw)\n",
" out['Kdp'] = radar.Kdp(x); out['Ai'] = radar.Ai(x)\n",
" return out\n",
"\n",
"mix_lightsnow = HydroMix([\n",
" MixtureComponent(rain, rain_psd, 'rain'),\n",
" MixtureComponent(snow, snow_psd_light, 'snow'),\n",
"])\n",
"mix_heavysnow = HydroMix([\n",
" MixtureComponent(rain, rain_psd, 'rain'),\n",
" MixtureComponent(snow, snow_psd_heavy, 'snow'),\n",
"])\n",
"mix_graupel = HydroMix([\n",
" MixtureComponent(rain, rain_psd, 'rain'),\n",
" MixtureComponent(graupel, graupel_psd, 'graupel'),\n",
"])\n",
"\n",
"# Evaluate each case eagerly — the scatterer objects are shared by\n",
"# mutation (snow used by both light- and heavy-snow cases), so we read\n",
"# out observables immediately after setting each psd.\n",
"results = {}\n",
"results['rain only'] = obs(rain)\n",
"snow.psd = snow_psd_light; results['light snow only'] = obs(snow)\n",
"snow.psd = snow_psd_heavy; results['heavy snow only'] = obs(snow)\n",
"graupel.psd = graupel_psd; results['graupel only'] = obs(graupel)\n",
"results['rain + light snow'] = obs(mix_lightsnow)\n",
"results['rain + heavy snow'] = obs(mix_heavysnow)\n",
"results['rain + graupel'] = obs(mix_graupel)\n",
"\n",
"hdr = f'{\"case\":<22} {\"Z_h\":>7} {\"Z_dr\":>7} {\"rho_hv\":>9} {\"delta_hv\":>9} {\"K_dp\":>9}'\n",
"print(hdr)\n",
"print(f'{\"\":<22} {\"[dBZ]\":>7} {\"[dB]\":>7} {\"\":>9} {\"[deg]\":>9} {\"[deg/km]\":>9}')\n",
"print('-' * len(hdr))\n",
"for name, r in results.items():\n",
" print(f'{name:<22} {r[\"Zh\"]:>7.2f} {r[\"Zdr\"]:>+7.3f} '\n",
" f'{r[\"rho\"]:>9.5f} {r[\"delta\"]:>+9.4f} {r[\"Kdp\"]:>+9.4f}')\n",
),
md("## Bar-chart comparison\n",
"\n",
"The heavy-snow and graupel mixtures both push ρ_hv below the\n",
"rain-only and light-snow baselines. Z_dr drops because the\n",
"horizontally-oriented rain signal is diluted by near-isotropic\n",
"tumbling-ice contributions. K_dp rises in the heavy-snow mix because\n",
"the snow itself contributes a small forward phase shift.\n"),
code(
"mix_names = ['rain only', 'rain + light snow',\n",
" 'rain + heavy snow', 'rain + graupel']\n",
"colors = ['C0', 'C2', 'C4', 'C1']\n",
"\n",
"fig, axes = plt.subplots(1, 4, figsize=(12, 3.4))\n",
"for ax, key, ylab in zip(\n",
" axes,\n",
" ['Zh', 'Zdr', 'rho', 'Kdp'],\n",
" ['Z_h [dBZ]', 'Z_dr [dB]', 'ρ_hv', 'K_dp [°/km]'],\n",
"):\n",
" vals = [results[n][key] for n in mix_names]\n",
" ax.bar(range(len(mix_names)), vals, color=colors)\n",
" ax.set_xticks(range(len(mix_names)))\n",
" ax.set_xticklabels(mix_names, rotation=20, ha='right', fontsize=8)\n",
" ax.set_ylabel(ylab)\n",
" ax.grid(True, axis='y', alpha=0.3)\n",
"axes[2].set_ylim(min(results[n]['rho'] for n in mix_names) - 0.003, 1.0005)\n",
"fig.suptitle('C-band mixtures: heavy-snow and graupel pull ρ_hv below rain-only baseline')\n",
"fig.tight_layout();\n",
),
md("## Sweep the snow fraction\n",
"\n",
"Scale the snow PSD's N0 from 0 up to well beyond the rain-matching\n",
"value and watch Z_dr and ρ_hv interpolate between the rain-only\n",
"limit (left edge) and the snow-dominated limit (right edge). The\n",
"minimum ρ_hv sits near the crossover where Z_h_snow ≈ Z_h_rain —\n",
"exactly where the Z_dr mismatch has maximum leverage.\n"),
code(
"N0_snow_sweep = np.geomspace(1e3, 5e5, 16)\n",
"snow.psd = snow_psd_light\n",
"Zh_sw, Zdr_sw, rho_sw, Kdp_sw = [], [], [], []\n",
"for N0 in N0_snow_sweep:\n",
" psd_i = rs_psd.ExponentialPSD(N0=float(N0), Lambda=2.0, D_max=8.0)\n",
" m_sw = HydroMix([\n",
" MixtureComponent(rain, rain_psd, 'rain'),\n",
" MixtureComponent(snow, psd_i, 'snow'),\n",
" ])\n",
" o = obs(m_sw)\n",
" Zh_sw.append(o['Zh']); Zdr_sw.append(o['Zdr'])\n",
" rho_sw.append(o['rho']); Kdp_sw.append(o['Kdp'])\n",
"\n",
"fig, axes = plt.subplots(2, 2, figsize=(9, 6), sharex=True)\n",
"axes[0, 0].semilogx(N0_snow_sweep, Zh_sw, 'C0-o'); axes[0, 0].set_ylabel('Z_h [dBZ]')\n",
"axes[0, 1].semilogx(N0_snow_sweep, Zdr_sw, 'C1-o'); axes[0, 1].set_ylabel('Z_dr [dB]')\n",
"axes[1, 0].semilogx(N0_snow_sweep, rho_sw, 'C2-o'); axes[1, 0].set_ylabel('ρ_hv')\n",
"axes[1, 1].semilogx(N0_snow_sweep, Kdp_sw, 'C3-o'); axes[1, 1].set_ylabel('K_dp [°/km]')\n",
"for ax in axes.flat:\n",
" ax.set_xlabel('snow N₀ [m⁻³ mm⁻¹]')\n",
" ax.grid(True, alpha=0.3, which='both')\n",
"fig.suptitle('C-band rain + tumbling snow — sweeping snow N₀ across five decades')\n",
"fig.tight_layout();\n",
),
]
NB07 = [
md(
"# Tutorial 07 — W-band Mie Doppler spectrum in convective rain (Kollias 2002)\n",
"\n",
"Reference: Kollias, P., Albrecht, B. A., and Marks Jr., F. D., 2002.\n",
"*Cloud radar observations of vertical drafts and microphysics in\n",
"convective rain*, J. Geophys. Res., 107 (doi:10.1029/2001JD002033).\n",
"\n",
"At 94 GHz the raindrop backscattering cross-section σ_b(D) is no\n",
"longer Rayleigh — it rings through a sequence of Mie maxima and\n",
"minima. Because each drop falls at a deterministic terminal velocity\n",
"v_t(D), these Mie features map directly onto the observed Doppler\n",
"spectrum: the first Mie minimum appears at v ≈ 5.9 m/s in still air,\n",
"and any displacement of that feature from its theoretical location\n",
"is a direct measurement of the mean vertical air motion — independent\n",
"of the drop-size distribution. This notebook reproduces the paper's\n",
"Figures 1–3 and demonstrates the air-motion retrieval.\n",
),
code(
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"\n",
"from rustmatrix import Scatterer, SpectralIntegrator, radar, spectra\n",
"from rustmatrix.psd import ExponentialPSD, PSDIntegrator\n",
"from rustmatrix.refractive import m_w_20C\n",
"from rustmatrix.tmatrix_aux import (K_w_sqr, dsr_thurai_2007,\n",
" geom_vert_back, wl_W)\n",
),
md("## Figure 1 — single-drop σ_b(D) at 94 GHz\n",
"\n",
"Each drop's backscatter cross-section at vertical incidence, as a\n",
"function of its equivalent diameter D. The Rayleigh D⁶ regime holds\n",
"only below about 0.7 mm; beyond that the curve oscillates through\n",
"distinct Mie resonances.\n"),
code(
"def build_drop(D_mm, oblate=True):\n",
" axis_ratio = 1.0 / dsr_thurai_2007(D_mm) if oblate else 1.0\n",
" s = Scatterer(radius=D_mm/2, wavelength=wl_W, m=m_w_20C[wl_W],\n",
" Kw_sqr=K_w_sqr[wl_W], axis_ratio=axis_ratio,\n",
" ddelt=1e-4, ndgs=2)\n",
" s.set_geometry(geom_vert_back)\n",
" return s\n",
"\n",
"D_grid = np.linspace(0.05, 5.0, 400) # mm\n",
"sigma_b_obl = np.array([radar.radar_xsect(build_drop(D, True))\n",
" for D in D_grid])\n",
"sigma_b_sph = np.array([radar.radar_xsect(build_drop(D, False))\n",
" for D in D_grid])\n",
"\n",
"fig, ax = plt.subplots(figsize=(7, 4))\n",
"ax.semilogy(D_grid, sigma_b_obl, 'C0-', label='oblate (Thurai 2007)')\n",
"ax.semilogy(D_grid, sigma_b_sph, 'k--', label='sphere')\n",
"ax.set_xlabel('diameter D [mm]')\n",
"ax.set_ylabel(r'$\\sigma_b$ [mm²]')\n",
"ax.set_title('Figure 1 — single-drop backscatter at 94 GHz, vertical incidence')\n",
"ax.legend(); ax.grid(True, which='both', alpha=0.3)\n",
"fig.tight_layout();\n",
),
md("## Figure 3 — Mie oscillations mapped to terminal velocity\n",
"\n",
"Plotting σ_b / (π r²) versus the drop's terminal fall speed collapses\n",
"the Mie ripple structure into the velocity coordinate the radar\n",
"actually measures. Oblate drops fall slightly faster than the spheres\n",
"of the same equivalent diameter, so their first Mie minimum lands at\n",
"v ≈ 5.95 m/s instead of 5.88 m/s — a 7 cm/s shift that matters when\n",
"you're using it as a zero-air-motion fiducial.\n"),
code(
"fall = spectra.fall_speed.atlas_srivastava_sekhon_1973\n",
"v_t = fall(D_grid)\n",
"r_mm = D_grid / 2.0\n",
"norm = np.pi * r_mm ** 2\n",
"\n",
"def first_min(D, sigma_b):\n",
" mask = (D > 1.3) & (D < 2.2)\n",
" i = np.argmin(sigma_b[mask])\n",
" return D[mask][i], v_t[mask][i]\n",
"\n",
"D_min_obl, v_min_obl = first_min(D_grid, sigma_b_obl)\n",
"D_min_sph, v_min_sph = first_min(D_grid, sigma_b_sph)\n",
"print(f'First Mie minimum, oblate: D = {D_min_obl:.3f} mm, '\n",
" f'v = {v_min_obl:.3f} m/s')\n",
"print(f'First Mie minimum, sphere: D = {D_min_sph:.3f} mm, '\n",
" f'v = {v_min_sph:.3f} m/s')\n",
"print(f'Shift: {100 * (v_min_obl - v_min_sph):+.1f} cm/s '\n",
" '(oblate vs sphere)')\n",
"\n",
"fig, ax = plt.subplots(figsize=(7, 4))\n",
"ax.semilogy(v_t, sigma_b_obl / norm, 'C0-', label='oblate spheroids')\n",
"ax.semilogy(v_t, sigma_b_sph / norm, 'k--', label='spheres')\n",
"ax.axvline(v_min_obl, color='C3', alpha=0.6,\n",
" label=f'oblate 1st min: {v_min_obl:.2f} m/s')\n",
"ax.axvline(v_min_sph, color='C2', alpha=0.6, ls=':',\n",
" label=f'sphere 1st min: {v_min_sph:.2f} m/s')\n",
"ax.set_xlim(0, 10); ax.set_ylim(1e-5, 1e1)\n",
"ax.set_xlabel('terminal velocity $v_t(D)$ [m/s]')\n",
"ax.set_ylabel(r'$\\sigma_b / (\\pi r^2)$')\n",
"ax.set_title('Figure 3 — Mie ripples in velocity space (94 GHz)')\n",
"ax.legend(); ax.grid(True, which='both', alpha=0.3)\n",
"fig.tight_layout();\n",
),
md("## Full Doppler spectrum at 94 GHz (cf. Figure 2)\n",
"\n",
"Running an exponential rain DSD through `SpectralIntegrator` yields\n",
"the Mie-modulated Doppler spectrum the paper observes. With enough\n",
"PSD diameter samples, the Mie minima appear as *notches* in sZ_h(v)\n",
"at the velocities we just identified.\n",
"\n",
"### Note on diameter sampling\n",
"With zero turbulence and too few PSD diameters, the spectrum\n",
"*looks* spiky — every tabulated drop lands in a single velocity bin\n",
"and the bins between stay at zero. This is a sampling artifact,\n",
"**not** a physical feature: densify the diameter grid (bump\n",
"`num_points` from 64 to 256) or add a hair of turbulence and the\n",
"spikes collapse into the smooth Mie-modulated spectrum.\n"),
code(
"def build_rain_W(num_points=256):\n",
" s = Scatterer(wavelength=wl_W, m=m_w_20C[wl_W],\n",
" Kw_sqr=K_w_sqr[wl_W], ddelt=1e-4, ndgs=2)\n",
" integ = PSDIntegrator()\n",
" integ.D_max = 5.0\n",
" integ.num_points = num_points\n",
" integ.axis_ratio_func = lambda D: 1.0 / dsr_thurai_2007(D)\n",
" integ.geometries = (geom_vert_back,)\n",
" s.psd_integrator = integ\n",
" s.psd_integrator.init_scatter_table(s)\n",
" s.psd = ExponentialPSD(N0=8e3, Lambda=2.2, D_max=5.0)\n",
" return s\n",
"\n",
"rain64 = build_rain_W(num_points=64)\n",
"rain256 = build_rain_W(num_points=256)\n",
"\n",
"def run(sc, turb, w=0.0):\n",
" return SpectralIntegrator(\n",
" sc, fall_speed=fall, turbulence=turb, w=w,\n",
" v_min=-1.0, v_max=12.0, n_bins=1024,\n",
" geometry_backscatter=geom_vert_back,\n",
" ).run()\n",
"\n",
"r_sparse = run(rain64, spectra.NoTurbulence())\n",
"r_dense = run(rain256, spectra.NoTurbulence())\n",
"r_turb = run(rain256, spectra.GaussianTurbulence(0.1))\n",
"\n",
"def dBZ(x):\n",
" return 10 * np.log10(np.maximum(x, 1e-10))\n",
"\n",
"fig, ax = plt.subplots(figsize=(8, 4))\n",
"ax.plot(r_sparse.v, dBZ(r_sparse.sZ_h), 'C1-', alpha=0.6,\n",
" label='num_points=64, no turb (spiky — sampling artifact)')\n",
"ax.plot(r_dense.v, dBZ(r_dense.sZ_h), 'C0-',\n",
" label='num_points=256, no turb (Mie notches resolved)')\n",
"ax.plot(r_turb.v, dBZ(r_turb.sZ_h), 'C3-', alpha=0.8,\n",
" label='num_points=256, σ_t=0.1 m/s')\n",
"ax.axvline(v_min_obl, color='k', ls=':', alpha=0.6,\n",
" label=f'1st Mie min ({v_min_obl:.2f} m/s)')\n",
"ax.set_xlim(0, 10); ax.set_ylim(-50, 25)\n",
"ax.set_xlabel('Doppler velocity v [m/s]')\n",
"ax.set_ylabel('sZ_h [dBZ / (m/s)]')\n",
"ax.set_title('W-band Doppler spectrum — Mie notches and sampling artifact')\n",
"ax.legend(fontsize=9); ax.grid(True, alpha=0.3)\n",
"fig.tight_layout();\n",
),
md("## Air-motion retrieval\n",
"\n",
"Simulate a 1 m/s downward air motion (positive *w* under our\n",
"convention): the whole spectrum shifts by *w*, so the first Mie\n",
"minimum moves from 5.95 m/s to ~6.95 m/s. Subtracting the theoretical\n",
"still-air Mie-minimum location from the observed one recovers the air\n",
"motion — this is the Kollias 2002 technique in a nutshell, and it\n",
"doesn't depend on knowing the DSD.\n"),
code(
"r_w = run(rain256, spectra.GaussianTurbulence(0.1), w=1.0)\n",
"\n",
"band = (r_w.v > v_min_obl) & (r_w.v < v_min_obl + 2.5)\n",
"v_obs = r_w.v[band][int(np.argmin(r_w.sZ_h[band]))]\n",
"w_retrieved = v_obs - v_min_obl\n",
"print(f'Prescribed air motion w = +1.00 m/s (downward)')\n",
"print(f'Observed 1st Mie min v = {v_obs:.2f} m/s')\n",
"print(f'Still-air 1st Mie min v = {v_min_obl:.2f} m/s')\n",
"print(f'Retrieved w = {w_retrieved:+.2f} m/s')\n",
"\n",
"fig, ax = plt.subplots(figsize=(8, 4))\n",
"ax.plot(r_turb.v, dBZ(r_turb.sZ_h), 'C0-', label='w = 0')\n",
"ax.plot(r_w.v, dBZ(r_w.sZ_h), 'C3-', label='w = +1 m/s')\n",
"ax.axvline(v_min_obl, color='k', ls=':', alpha=0.6,\n",
" label=f'still-air min ({v_min_obl:.2f})')\n",
"ax.axvline(v_obs, color='C3', ls=':', alpha=0.6,\n",
" label=f'shifted min ({v_obs:.2f})')\n",
"ax.set_xlim(0, 10); ax.set_ylim(-50, 25)\n",
"ax.set_xlabel('Doppler velocity v [m/s]')\n",
"ax.set_ylabel('sZ_h [dBZ / (m/s)]')\n",
"ax.set_title('Kollias 2002 air-motion retrieval: Mie minimum as a fiducial')\n",
"ax.legend(fontsize=9); ax.grid(True, alpha=0.3)\n",
"fig.tight_layout();\n",
),
]
NB08 = [
md(
"# Tutorial 08 — Dual-frequency non-Rayleigh snowfall spectra (Billault-Roux 2023)\n",
"\n",
"Reference: Billault-Roux, A.-C., Ghiggi, G., Jaffeux, L., Martini, A.,\n",
"Viltard, N., and Berne, A., 2023. *Dual-frequency spectral radar\n",
"retrieval of snowfall microphysics: a physics-driven deep-learning\n",
"approach*, Atmos. Meas. Tech., 16, 911–931\n",
"(doi:10.5194/amt-16-911-2023).\n",
"\n",
"The paper relies on a simple but powerful asymmetry: at cloud-radar\n",
"frequencies small, slow-falling snow particles are still Rayleigh, so\n",
"their X-band and W-band reflectivities agree — but the large,\n",
"fast-falling ones are non-Rayleigh at W-band, so their W-band\n",
"spectral reflectivity is *smaller* than X-band. The **spectral\n",
"dual-wavelength ratio** sDWR(v) = 10·log₁₀(sZ_X / sZ_W) stays near 0\n",
"at low velocities and rises to several dB at high velocities — a\n",
"direct fingerprint of the large-particle size distribution tail,\n",
"untangled from turbulence and wind offsets.\n",
"\n",
"This notebook reproduces that signature by running the same snow\n",
"scatterer, PSD, fall-speed, and turbulence through\n",
"`SpectralIntegrator` at X-band (9.5 GHz) and W-band (94 GHz).\n",
),
code(
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"\n",
"from rustmatrix import Scatterer, SpectralIntegrator, radar, spectra\n",
"from rustmatrix.psd import ExponentialPSD, PSDIntegrator\n",
"from rustmatrix.refractive import mi\n",
"from rustmatrix.tmatrix_aux import (K_w_sqr, geom_vert_back,\n",
" wl_X, wl_W)\n",
),
md("## Build matched snow scatterers at X-band and W-band\n",
"\n",
"Both scatterers share the same oblate low-density aggregate habit,\n",
"PSD, and velocity grid. The *only* change across the two runs is the\n",
"radar wavelength.\n",
"\n",
"PSD and habit parameters are grounded in the ICE-GENESIS 23 January\n",
"2021 case of Billault-Roux et al. 2023 (Fig. 5 snowfall layer):\n",
"\n",
"* oblate aggregates, ρ_ice = 0.1 g/cm³ (low-density, typical of\n",
" mixed-habit aggregates), axis ratio 0.6, D_max = 5 mm\n",
"* exponential PSD, N₀ = 2×10⁴ m⁻³ mm⁻¹, Λ = 2.5 mm⁻¹\n",
"* median-volume diameter D₀ = 3.67/Λ ≈ 1.5 mm\n",
"* IWC = (π ρ_ice / Λ⁴) · N₀ ≈ 0.16 g/m³ — a moderate aggregation\n",
" layer, well within the synthetic training range the authors drew\n",
" from their MASCDB disdrometer database.\n",
"\n",
"These numbers land bulk Z_h(X) around 17 dBZ (moderate snowfall,\n",
"consistent with the paper's Fig. 5 observations at ~15:10 UTC) and a\n",
"non-Rayleigh-driven bulk DWR on the order of 15–20 dB — *not* the\n",
"50 dBZ / 25 dB outlier you get if you blindly pick Λ = 0.8 mm⁻¹ and\n",
"D_max = 10 mm (D₀ ≈ 5 mm, IWC > 10 g/m³ — a reflectivity-saturating\n",
"deep convective snow cell, not an aggregation layer).\n"),
code(
"RHO_ICE = 0.1\n",
"AXIS_RATIO = 0.6\n",
"D_MAX = 5.0\n",
"N0 = 2e4\n",
"LAMBDA = 2.5\n",
"\n",
"def build_snow(wl):\n",
" s = Scatterer(wavelength=wl, m=mi(wl, RHO_ICE),\n",
" Kw_sqr=K_w_sqr[wl], axis_ratio=AXIS_RATIO,\n",
" ddelt=1e-4, ndgs=2)\n",
" integ = PSDIntegrator()\n",
" integ.D_max = D_MAX\n",
" integ.num_points = 256\n",
" integ.geometries = (geom_vert_back,)\n",
" s.psd_integrator = integ\n",
" s.psd_integrator.init_scatter_table(s)\n",
" s.psd = ExponentialPSD(N0=N0, Lambda=LAMBDA, D_max=D_MAX)\n",
" return s\n",
"\n",
"snow_X = build_snow(wl_X)\n",
"snow_W = build_snow(wl_W)\n",
),
md("## σ_b(D) at X-band vs W-band — the non-Rayleigh onset\n",
"\n",
"Up to about 2 mm equivalent diameter the two σ_b(D) curves are\n",
"essentially parallel D⁶ power laws (Rayleigh). Above ~3 mm the\n",
"W-band curve rolls over and develops Mie structure while the X-band\n",
"curve keeps rising — this is the imprint the dual-frequency spectrum\n",
"will show up as a velocity-resolved sDWR.\n"),
code(
"def sigma_b(wl, D_mm):\n",
" s = Scatterer(radius=D_mm/2, wavelength=wl,\n",
" m=mi(wl, RHO_ICE), Kw_sqr=K_w_sqr[wl],\n",
" axis_ratio=AXIS_RATIO, ddelt=1e-4, ndgs=2)\n",
" s.set_geometry(geom_vert_back)\n",
" return radar.radar_xsect(s)\n",
"\n",
"D = np.linspace(0.1, D_MAX, 150)\n",
"sb_X = np.array([sigma_b(wl_X, d) for d in D])\n",
"sb_W = np.array([sigma_b(wl_W, d) for d in D])\n",
"\n",
"fig, ax = plt.subplots(figsize=(7, 4))\n",
"ax.loglog(D, sb_X, 'C0-', label='X-band (9 GHz)')\n",
"ax.loglog(D, sb_W, 'C3-', label='W-band (94 GHz)')\n",
"ax.set_xlabel('equivalent diameter D [mm]')\n",
"ax.set_ylabel(r'$\\sigma_b$ [mm²]')\n",
"ax.set_title('Snow backscatter at vertical incidence')\n",
"ax.legend(); ax.grid(True, which='both', alpha=0.3)\n",
"fig.tight_layout();\n",
),
md("## Dual-frequency Doppler spectra\n",
"\n",
"Same fall-speed model (Locatelli–Hobbs aggregates), same turbulence,\n",
"same velocity grid — only the radar wavelength changes.\n"),
code(
"fall = spectra.fall_speed.locatelli_hobbs_1974_aggregates\n",
"turb = spectra.GaussianTurbulence(0.2)\n",
"V_MIN, V_MAX = -2.0, 4.0\n",
"N_BINS = 1024\n",
"\n",
"def run(sc):\n",
" return SpectralIntegrator(\n",
" sc, fall_speed=fall, turbulence=turb,\n",
" v_min=V_MIN, v_max=V_MAX, n_bins=N_BINS,\n",
" geometry_backscatter=geom_vert_back,\n",
" ).run()\n",
"\n",
"r_X = run(snow_X)\n",
"r_W = run(snow_W)\n",
"\n",
"def dBZ(x):\n",
" return 10 * np.log10(np.maximum(x, 1e-12))\n",
"\n",
"fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(8, 6), sharex=True)\n",
"ax1.plot(r_X.v, dBZ(r_X.sZ_h), 'C0-', label='X-band')\n",
"ax1.plot(r_W.v, dBZ(r_W.sZ_h), 'C3-', label='W-band')\n",
"ax1.set_ylabel('sZ_h [dBZ / (m/s)]')\n",
"ax1.set_title('Dual-frequency spectra of snowfall — vertical pointing')\n",
"ax1.legend(); ax1.grid(True, alpha=0.3)\n",
"\n",
"sDWR = dBZ(r_X.sZ_h) - dBZ(r_W.sZ_h)\n",
"good = (r_X.sZ_h > 1e-8) & (r_W.sZ_h > 1e-10)\n",
"ax2.plot(r_X.v[good], sDWR[good], 'C2-')\n",
"ax2.axhline(0.0, color='k', ls=':', alpha=0.6)\n",
"ax2.set_xlabel('Doppler velocity v [m/s]')\n",
"ax2.set_ylabel('sDWR$_{X-W}$ [dB]')\n",
"ax2.set_title('Spectral dual-wavelength ratio: Rayleigh at low v, '\n",
" 'non-Rayleigh at high v')\n",
"ax2.grid(True, alpha=0.3)\n",
"ax2.set_xlim(V_MIN, V_MAX)\n",
"fig.tight_layout();\n",
),
md("## Reading the signature\n",
"\n",
"At velocities below ~0.5 m/s (small aggregates, well inside the\n",
"Rayleigh regime at both bands) sDWR is essentially zero. It rises\n",
"monotonically with v as the fastest-falling particles move into the\n",
"W-band non-Rayleigh regime while remaining Rayleigh at X-band. The\n",
"magnitude of the rise at a given velocity is a direct proxy for the\n",
"*size* of the drops there — which is exactly the lever Billault-Roux\n",
"et al. 2023 use to retrieve the PSD tail, with deep learning\n",
"shouldering the joint dependence on turbulence and shape.\n"),
code(
"print(f'bulk Z_h (X-band) = '\n",
" f'{10*np.log10(np.trapezoid(r_X.sZ_h, r_X.v)):.2f} dBZ')\n",
"print(f'bulk Z_h (W-band) = '\n",
" f'{10*np.log10(np.trapezoid(r_W.sZ_h, r_W.v)):.2f} dBZ')\n",
"print(f'bulk DWR = '\n",
" f'{10*np.log10(np.trapezoid(r_X.sZ_h, r_X.v) / np.trapezoid(r_W.sZ_h, r_W.v)):.2f} dB')\n",
"print()\n",
"v_samples = [0.3, 0.5, 0.8, 1.0, 1.3, 1.6]\n",
"print('sDWR(v) at selected velocities:')\n",
"for vs in v_samples:\n",
" i = int(np.argmin(np.abs(r_X.v - vs)))\n",
" print(f' v = {r_X.v[i]:.2f} m/s sDWR = '\n",
" f'{dBZ(r_X.sZ_h[i]) - dBZ(r_W.sZ_h[i]):+.2f} dB')\n",
),
]
NB09 = [
md(
"# Tutorial 09 — Reproducing Zhu, Kollias, Yang (2023): particle-inertia spectra\n",
"\n",
"Reference: Zhu, Z., Kollias, P., and Yang, F., 2023. *Particle Inertia\n",
"Effects on Radar Doppler Spectra Simulation*, Atmos. Meas. Tech.\n",
"Discuss. MATLAB simulator source: https://zenodo.org/records/7897981\n",
"\n",
"The conventional way to get a turbulent radar Doppler spectrum from a\n",
"drop-size distribution is to convolve the drop reflectivity spectrum\n",
"with a single Gaussian of width σ_air. Zhu 2023 argues — and shows\n",
"with a per-drop drag-ODE simulator — that this is wrong: small drops\n",
"are passive tracers, but large drops are ballistic and under-respond\n",
"to the small-scale eddies of the wind field. The broadening kernel is\n",
"*diameter-dependent*, σ_t(D), and narrows toward the large-drop end.\n",
"\n",
"rustmatrix ships an analytical Stokes-number low-pass response,\n",
"`InertialZeng2023`, implementing this physics family (the Zeng et al.\n",
"2023 companion paper):\n",
"\n",
"$$\\sigma_t(D) = \\frac{\\sigma_{\\rm air}}{\\sqrt{1 + \\mathrm{St}(D)^2}}, \\qquad \\mathrm{St} = \\tau_p(D) / \\tau_{\\rm eddy}.$$\n",
"\n",
"This notebook reproduces Zhu 2023's configuration — W-band, ±12 m/s\n",
"Nyquist, 1024 bins, SNR = 40 dB, exponential warm-rain PSD — and\n",
"compares conventional Gaussian broadening to the inertia-aware kernel.\n",
),
code(
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"\n",
"from rustmatrix import Scatterer, SpectralIntegrator, spectra\n",
"from rustmatrix.psd import ExponentialPSD, PSDIntegrator\n",
"from rustmatrix.refractive import m_w_10C\n",
"from rustmatrix.tmatrix_aux import K_w_sqr, geom_vert_back, wl_W\n",
"\n",
"# --- Zhu 2023 configuration ---\n",
"V_MIN, V_MAX = -12.0, 12.0\n",
"NFFT = 1024\n",
"SNR_DB = 40.0\n",
"Z_DBZ = 20.0\n",
"SIGMA_AIR = 0.289 # std of uniform [-0.5, 0.5] wind field\n",
"EPSILON = 1e-2 # dissipation rate (m²/s³)\n",
"L_O = 0.1 # integral length scale (m)\n",
),
md("## Build the W-band rain scatterer and PSD\n",
"\n",
"Zhu 2023's MATLAB code uses N(D) = N₀·exp(shape·D_cm)·1e4 m⁻³ µm⁻¹\n",
"with N₀ = 0.08 and shape = -20. In rustmatrix's mm-based units this\n",
"is N(D_mm) = 8×10⁵ · exp(-2·D_mm) m⁻³ mm⁻¹, truncated at D_max=4 mm.\n"),
code(
"s = Scatterer(wavelength=wl_W, m=m_w_10C[wl_W],\n",
" Kw_sqr=K_w_sqr[wl_W], ddelt=1e-4, ndgs=2)\n",
"integ = PSDIntegrator()\n",
"integ.D_max = 4.0; integ.num_points = 128\n",
"integ.geometries = (geom_vert_back,)\n",
"s.psd_integrator = integ\n",
"s.psd_integrator.init_scatter_table(s)\n",
"s.psd = ExponentialPSD(N0=8e5, Lambda=2.0, D_max=4.0)\n",
),
md("## The inertial reduction factor σ_t(D) / σ_air\n",
"\n",
"At fixed σ_air and ε, the reduction factor depends only on D through\n",
"the particle relaxation time τ_p(D) = v_t(D)/g. This plot is the\n",
"physical crux of Zhu 2023 — large drops experience less-effective\n",
"turbulence broadening than small drops.\n"),
code(
"zeng = spectra.InertialZeng2023(sigma_air=SIGMA_AIR,\n",
" epsilon=EPSILON, L_o=L_O)\n",
"D_plot = np.linspace(0.1, 4.0, 200)\n",
"fig, ax = plt.subplots(figsize=(7, 4))\n",
"ax.plot(D_plot, zeng(D_plot) / SIGMA_AIR, 'C0-', lw=2)\n",
"ax.axhline(1.0, color='k', ls='--', alpha=0.5,\n",
" label='conventional: σ_t / σ_air = 1')\n",
"ax.set_xlabel('drop diameter D [mm]')\n",
"ax.set_ylabel('σ_t(D) / σ_air')\n",
"ax.set_title(f'Inertial reduction (σ_air={SIGMA_AIR}, '\n",
" f'ε={EPSILON}, L_o={L_O})')\n",
"ax.legend(); ax.grid(True, alpha=0.3)\n",
"ax.set_ylim(0.0, 1.05)\n",
"fig.tight_layout();\n",
),
md("## Run two spectra: conventional vs. inertia-aware\n",
"\n",
"Both use the same fall-speed model, velocity grid, and noise floor.\n",
"The only difference is the turbulence kernel — constant σ_air vs.\n",
"σ_t(D) from Zeng 2023.\n"),
code(
"# Receiver noise: SNR = 40 dB over 20 dBZ reference.\n",
"noise = 10 ** (Z_DBZ / 10.0) / 10 ** (SNR_DB / 10.0)\n",
"\n",
"common = dict(\n",
" fall_speed=spectra.fall_speed.atlas_srivastava_sekhon_1973,\n",
" v_min=V_MIN, v_max=V_MAX, n_bins=NFFT,\n",
" geometry_backscatter=geom_vert_back,\n",
" noise=noise,\n",
")\n",
"\n",
"conv = SpectralIntegrator(s,\n",
" turbulence=spectra.GaussianTurbulence(SIGMA_AIR),\n",
" **common,\n",
").run()\n",
"\n",
"iner = SpectralIntegrator(s,\n",
" turbulence=spectra.InertialZeng2023(\n",
" sigma_air=SIGMA_AIR, epsilon=EPSILON, L_o=L_O),\n",
" **common,\n",
").run()\n",
"\n",
"print(f'Conventional integrated Z_h = '\n",
" f'{10*np.log10(np.trapezoid(conv.sZ_h, conv.v)):.2f} dBZ')\n",
"print(f'Inertia-aware integrated Z_h = '\n",
" f'{10*np.log10(np.trapezoid(iner.sZ_h, iner.v)):.2f} dBZ')\n",
"print(f'Noise floor = '\n",
" f'{noise:.4f} mm⁶/m³ ('\n",
" f'{10*np.log10(noise / (V_MAX - V_MIN)):.2f} '\n",
" f'dBZ / (m/s))')\n",
),
md("## Spectrum comparison in dBZ\n",
"\n",
"On a log scale the inertia correction visibly narrows the spectrum\n",
"near the fall-speed of the largest drops (the right-hand tail). The\n",
"total power integrates to the same reflectivity — this is a\n",
"redistribution of spectral shape, not a gain change.\n"),
code(
"def dBZ(x):\n",
" return 10 * np.log10(np.maximum(x, 1e-12))\n",
"\n",
"fig, ax = plt.subplots(figsize=(8, 4))\n",
"ax.plot(conv.v, dBZ(conv.sZ_h), 'C0-',\n",
" label='conventional Gaussian')\n",
"ax.plot(iner.v, dBZ(iner.sZ_h), 'C3-',\n",
" label='inertia-aware (Zeng 2023)')\n",
"ax.axhline(dBZ(noise / (V_MAX - V_MIN)), color='k',\n",
" ls='--', alpha=0.5, label='noise floor')\n",
"peak = max(dBZ(conv.sZ_h).max(), dBZ(iner.sZ_h).max())\n",
"ax.set_xlim(0, 10)\n",
"ax.set_ylim(-50, np.ceil(peak / 5.0) * 5.0 + 5)\n",
"ax.set_xlabel('Doppler velocity v [m/s]')\n",
"ax.set_ylabel('sZ_h [dBZ / (m/s)]')\n",
"ax.set_title('Zhu 2023 setup — W-band warm-rain Doppler spectrum')\n",
"ax.legend(); ax.grid(True, alpha=0.3)\n",
"fig.tight_layout();\n",
),
md("## Zoom on the large-drop tail\n",
"\n",
"Narrowing is most visible between v ≈ 7–9 m/s, where the fastest\n",
"drops land. Zhu 2023's Lagrangian drop-tracking simulator produces\n",
"the same qualitative signature (their Fig. 3–5).\n"),
code(
"fig, ax = plt.subplots(figsize=(8, 4))\n",
"ax.plot(conv.v, dBZ(conv.sZ_h), 'C0-', label='conventional')\n",
"ax.plot(iner.v, dBZ(iner.sZ_h), 'C3-', label='inertia-aware')\n",
"ax.set_xlim(6, 10)\n",
"tail_peak = max(dBZ(conv.sZ_h[(conv.v > 6) & (conv.v < 10)]).max(),\n",
" dBZ(iner.sZ_h[(iner.v > 6) & (iner.v < 10)]).max())\n",
"ax.set_ylim(-25, np.ceil(tail_peak / 5.0) * 5.0 + 5)\n",
"ax.set_xlabel('Doppler velocity v [m/s]')\n",
"ax.set_ylabel('sZ_h [dBZ / (m/s)]')\n",
"ax.set_title('Large-drop tail: inertia narrows the spectrum')\n",
"ax.legend(); ax.grid(True, alpha=0.3)\n",
"fig.tight_layout();\n",
),
]
NB10 = [
md(
"# Tutorial 10 — Supercooled liquid water vs. snow at cloud-radar frequencies\n",
"\n",
"Reference: Billault-Roux, A.-C., Georgakaki, P., Grazioli, J.,\n",
"Romanens, G., Sotiropoulou, G., Phillips, V., Nenes, A., and Berne,\n",
"A., 2023. *Distinct secondary ice production processes observed in\n",
"radar Doppler spectra: insights from a case study*, Atmos. Chem.\n",
"Phys., 23, 10207–10234 (doi:10.5194/acp-23-10207-2023).\n",
"\n",
"Mixed-phase cloud volumes contain supercooled liquid droplets (SLW)\n",
"and ice particles side by side. The two populations have completely\n",
"different scattering and fall-speed signatures:\n",
"\n",
"* **SLW**: small (≲200 µm), spherical, cold refractive index of\n",
" water, falls at a few cm/s. Rayleigh at both X- and W-band ⇒ very\n",
" low Z_h, Z_dr ≈ 0 dB, no polarimetric signature.\n",
"* **Snow**: millimetre-scale oblate low-density aggregates, falls at\n",
" ~0.5–1 m/s. Moderately non-Rayleigh at W-band ⇒ higher Z_h, small\n",
" positive Z_dr from the habit anisotropy.\n",
"\n",
"This notebook builds dedicated scatterers for each population, shows\n",
"their σ_b(D) and bulk radar observables at X and W band, and then\n",
"combines them in a `HydroMix` to produce the bimodal Doppler spectrum\n",
"that motivates the Billault-Roux et al. case-study analysis.\n",
),
code(
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"\n",
"from rustmatrix import (HydroMix, MixtureComponent, Scatterer,\n",
" SpectralIntegrator, radar, spectra)\n",
"from rustmatrix.psd import ExponentialPSD, GammaPSD, PSDIntegrator\n",
"from rustmatrix.refractive import m_w_0C, mi\n",
"from rustmatrix.tmatrix_aux import (K_w_sqr, geom_vert_back,\n",
" wl_X, wl_W)\n",
),
md("## Build a supercooled-liquid-water scatterer\n",
"\n",
"Small spherical drops (D in 10–200 µm), refractive index at 0°C,\n",
"log-normal-ish gamma PSD centred around 30 µm. We use mm units\n",
"throughout to match `rustmatrix` conventions.\n"),
code(
"SLW_DMAX = 0.2 # mm\n",
"\n",
"def build_slw(wl):\n",
" s = Scatterer(wavelength=wl, m=m_w_0C[wl],\n",
" Kw_sqr=K_w_sqr[wl], axis_ratio=1.0,\n",
" ddelt=1e-4, ndgs=2)\n",
" integ = PSDIntegrator()\n",
" integ.D_max = SLW_DMAX\n",
" integ.num_points = 128\n",
" integ.geometries = (geom_vert_back,)\n",
" s.psd_integrator = integ\n",
" s.psd_integrator.init_scatter_table(s)\n",
" # Gamma PSD: D0 = 0.03 mm (30 µm), Nw = 1e11 m⁻³ mm⁻¹, μ = 4\n",
" s.psd = GammaPSD(D0=0.03, Nw=1e11, mu=4, D_max=SLW_DMAX)\n",
" return s\n",
"\n",
"slw_X = build_slw(wl_X)\n",
"slw_W = build_slw(wl_W)\n",
),
md("## Build a snow scatterer (low-density oblate aggregates)\n"),
code(
"SNOW_DMAX = 8.0\n",
"RHO_SNOW = 0.2\n",
"AXIS_SNOW = 0.6\n",
"\n",
"def build_snow(wl):\n",
" s = Scatterer(wavelength=wl, m=mi(wl, RHO_SNOW),\n",
" Kw_sqr=K_w_sqr[wl], axis_ratio=AXIS_SNOW,\n",
" ddelt=1e-4, ndgs=2)\n",
" integ = PSDIntegrator()\n",
" integ.D_max = SNOW_DMAX\n",
" integ.num_points = 256\n",
" integ.geometries = (geom_vert_back,)\n",
" s.psd_integrator = integ\n",
" s.psd_integrator.init_scatter_table(s)\n",
" s.psd = ExponentialPSD(N0=5e3, Lambda=1.0, D_max=SNOW_DMAX)\n",
" return s\n",
"\n",
"snow_X = build_snow(wl_X)\n",
"snow_W = build_snow(wl_W)\n",
),
md("## σ_b(D) — Rayleigh vs. Mie across the two populations\n"),
code(
"def sigma_b(wl, D_mm, m, axis_ratio):\n",
" s = Scatterer(radius=D_mm/2, wavelength=wl, m=m,\n",
" Kw_sqr=K_w_sqr[wl], axis_ratio=axis_ratio,\n",
" ddelt=1e-4, ndgs=2)\n",
" s.set_geometry(geom_vert_back)\n",
" return radar.radar_xsect(s)\n",
"\n",
"D_slw = np.linspace(0.005, SLW_DMAX, 60)\n",
"D_snow = np.linspace(0.1, SNOW_DMAX, 200)\n",
"\n",
"sb_slw_X = np.array([sigma_b(wl_X, d, m_w_0C[wl_X], 1.0) for d in D_slw])\n",
"sb_slw_W = np.array([sigma_b(wl_W, d, m_w_0C[wl_W], 1.0) for d in D_slw])\n",
"sb_snow_X = np.array([sigma_b(wl_X, d, mi(wl_X, RHO_SNOW), AXIS_SNOW)\n",
" for d in D_snow])\n",
"sb_snow_W = np.array([sigma_b(wl_W, d, mi(wl_W, RHO_SNOW), AXIS_SNOW)\n",
" for d in D_snow])\n",
"\n",
"fig, ax = plt.subplots(figsize=(7.5, 4.5))\n",
"ax.loglog(D_slw, sb_slw_X, 'C0-', label='SLW, X-band')\n",
"ax.loglog(D_slw, sb_slw_W, 'C0--', label='SLW, W-band')\n",
"ax.loglog(D_snow, sb_snow_X, 'C3-', label='snow, X-band')\n",
"ax.loglog(D_snow, sb_snow_W, 'C3--', label='snow, W-band')\n",
"ax.set_xlabel('diameter D [mm]')\n",
"ax.set_ylabel(r'$\\sigma_b$ [mm²]')\n",
"ax.set_title('Single-particle backscatter: SLW droplets vs. snow aggregates')\n",
"ax.legend(fontsize=9); ax.grid(True, which='both', alpha=0.3)\n",
"fig.tight_layout();\n",
),
md("## Bulk radar observables at X and W band\n",
"\n",
"Z_h(SLW) is ~40+ dB below Z_h(snow) even though the droplet\n",
"concentration is high — the D⁶ Rayleigh penalty at 30-µm scale kills\n",
"the signal. Z_dr is zero for SLW (spheres) but positive for snow\n",
"(oblate). Dual-wavelength ratio DWR = Z_X − Z_W is near zero for SLW\n",
"and positive for snow (non-Rayleigh at W-band).\n"),
code(
"def bulk(s):\n",
" s.set_geometry(geom_vert_back)\n",
" return 10 * np.log10(radar.refl(s))\n",
"\n",
"print(f'{\"\":<10} {\"Z_h X [dBZ]\":>12} {\"Z_h W [dBZ]\":>12} '\n",
" f'{\"DWR [dB]\":>10}')\n",
"for name, sX, sW in [('SLW', slw_X, slw_W),\n",
" ('snow', snow_X, snow_W)]:\n",
" zX = bulk(sX); zW = bulk(sW)\n",
" print(f'{name:<10} {zX:>12.2f} {zW:>12.2f} {zX - zW:>10.2f}')\n",
),
md("## Bimodal Doppler spectrum of the SLW + snow mixture\n",
"\n",
"A `HydroMix` with per-component fall-speed and turbulence captures\n",
"the two populations cleanly: SLW pins to v ≈ 0 with a narrow spread,\n",
"snow sits near 0.5–1.0 m/s with a broader tail. Plotted in dBZ space\n",
"the two modes are clearly separable even with realistic noise.\n"),
code(
"# Use W-band for the spectrum — best velocity-wise resolution of the\n",
"# SLW mode, and non-Rayleigh snow signature is most evident here.\n",
"mix = HydroMix([\n",
" MixtureComponent(slw_W, slw_W.psd, 'slw'),\n",
" MixtureComponent(snow_W, snow_W.psd, 'snow'),\n",
"])\n",
"\n",
"integ = SpectralIntegrator(\n",
" mix,\n",
" component_kinematics={\n",
" 'slw': (lambda D: 0.003 * (D / 0.01) ** 2,\n",
" spectra.GaussianTurbulence(0.1)),\n",
" 'snow': (spectra.fall_speed.locatelli_hobbs_1974_aggregates,\n",
" spectra.GaussianTurbulence(0.2)),\n",
" },\n",
" v_min=-1.0, v_max=3.0, n_bins=1024,\n",
" geometry_backscatter=geom_vert_back,\n",
" noise='realistic',\n",
").run()\n",
"\n",
"def dBZ(x):\n",
" return 10 * np.log10(np.maximum(x, 1e-12))\n",
"\n",
"# Run each population on its own (no noise) so each mode's true peak\n",
"# is unambiguous. The HydroMix spectrum is dominated by snow (~50 dB\n",
"# louder than SLW), so looking for the SLW peak on the combined,\n",
"# noise-added spectrum would just pick up the snow tail crossing\n",
"# through the near-zero velocity range.\n",
"slw_only = SpectralIntegrator(\n",
" slw_W,\n",
" fall_speed=lambda D: 0.003 * (D / 0.01) ** 2,\n",
" turbulence=spectra.GaussianTurbulence(0.1),\n",
" v_min=-1.0, v_max=3.0, n_bins=1024,\n",
" geometry_backscatter=geom_vert_back,\n",
").run()\n",
"snow_only = SpectralIntegrator(\n",
" snow_W,\n",
" fall_speed=spectra.fall_speed.locatelli_hobbs_1974_aggregates,\n",
" turbulence=spectra.GaussianTurbulence(0.2),\n",
" v_min=-1.0, v_max=3.0, n_bins=1024,\n",
" geometry_backscatter=geom_vert_back,\n",
").run()\n",
"\n",
"v_slw = slw_only.v[int(np.argmax(slw_only.sZ_h))]\n",
"v_snow = snow_only.v[int(np.argmax(snow_only.sZ_h))]\n",
"print(f'SLW mode peak v = {v_slw:+.3f} m/s')\n",
"print(f'snow mode peak v = {v_snow:+.3f} m/s')\n",
"\n",
"fig, ax = plt.subplots(figsize=(8, 4))\n",
"ax.plot(slw_only.v, dBZ(slw_only.sZ_h), 'C2-', lw=1.3, label='SLW only')\n",
"ax.plot(snow_only.v, dBZ(snow_only.sZ_h), 'C3-', lw=1.3, label='snow only')\n",
"ax.plot(integ.v, dBZ(integ.sZ_h), 'C0-', lw=2.0,\n",
" label='SLW + snow (HydroMix, with noise)')\n",
"ax.axvline(v_slw, color='C2', alpha=0.4, ls=':',\n",
" label=f'SLW peak ({v_slw:+.2f})')\n",
"ax.axvline(v_snow, color='C3', alpha=0.4, ls=':',\n",
" label=f'snow peak ({v_snow:+.2f})')\n",
"ax.axhline(dBZ(integ.noise_h / (integ.v[-1] - integ.v[0])),\n",
" color='k', ls='--', alpha=0.5, label='noise floor')\n",
"ax.set_xlabel('Doppler velocity v [m/s]')\n",
"ax.set_ylabel('sZ_h [dBZ / (m/s)]')\n",
"ax.set_title('W-band SLW vs. snow vs. combined Doppler spectrum')\n",
"ax.legend(fontsize=9, loc='upper right'); ax.grid(True, alpha=0.3)\n",
"fig.tight_layout();\n",
),
md("## Takeaway\n",
"\n",
"The Billault-Roux et al. 2023 case study relies on the *spectral*\n",
"separability demonstrated here: a narrow near-zero Doppler peak is\n",
"SLW, a broader fall-velocity peak is ice; the polarimetric spectral\n",
"variables further distinguish columnar (high SLDR) from planar and\n",
"aggregated ice. With the two populations built as independent\n",
"rustmatrix scatterers you can plug them into `HydroMix` and recover\n",
"bulk radar observables that match what a vertically-pointing radar\n",
"would measure, or analyse the spectral signatures by looking at the\n",
"individual modes.\n"),
]
NB11 = [
md(
"# Tutorial 11 — Dual-frequency radar signatures across hydrometeor classes\n",
"\n",
"Reference: Honeyager, R., 2013. *Investigating the use of the T-matrix\n",
"method as a proxy for the discrete dipole approximation*, M.S. thesis,\n",
"Florida State University.\n",
"\n",
"Honeyager's thesis argues that the T-matrix applied to a size-matched\n",
"spheroid, with a dielectric built from an ice / air mixing formula at\n",
"the right *volume fraction*, reproduces the single-scattering\n",
"properties of a geometrically-complex hydrometeor (bullet rosettes,\n",
"plates, dendrites, aggregates) computed with DDA — provided the size\n",
"parameter χ = 2π r / λ and the effective density are right.\n",
"\n",
"That framing is directly useful for dual-frequency radar work: the\n",
"whole zoo of ice habits collapses, to first order, onto a\n",
"two-parameter family (*effective density ρ_eff*, *axis ratio*) that\n",
"`rustmatrix` already supports via `refractive.mi(wl, rho)` and the\n",
"`axis_ratio` keyword. This notebook walks through four representative\n",
"hydrometeor classes:\n",
"\n",
"| class | ρ_eff [g/cm³] | axis ratio |\n",
"|----------------------|---------------|------------|\n",
"| rain | 1.00 (water) | Thurai 2007|\n",
"| low-density aggregate| 0.10 | 0.70 |\n",
"| graupel | 0.50 | 0.90 |\n",
"| high-density ice | 0.90 | 1.00 |\n",
"\n",
"and shows how each class signs itself into single-particle σ_b(D),\n",
"bulk DWR vs. D₀, and the spectral DWR profile sDWR(v).\n",
),
code(
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"\n",
"from rustmatrix import Scatterer, SpectralIntegrator, radar, spectra\n",
"from rustmatrix.psd import ExponentialPSD, PSDIntegrator\n",
"from rustmatrix.refractive import m_w_10C, mi\n",
"from rustmatrix.tmatrix_aux import (K_w_sqr, dsr_thurai_2007,\n",
" geom_vert_back, wl_X, wl_W)\n",
"\n",
"CLASSES = {\n",
" 'rain': dict(rho=1.00, axis_ratio=None), # Thurai below\n",
" 'low-ρ agg':dict(rho=0.10, axis_ratio=0.70),\n",
" 'graupel': dict(rho=0.50, axis_ratio=0.90),\n",
" 'high-ρ ice':dict(rho=0.90, axis_ratio=1.00),\n",
"}\n",
"\n",
"def refractive_index(wl, cls):\n",
" return m_w_10C[wl] if cls == 'rain' else mi(wl, CLASSES[cls]['rho'])\n",
"\n",
"def class_axis_ratio(cls, D_mm):\n",
" ar = CLASSES[cls]['axis_ratio']\n",
" return (1.0 / dsr_thurai_2007(D_mm)) if ar is None else ar\n",
"\n",
"COLORS = {'rain': 'C0', 'low-ρ agg': 'C2',\n",
" 'graupel': 'C1', 'high-ρ ice': 'C3'}\n",
),
md("## Single-particle σ_b(D) at X and W\n",
"\n",
"For each class, tabulate the single-particle backscatter cross-section\n",
"at vertical incidence across 0.2–8 mm equivalent diameter. The key\n",
"effects to watch:\n",
"\n",
"* **Rayleigh regime** — all classes follow σ_b ∝ D⁶ at small D, but\n",
" offset vertically by |K(ρ)|² differences. Low-ρ ice sits ~30 dB below\n",
" water at equal D.\n",
"* **Non-Rayleigh onset at W-band** — the first Mie minimum appears once\n",
" χ = πD/λ ≈ 1, i.e. D ≈ 1 mm at W-band. That happens for all four\n",
" classes because χ is set by size, not ρ.\n",
"* **Mie oscillation amplitude** — *this* is where ρ matters. Denser\n",
" classes (rain, high-ρ ice) show sharp Mie minima/maxima; low-density\n",
" aggregates barely oscillate because their refractive index contrast\n",
" with air is small.\n"),
code(
"def sigma_b(wl, cls, D_mm):\n",
" s = Scatterer(radius=D_mm/2, wavelength=wl,\n",
" m=refractive_index(wl, cls),\n",
" Kw_sqr=K_w_sqr[wl],\n",
" axis_ratio=class_axis_ratio(cls, D_mm),\n",
" ddelt=1e-4, ndgs=2)\n",
" s.set_geometry(geom_vert_back)\n",
" return radar.radar_xsect(s)\n",
"\n",
"D_grid = np.linspace(0.2, 8.0, 100)\n",
"sigma = {cls: {b: np.array([sigma_b(wl, cls, D) for D in D_grid])\n",
" for b, wl in [('X', wl_X), ('W', wl_W)]}\n",
" for cls in CLASSES}\n",
"\n",
"fig, (axX, axW) = plt.subplots(1, 2, figsize=(11, 4), sharey=True)\n",
"for cls in CLASSES:\n",
" axX.loglog(D_grid, sigma[cls]['X'], color=COLORS[cls], label=cls)\n",
" axW.loglog(D_grid, sigma[cls]['W'], color=COLORS[cls], label=cls)\n",
"axX.set_title('σ_b(D) at X-band (9 GHz)')\n",
"axW.set_title('σ_b(D) at W-band (94 GHz)')\n",
"for ax in (axX, axW):\n",
" ax.set_xlabel('D [mm]')\n",
" ax.grid(True, which='both', alpha=0.3)\n",
" ax.legend(fontsize=9)\n",
"axX.set_ylabel(r'$\\sigma_b$ [mm²]')\n",
"fig.tight_layout();\n",
),
md("## DWR(D) — single-particle dual-wavelength ratio\n",
"\n",
"Re-plotting as DWR(D) = 10·log₁₀(σ_b^X / σ_b^W) collapses the key\n",
"information: the *size where each class first walks out of Rayleigh*.\n",
"This is Honeyager's thesis Table 2.2 idea generalised to arbitrary\n",
"habit — the χ parameter that sets Rayleigh validity is hidden inside\n",
"ρ_eff via the refractive index contrast.\n"),
code(
"# Equivalent reflectivity normalisation: Ze = λ⁴ / (π⁵ |K_w|²) · σ_b.\n",
"# In the Rayleigh regime Ze_X = Ze_W, so DWR_{X-W} = 0 dB; any positive\n",
"# excursion is a direct signature of non-Rayleigh scattering at W-band.\n",
"def single_Ze(sig, wl, Kw2):\n",
" return wl**4 / (np.pi**5 * Kw2) * sig\n",
"\n",
"fig, ax = plt.subplots(figsize=(8, 4.5))\n",
"dwr_curves_1p = {}\n",
"for cls in CLASSES:\n",
" zeX = single_Ze(sigma[cls]['X'], wl_X, K_w_sqr[wl_X])\n",
" zeW = single_Ze(sigma[cls]['W'], wl_W, K_w_sqr[wl_W])\n",
" dwr = 10 * np.log10(zeX / zeW)\n",
" dwr_curves_1p[cls] = dwr\n",
" ax.plot(D_grid, dwr, color=COLORS[cls], lw=1.8, label=cls)\n",
"ax.axhline(0.0, color='k', ls=':', alpha=0.5, label='Rayleigh (DWR=0)')\n",
"ax.set_xlabel('diameter D [mm]')\n",
"ax.set_ylabel('single-particle DWR$_{X-W}$ [dB]')\n",
"ax.set_title('Where does each habit walk out of Rayleigh?')\n",
"ax.set_xlim(0, 8); ax.grid(True, alpha=0.3); ax.legend(fontsize=9)\n",
"fig.tight_layout();\n",
"\n",
"print('Approximate D at which single-particle DWR crosses +3 dB:')\n",
"for cls, dwr in dwr_curves_1p.items():\n",
" above = np.where(dwr > 3.0)[0]\n",
" label = f'{D_grid[above[0]]:.2f} mm' if len(above) else '> 8 mm (stays Rayleigh)'\n",
" print(f' {cls:<12} {label}')\n",
),
md("## Bulk DWR vs. median-volume diameter D₀\n",
"\n",
"The single-particle curves fold through a PSD in the bulk radar\n",
"measurement. Hold the habit fixed, sweep the exponential PSD slope Λ\n",
"(so D₀ = 3.67/Λ varies), and plot the resulting bulk DWR. The curve\n",
"tells you how sensitive each habit's dual-frequency signature is to\n",
"the *population-average* particle size.\n"),
code(
"def bulk_Zh(wl, cls, N0, lam, Dmax):\n",
" s = Scatterer(wavelength=wl, m=refractive_index(wl, cls),\n",
" Kw_sqr=K_w_sqr[wl], ddelt=1e-4, ndgs=2)\n",
" integ = PSDIntegrator()\n",
" integ.D_max = Dmax; integ.num_points = 128\n",
" integ.geometries = (geom_vert_back,)\n",
" if CLASSES[cls]['axis_ratio'] is None:\n",
" integ.axis_ratio_func = lambda D: 1.0 / dsr_thurai_2007(D)\n",
" else:\n",
" s.axis_ratio = CLASSES[cls]['axis_ratio']\n",
" s.psd_integrator = integ\n",
" s.psd_integrator.init_scatter_table(s)\n",
" s.psd = ExponentialPSD(N0=N0, Lambda=lam, D_max=Dmax)\n",
" s.set_geometry(geom_vert_back)\n",
" return radar.refl(s)\n",
"\n",
"# Hold N0 fixed, sweep Λ so D0 = 3.67/Λ spans 0.4–3 mm\n",
"Lambdas = np.linspace(1.2, 9.0, 10) # mm⁻¹\n",
"D0s = 3.67 / Lambdas\n",
"DMAX = 8.0; N0 = 8e3\n",
"\n",
"dwr_curves = {}\n",
"for cls in CLASSES:\n",
" dwr = []\n",
" for lam in Lambdas:\n",
" zX = bulk_Zh(wl_X, cls, N0, lam, DMAX)\n",
" zW = bulk_Zh(wl_W, cls, N0, lam, DMAX)\n",
" dwr.append(10 * np.log10(zX / zW))\n",
" dwr_curves[cls] = np.asarray(dwr)\n",
"\n",
"fig, ax = plt.subplots(figsize=(8, 4.5))\n",
"for cls, d in dwr_curves.items():\n",
" ax.plot(D0s, d, color=COLORS[cls], lw=1.8, marker='o',\n",
" markersize=4, label=cls)\n",
"ax.axhline(0.0, color='k', ls=':', alpha=0.5)\n",
"ax.set_xlabel('median-volume diameter $D_0$ [mm]')\n",
"ax.set_ylabel('bulk DWR$_{X-W}$ [dB]')\n",
"ax.set_title('Bulk DWR vs. PSD size — one line per hydrometeor class')\n",
"ax.grid(True, alpha=0.3); ax.legend(fontsize=9)\n",
"fig.tight_layout();\n",
),
md("## Dual-frequency Doppler spectra — aggregate vs. graupel\n",
"\n",
"Now pull the spectral version together. Two bulk populations that\n",
"*can't* be told apart by reflectivity alone — a low-density aggregate\n",
"PSD and a graupel PSD tuned to give the same X-band Z_h — show very\n",
"different sDWR(v) profiles. Aggregates have a mostly-monotonic rise\n",
"because their σ_b(D) curve stays smooth; graupel's spectrum carries\n",
"the Mie-resonance structure of its single-particle σ_b(D) through into\n",
"velocity space, since every drop size maps deterministically to its\n",
"own fall speed.\n"),
code(
"def build(cls, N0, lam, Dmax, wl):\n",
" s = Scatterer(wavelength=wl, m=refractive_index(wl, cls),\n",
" Kw_sqr=K_w_sqr[wl], ddelt=1e-4, ndgs=2)\n",
" integ = PSDIntegrator()\n",
" integ.D_max = Dmax; integ.num_points = 256\n",
" integ.geometries = (geom_vert_back,)\n",
" if CLASSES[cls]['axis_ratio'] is None:\n",
" integ.axis_ratio_func = lambda D: 1.0 / dsr_thurai_2007(D)\n",
" else:\n",
" s.axis_ratio = CLASSES[cls]['axis_ratio']\n",
" s.psd_integrator = integ\n",
" s.psd_integrator.init_scatter_table(s)\n",
" s.psd = ExponentialPSD(N0=N0, Lambda=lam, D_max=Dmax)\n",
" return s\n",
"\n",
"# PSDs chosen so both give ~moderate snowfall X-band reflectivity.\n",
"# Aggregate: low density, broader PSD to Dmax=6 mm.\n",
"# Graupel: medium density, narrower PSD (smaller max D).\n",
"agg_X = build('low-ρ agg', N0=2e4, lam=2.0, Dmax=6.0, wl=wl_X)\n",
"agg_W = build('low-ρ agg', N0=2e4, lam=2.0, Dmax=6.0, wl=wl_W)\n",
"gra_X = build('graupel', N0=5e3, lam=3.0, Dmax=4.0, wl=wl_X)\n",
"gra_W = build('graupel', N0=5e3, lam=3.0, Dmax=4.0, wl=wl_W)\n",
"\n",
"fall_agg = spectra.fall_speed.locatelli_hobbs_1974_aggregates\n",
"fall_gra = spectra.fall_speed.locatelli_hobbs_1974_graupel_hex\n",
"turb = spectra.GaussianTurbulence(0.2)\n",
"V_MIN, V_MAX = -0.5, 4.0\n",
"\n",
"def run(sc, fall):\n",
" return SpectralIntegrator(\n",
" sc, fall_speed=fall, turbulence=turb,\n",
" v_min=V_MIN, v_max=V_MAX, n_bins=1024,\n",
" geometry_backscatter=geom_vert_back,\n",
" ).run()\n",
"\n",
"r_agg_X = run(agg_X, fall_agg); r_agg_W = run(agg_W, fall_agg)\n",
"r_gra_X = run(gra_X, fall_gra); r_gra_W = run(gra_W, fall_gra)\n",
"\n",
"def dBZ(x):\n",
" return 10 * np.log10(np.maximum(x, 1e-12))\n",
"\n",
"print(f'{\"population\":<15} {\"Z_h X [dBZ]\":>12} {\"Z_h W [dBZ]\":>12} '\n",
" f'{\"DWR [dB]\":>10}')\n",
"for name, rX, rW in [('low-ρ agg', r_agg_X, r_agg_W),\n",
" ('graupel', r_gra_X, r_gra_W)]:\n",
" zX = 10 * np.log10(np.trapezoid(rX.sZ_h, rX.v))\n",
" zW = 10 * np.log10(np.trapezoid(rW.sZ_h, rW.v))\n",
" print(f'{name:<15} {zX:>12.2f} {zW:>12.2f} {zX - zW:>10.2f}')\n",
"\n",
"fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(8, 6.5), sharex=True)\n",
"ax1.plot(r_agg_X.v, dBZ(r_agg_X.sZ_h), 'C2-', label='agg, X-band')\n",
"ax1.plot(r_agg_W.v, dBZ(r_agg_W.sZ_h), 'C2--', label='agg, W-band')\n",
"ax1.plot(r_gra_X.v, dBZ(r_gra_X.sZ_h), 'C1-', label='graupel, X-band')\n",
"ax1.plot(r_gra_W.v, dBZ(r_gra_W.sZ_h), 'C1--', label='graupel, W-band')\n",
"ax1.set_ylabel('sZ_h [dBZ / (m/s)]')\n",
"ax1.set_title('Dual-frequency Doppler spectra')\n",
"ax1.legend(fontsize=9); ax1.grid(True, alpha=0.3)\n",
"\n",
"sdwr_agg = dBZ(r_agg_X.sZ_h) - dBZ(r_agg_W.sZ_h)\n",
"sdwr_gra = dBZ(r_gra_X.sZ_h) - dBZ(r_gra_W.sZ_h)\n",
"mask_agg = (r_agg_X.sZ_h > 1e-8) & (r_agg_W.sZ_h > 1e-10)\n",
"mask_gra = (r_gra_X.sZ_h > 1e-8) & (r_gra_W.sZ_h > 1e-10)\n",
"ax2.plot(r_agg_X.v[mask_agg], sdwr_agg[mask_agg], 'C2-',\n",
" lw=1.6, label='low-ρ agg')\n",
"ax2.plot(r_gra_X.v[mask_gra], sdwr_gra[mask_gra], 'C1-',\n",
" lw=1.6, label='graupel')\n",
"ax2.axhline(0.0, color='k', ls=':', alpha=0.5)\n",
"ax2.set_xlabel('Doppler velocity v [m/s]')\n",
"ax2.set_ylabel('sDWR$_{X-W}$ [dB]')\n",
"ax2.set_title('Spectral DWR — aggregates smooth, graupel resonant')\n",
"ax2.set_xlim(V_MIN, V_MAX)\n",
"ax2.legend(fontsize=9); ax2.grid(True, alpha=0.3)\n",
"fig.tight_layout();\n",
),
md("## Takeaway\n",
"\n",
"Following Honeyager 2013's framing, we've used a single T-matrix\n",
"spheroid code path with four choices of (ρ_eff, axis ratio) to span\n",
"four physically distinct hydrometeor classes. Two observations pop\n",
"out that matter for interpreting dual-frequency radar data:\n",
"\n",
"* **The *size* at which DWR departs from zero is set by χ = πD/λ**\n",
" (≈ 1 mm at W-band), not by ρ_eff. But **the amplitude and shape**\n",
" of the DWR(D) curve beyond that threshold is a direct fingerprint\n",
" of ρ_eff: dense classes oscillate strongly, low-ρ aggregates rise\n",
" smoothly and slowly.\n",
"* **Bulk DWR sweeps out very different trajectories with D₀**, so a\n",
" single (Z_h, DWR) measurement can be mapped to a habit-dependent\n",
" D₀ estimate provided you commit to a class — this is exactly the\n",
" lookup-table strategy that papers like Mason et al. 2019 and\n",
" Tridon et al. 2019 build their retrievals on.\n",
"* **Spectral DWR resolves habit ambiguity** that bulk DWR can't:\n",
" graupel's Mie resonances appear as wiggles in sDWR(v), while\n",
" low-ρ aggregates produce a smooth monotonic rise.\n",
"\n",
"Swapping (ρ_eff, axis_ratio) for any other habit class (dendrites,\n",
"rimed aggregates, hail) is one `CLASSES` dict entry away.\n"),
]
NB12 = [
md(
"# Tutorial 12 — Spectral polarimetry of a SLW / rain / hail mixture at 500 hPa\n",
"\n",
"**Reference.** Lakshmi K., K., Sahoo, S., Biswas, S. K., and\n",
"Chandrasekar, V., 2024: Study of Microphysical Signatures Based on\n",
"Spectral Polarimetry during the RELAMPAGO Field Experiment in\n",
"Argentina. *J. Atmos. Oceanic Technol.*, 41, 235–256,\n",
"[doi:10.1175/JTECH-D-22-0113.1](https://doi.org/10.1175/JTECH-D-22-0113.1).\n",
"\n",
"Lakshmi et al. (2024) used C-band Doppler spectral polarimetry from\n",
"the CSU-CHIVO radar during RELAMPAGO to dissect mixed-phase and\n",
"convective precipitation volumes. Their central point is that\n",
"*spectral* $Z_h$, $Z_\\mathrm{dr}$, $K_\\mathrm{dp}$, $\\rho_\\mathrm{hv}$\n",
"separate hydrometeor populations that the *bulk* observables mash\n",
"into one number — because each species occupies a distinct\n",
"Doppler-velocity window set by its terminal fall speed.\n",
"\n",
"Here we build a synthetic C-band radar resolution volume at altitude\n",
"where $P = 500$ hPa ($T \\approx 252$ K, $\\rho \\approx 0.69$ kg/m³,\n",
"reached near 6 km MSL in a convective updraft). Three populations\n",
"coexist:\n",
"\n",
"* **Supercooled cloud liquid water** (SLW) — spherical droplets,\n",
" $D \\lesssim 0.2$ mm, $v_t \\lesssim 0.1$ m/s.\n",
"* **Rain** — oblate drops following the Thurai (2007) shape,\n",
" $D \\lesssim 5$ mm.\n",
"* **Small wet (melting) hail** — water-coated ice, 30 %% meltwater\n",
" by volume via Maxwell-Garnett EMA, axis ratio 0.75, canting σ = 40°,\n",
" $D$ up to 12 mm. The water coating drives a strong C-band Mie\n",
" resonance near $D \\approx 8$–10 mm with large $|\\delta_\\mathrm{hv}|$\n",
" — the rain-hail mixing regime flagged by Lakshmi et al. (2024) at the\n",
" melting layer.\n",
"\n",
"## Slant geometry\n",
"Polarimetric observables require a slant beam: oblate raindrops at\n",
"nadir project to circles, collapsing $Z_\\mathrm{dr}$ and\n",
"$K_\\mathrm{dp}$ to zero. We use an elevation angle $\\phi = 30°$,\n",
"so the radial velocity of a particle of diameter $D$ is\n",
"$v_r(D) = v_t(D)\\,\\sin\\phi$ and the scattering matrices come from\n",
"the slant back-/forward-scatter geometries $(\\theta_0, \\theta) =\n",
"(60°, 120°)$ and $(60°, 60°)$. Lakshmi et al. (2024) work with\n",
"CSU-CHIVO RHI scans at elevation angles that span $0°$–$45°$; our\n",
"$\\phi = 30°$ sits in their mid-elevation band where horizontal-wind\n",
"and fall-speed contributions to the radial velocity are both\n",
"meaningful.\n"),
code(
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"\n",
"from rustmatrix import (HydroMix, MixtureComponent, Scatterer,\n",
" SpectralIntegrator, orientation, radar, spectra)\n",
"from rustmatrix.psd import ExponentialPSD, GammaPSD, PSDIntegrator\n",
"from rustmatrix.refractive import m_w_0C, mg_refractive, mi\n",
"from rustmatrix.tmatrix_aux import K_w_sqr, dsr_thurai_2007, wl_C\n",
"\n",
"# Ambient state at 500 hPa.\n",
"P_HPA = 500.0\n",
"T_K = 252.0\n",
"R_D = 287.05\n",
"RHO_AIR = P_HPA * 100.0 / (R_D * T_K)\n",
"RHO_0 = 1.225\n",
"RHO_RATIO = RHO_0 / RHO_AIR\n",
"DENS_CORR_POW4 = RHO_RATIO ** 0.4 # Foote-duToit for rain\n",
"DENS_CORR_SQRT = RHO_RATIO ** 0.5 # hail / large particles\n",
"\n",
"# 30° slant geometry.\n",
"PHI_DEG = 30.0\n",
"SIN_PHI = np.sin(np.deg2rad(PHI_DEG))\n",
"THETA0 = 90.0 - PHI_DEG\n",
"GEOM_BACK = (THETA0, 180.0 - THETA0, 0.0, 180.0, 0.0, 0.0)\n",
"GEOM_FORW = (THETA0, THETA0, 0.0, 0.0, 0.0, 0.0)\n",
"\n",
"V_MIN, V_MAX, N_BINS = -2.0, 15.0, 512\n",
"\n",
"print(f'ρ_air = {RHO_AIR:.3f} kg/m³, ρ₀/ρ = {RHO_RATIO:.3f}')\n",
"print(f'rain fall-speed ×{DENS_CORR_POW4:.3f}, hail ×{DENS_CORR_SQRT:.3f}')\n",
"print(f'φ = {PHI_DEG:.0f}°, sin φ = {SIN_PHI:.3f}')\n",
),
md("## Build the three hydrometeor scatterers\n",
"\n",
"Each species uses the slant geometry tables for its scattering matrix.\n"),
code(
"def build_rain():\n",
" s = Scatterer(wavelength=wl_C, m=m_w_0C[wl_C], Kw_sqr=K_w_sqr[wl_C],\n",
" ddelt=1e-4, ndgs=2)\n",
" integ = PSDIntegrator()\n",
" integ.D_max = 5.0\n",
" integ.num_points = 64\n",
" integ.axis_ratio_func = lambda D: 1.0 / dsr_thurai_2007(D)\n",
" integ.geometries = (GEOM_BACK, GEOM_FORW)\n",
" s.psd_integrator = integ\n",
" s.psd_integrator.init_scatter_table(s)\n",
" # Heavy convective rain: D0 = 2.0 mm, Nw = 8e4 ⇒ Z_h ≈ 52 dBZ.\n",
" s.psd = GammaPSD(D0=2.0, Nw=8e4, mu=2, D_max=5.0)\n",
" return s\n",
"\n",
"def build_slw():\n",
" s = Scatterer(wavelength=wl_C, m=m_w_0C[wl_C], Kw_sqr=K_w_sqr[wl_C],\n",
" axis_ratio=1.0, ddelt=1e-4, ndgs=2)\n",
" integ = PSDIntegrator()\n",
" integ.D_max = 0.2\n",
" integ.num_points = 64\n",
" integ.geometries = (GEOM_BACK, GEOM_FORW)\n",
" s.psd_integrator = integ\n",
" s.psd_integrator.init_scatter_table(s)\n",
" # D₀ ≈ 30 µm, LWC ≈ 0.5 g/m³.\n",
" s.psd = GammaPSD(D0=0.03, Nw=1e11, mu=4, D_max=0.2)\n",
" return s\n",
"\n",
"# Wet-hail effective refractive index: Maxwell-Garnett with water as\n",
"# the matrix (30% meltwater by volume) and ice as the inclusion.\n",
"m_wet_hail = mg_refractive((m_w_0C[wl_C], mi(wl_C, 0.917)),\n",
" (0.30, 0.70))\n",
"\n",
"def build_hail():\n",
" s = Scatterer(wavelength=wl_C, m=m_wet_hail, Kw_sqr=K_w_sqr[wl_C],\n",
" axis_ratio=0.75, ddelt=1e-4, ndgs=2)\n",
" s.orient = orientation.orient_averaged_fixed\n",
" s.or_pdf = orientation.gaussian_pdf(std=40.0, mean=90.0)\n",
" s.n_alpha = 6; s.n_beta = 12\n",
" integ = PSDIntegrator()\n",
" integ.D_max = 12.0\n",
" integ.num_points = 64\n",
" integ.geometries = (GEOM_BACK, GEOM_FORW)\n",
" s.psd_integrator = integ\n",
" s.psd_integrator.init_scatter_table(s)\n",
" # Small-wet-hail PSD: N0 = 150 m⁻³ mm⁻¹, Λ = 0.6 mm⁻¹, D_max = 12 mm.\n",
" s.psd = ExponentialPSD(N0=150.0, Lambda=0.6, D_max=12.0)\n",
" return s\n",
"\n",
"rain = build_rain()\n",
"slw = build_slw()\n",
"hail = build_hail()\n",
),
md("## Fall-speed callables\n",
"\n",
"Each callable returns the radial velocity (already projected onto the\n",
"slant beam by $\\sin\\varepsilon$) and includes the $(\\rho_0/\\rho)$\n",
"correction appropriate for its size regime.\n"),
code(
"def v_rain(D):\n",
" # Beard 1976 handles (T, P) density correction itself.\n",
" v_t = spectra.fall_speed.beard_1976(D, T=T_K, P=P_HPA * 100.0)\n",
" return v_t * SIN_PHI\n",
"\n",
"_hail_power = spectra.fall_speed.power_law(a=9.0, b=0.64, D_ref=10.0)\n",
"\n",
"def v_hail(D):\n",
" # Matson-Huggins v ≈ 9 (D/1cm)^0.64, density-corrected.\n",
" return _hail_power(D) * DENS_CORR_SQRT * SIN_PHI\n",
"\n",
"def v_slw(D):\n",
" # Stokes-regime drag scaled by (ρ₀/ρ)^0.4.\n",
" return 3.0 * np.asarray(D, dtype=float) ** 2 * DENS_CORR_POW4 * SIN_PHI\n",
"\n",
"D_probe = np.array([0.03, 0.1, 1.0, 3.0, 5.0, 10.0, 20.0])\n",
"print(f'{\"D [mm]\":>8} {\"v_slw\":>8} {\"v_rain\":>8} {\"v_hail\":>8} [m/s]')\n",
"for D in D_probe:\n",
" vs, vr, vh = np.atleast_1d(v_slw(D))[0], np.atleast_1d(v_rain(D))[0], np.atleast_1d(v_hail(D))[0]\n",
" print(f'{D:>8.3f} {vs:>8.3f} {vr:>8.3f} {vh:>8.3f}')\n",
),
md("## Bulk single-species observables\n",
"\n",
"Before slicing spectra, print the bulk $Z_h$, $Z_\\mathrm{dr}$,\n",
"$\\rho_\\mathrm{hv}$, $K_\\mathrm{dp}$ for each species alone. Note\n",
"three things in the printout below:\n",
"\n",
"* **Rain dominates $K_\\mathrm{dp}$** ($\\sim$17 °/km vs $\\sim$1.5 for\n",
" hail). Wet hail has a broad canting distribution (σ = 40°), which\n",
" smears its forward-scatter differential phase toward zero, while\n",
" Thurai-shaped raindrops aligned with gravity produce strong,\n",
" coherent $K_\\mathrm{dp}$.\n",
"* **Hail's reflectivity is comparable to rain's** at these\n",
" concentrations (62.9 vs 57.8 dBZ) — the C-band resonance on the\n",
" wet-hail PSD tail is loud.\n",
"* **Hail's $\\rho_\\mathrm{hv}$ sits well below unity** (≈ 0.86)\n",
" because the resonant oscillations in $f_h - f_v$ span a wide\n",
" diameter range.\n"),
code(
"def bulk(sc):\n",
" sc.set_geometry(GEOM_BACK)\n",
" Z = 10 * np.log10(radar.refl(sc))\n",
" Zdr = 10 * np.log10(radar.Zdr(sc))\n",
" rho = radar.rho_hv(sc)\n",
" sc.set_geometry(GEOM_FORW)\n",
" Kdp = radar.Kdp(sc)\n",
" return Z, Zdr, rho, Kdp\n",
"\n",
"print(f\" {'species':<8} {'Z_h':>8} {'Z_dr':>7} {'ρ_hv':>8} {'K_dp':>9}\")\n",
"print(f\" {'':<8} {'[dBZ]':>8} {'[dB]':>7} {'':>8} {'[°/km]':>9}\")\n",
"for name, sc in (('SLW', slw), ('rain', rain), ('hail', hail)):\n",
" Z, Zdr, rho, Kdp = bulk(sc)\n",
" print(f' {name:<8} {Z:>8.2f} {Zdr:>+7.3f} {rho:>8.5f} {Kdp:>+9.4f}')\n",
),
md("## Spectral integrators\n",
"\n",
"Run one `SpectralIntegrator` per species (so we can plot the\n",
"contribution of each component) plus one on a `HydroMix` with\n",
"per-component kinematics (the combined spectrum).\n"),
code(
"turb = spectra.GaussianTurbulence(0.5)\n",
"\n",
"def run_single(sc, fall):\n",
" return SpectralIntegrator(\n",
" sc, fall_speed=fall, turbulence=turb,\n",
" v_min=V_MIN, v_max=V_MAX, n_bins=N_BINS,\n",
" geometry_backscatter=GEOM_BACK,\n",
" geometry_forward=GEOM_FORW,\n",
" ).run()\n",
"\n",
"r_slw = run_single(slw, v_slw)\n",
"r_rain = run_single(rain, v_rain)\n",
"r_hail = run_single(hail, v_hail)\n",
"\n",
"mix = HydroMix([\n",
" MixtureComponent(slw, slw.psd, 'slw'),\n",
" MixtureComponent(rain, rain.psd, 'rain'),\n",
" MixtureComponent(hail, hail.psd, 'hail'),\n",
"])\n",
"r_mix = SpectralIntegrator(\n",
" mix, component_kinematics={\n",
" 'slw': (v_slw, turb),\n",
" 'rain': (v_rain, turb),\n",
" 'hail': (v_hail, turb),\n",
" },\n",
" v_min=V_MIN, v_max=V_MAX, n_bins=N_BINS,\n",
" geometry_backscatter=GEOM_BACK,\n",
" geometry_forward=GEOM_FORW,\n",
").run()\n",
"\n",
"# Second scenario: hail concentration halved (N0 = 75) so we can\n",
"# compare the spectrum of each observable against the baseline mix.\n",
"hail_psd_half = ExponentialPSD(N0=75.0, Lambda=0.6, D_max=12.0)\n",
"mix_half = HydroMix([\n",
" MixtureComponent(slw, slw.psd, 'slw'),\n",
" MixtureComponent(rain, rain.psd, 'rain'),\n",
" MixtureComponent(hail, hail_psd_half, 'hail'),\n",
"])\n",
"r_mix_half = SpectralIntegrator(\n",
" mix_half, component_kinematics={\n",
" 'slw': (v_slw, turb),\n",
" 'rain': (v_rain, turb),\n",
" 'hail': (v_hail, turb),\n",
" },\n",
" v_min=V_MIN, v_max=V_MAX, n_bins=N_BINS,\n",
" geometry_backscatter=GEOM_BACK,\n",
" geometry_forward=GEOM_FORW,\n",
").run()\n",
"\n",
"v = r_mix.v\n",
"print(f'v-grid: {v[0]:+.2f} … {v[-1]:+.2f} m/s, N={len(v)}')\n",
),
md("## Spectral reflectivity $sZ_h(v)$\n",
"\n",
"Each species occupies its own Doppler-velocity window. SLW sits at\n",
"$v \\approx 0$, rain peaks near 3–5 m/s (with its large-drop tail\n",
"extending a bit further), and hail's spectral peak is near $v \\approx\n",
"5$ m/s — its exponential PSD front-loads the small-diameter (slow)\n",
"end, so most hail mass moves at modest velocities, with the\n",
"large-hail tail extending out to $v \\approx 7$ m/s. Hail *dominates*\n",
"the spectrum past $v \\approx 7$ m/s simply because rain has run out of\n",
"drops by then ($D_\\mathrm{max}^\\mathrm{rain} = 5$ mm). The mixture\n",
"spectrum is the incoherent sum — the loudest species wins each bin.\n",
"\n",
"**Compared with Lakshmi et al. (2024).** Their Fig. 8 (the 14 Dec\n",
"2018 convective case, altitudes 1.5–2.5 km below the melting layer)\n",
"shows broad, sometimes bimodal $sZ_h$ spectra in the 22–30 dB range\n",
"whenever rain and partially melted hail share the volume. Their\n",
"Fig. 2 reports bulk $Z_h \\ge 35$–40 dBZ in the rain layer below the\n",
"melting band. Our synthetic bulk $Z_h$ in the mix is $\\sim$63 dBZ\n",
"(very heavy convection, brighter than either of their case studies)\n",
"but the *shape* of the spectrum — rain peak near 3–5 m/s, broad\n",
"continuation into a hail-dominated fast tail — reproduces the\n",
"rain+hail bimodality they report at 5 km altitude in their 240°\n",
"RHI scan (paper text, p. 242).\n"),
code(
"def dB(x):\n",
" return 10 * np.log10(np.maximum(x, 1e-12))\n",
"\n",
"fig, ax = plt.subplots(figsize=(8, 4))\n",
"ax.plot(v, dB(r_slw.sZ_h), color='tab:blue', lw=1.2, label='SLW')\n",
"ax.plot(v, dB(r_rain.sZ_h), color='tab:green', lw=1.2, label='rain')\n",
"ax.plot(v, dB(r_hail.sZ_h), color='tab:red', lw=1.2, label='hail')\n",
"ax.plot(v, dB(r_mix.sZ_h), color='black', lw=1.8,\n",
" label='mixture', linestyle='--')\n",
"ax.plot(v, dB(r_mix_half.sZ_h), color='dimgray', lw=1.8,\n",
" label='mixture (½ hail)', linestyle=':')\n",
"ax.set_xlabel('Doppler velocity [m/s]')\n",
"ax.set_ylabel(r'$sZ_h$ [dBZ / (m s$^{-1}$)]')\n",
"ax.set_title(r'Spectral reflectivity at C-band, $\\phi$ = 30°, 500 hPa')\n",
"ax.set_xlim(V_MIN, V_MAX)\n",
"ax.set_ylim(-40, 80)\n",
"ax.grid(True, alpha=0.3)\n",
"ax.legend(loc='upper right')\n",
"plt.tight_layout()\n",
"plt.show()\n",
),
md("## Spectral differential reflectivity $sZ_\\mathrm{dr}(v)$\n",
"\n",
"Rain's $sZ_\\mathrm{dr}$ climbs from 0 dB at the small-drop end to\n",
"$\\sim$+2.5 dB at the fast end — a clean size-sorting signal, since\n",
"larger drops fall faster and are more oblate.\n",
"\n",
"Wet hail shows the hallmark **C-band Mie resonance notch**: slightly\n",
"positive at small $v$ (small D, Rayleigh), then diving to ≈ −1 dB as\n",
"the $D \\sim 8$–10 mm wet-hail population hits the C-band resonance\n",
"(the water coating pushes the resonance to smaller $D$ than dry hail;\n",
"see Kumjian, Ryzhkov et al. for the phenomenology).\n",
"\n",
"The **mixture** curve is the interesting one. Around $v \\approx 5$ m/s\n",
"rain's positive $sZ_\\mathrm{dr}$ (≈ +1 dB) and hail's negative swing\n",
"(≈ −0.4 dB) carry comparable power, and the mixture $sZ_\\mathrm{dr}$\n",
"collapses toward zero — a velocity bin where neither species wins.\n",
"\n",
"**Compared with Lakshmi et al. (2024).** Their Fig. 8 shows\n",
"$sZ_\\mathrm{dr}$ climbing monotonically with radial velocity in the\n",
"liquid-phase rain layer — slopes of +0.455 and +0.57 dB (m s⁻¹)⁻¹ at\n",
"2.5- and 2-km altitudes, interpreted explicitly as shear-induced\n",
"size sorting where larger/more-oblate drops sit at higher $v$ (paper\n",
"pp. 242–243). Our rain-alone curve rises from $\\sim$0 to +2.5 dB\n",
"across a $\\sim$12 m/s window — the same qualitative signature,\n",
"offset by our 30° elevation and without the paper's strong shear\n",
"advecting small drops. Lakshmi et al. also report $Z_\\mathrm{dr}$\n",
"values of 6–8 dB in the convective core at 2 km, attributed\n",
"explicitly to *partially melted hail* depolarising horizontally —\n",
"that is the regime our wet-hail EMA is modelling, although our\n",
"narrower wet-hail PSD produces a spectral *notch* rather than a\n",
"broad +6 dB plateau.\n"),
code(
"fig, ax = plt.subplots(figsize=(8, 4))\n",
"ax.plot(v, dB(r_slw.sZ_dr), color='tab:blue', lw=1.2, label='SLW')\n",
"ax.plot(v, dB(r_rain.sZ_dr), color='tab:green', lw=1.2, label='rain')\n",
"ax.plot(v, dB(r_hail.sZ_dr), color='tab:red', lw=1.2, label='hail')\n",
"ax.plot(v, dB(r_mix.sZ_dr), color='black', lw=1.8,\n",
" label='mixture', linestyle='--')\n",
"ax.plot(v, dB(r_mix_half.sZ_dr), color='dimgray', lw=1.8,\n",
" label='mixture (½ hail)', linestyle=':')\n",
"ax.set_xlabel('Doppler velocity [m/s]')\n",
"ax.set_ylabel(r'$sZ_\\mathrm{dr}$ [dB]')\n",
"ax.set_title(r'Spectral differential reflectivity')\n",
"ax.set_xlim(V_MIN, V_MAX)\n",
"ax.set_ylim(-10, 4)\n",
"ax.axhline(0.0, color='gray', lw=0.6)\n",
"ax.grid(True, alpha=0.3)\n",
"ax.legend(loc='lower left')\n",
"plt.tight_layout()\n",
"plt.show()\n",
),
md("## Spectral specific differential phase $sK_\\mathrm{dp}(v)$\n",
"\n",
"The spectral proxy for $\\phi_\\mathrm{dp}$ is $sK_\\mathrm{dp}$: the\n",
"forward-scatter differential phase contributed by particles in each\n",
"velocity bin. In this mixture it is **rain-dominated**. The bulk\n",
"$K_\\mathrm{dp}$ was $\\sim$17 °/km for rain versus $\\sim$1.5 °/km for\n",
"wet hail — Thurai-shaped, gravity-aligned raindrops are nearly ideal\n",
"oblate forward-scatterers, while wet hail's broad canting\n",
"distribution (σ = 40°) smears its differential phase toward zero.\n",
"\n",
"The rain peak in $sK_\\mathrm{dp}(v)$ lands at $v \\approx 5$ m/s —\n",
"that is where the large oblate raindrops ($D \\sim 3$–5 mm, the ones\n",
"that carry almost all the differential phase shift) sit in Doppler\n",
"space. Hail contributes only a small bump at similar velocities and\n",
"is essentially invisible against rain in the mixture curve. Halving\n",
"hail's concentration (gray dotted) barely moves $sK_\\mathrm{dp}$ at\n",
"all.\n",
"\n",
"**Compared with Lakshmi et al. (2024).** The paper does not show\n",
"$sK_\\mathrm{dp}$ spectra directly — their spectral analysis is\n",
"restricted to $sZ_h$, $sZ_\\mathrm{dr}$, and $s\\rho_\\mathrm{hv}$\n",
"(Eqs. 3–5, p. 240). The physics our curves illustrate — differential\n",
"phase is *rain-dominated* whenever rain coexists with tumbling\n",
"ice-phase scatterers — is why Lakshmi et al. rely on\n",
"$sZ_\\mathrm{dr}$ and $s\\rho_\\mathrm{hv}$ rather than $K_\\mathrm{dp}$\n",
"to fingerprint the ice fraction of a mixed-phase volume.\n"),
code(
"fig, ax = plt.subplots(figsize=(8, 4))\n",
"ax.plot(v, r_slw.sKdp, color='tab:blue', lw=1.2, label='SLW')\n",
"ax.plot(v, r_rain.sKdp, color='tab:green', lw=1.2, label='rain')\n",
"ax.plot(v, r_hail.sKdp, color='tab:red', lw=1.2, label='hail')\n",
"ax.plot(v, r_mix.sKdp, color='black', lw=1.8,\n",
" label='mixture', linestyle='--')\n",
"ax.plot(v, r_mix_half.sKdp, color='dimgray', lw=1.8,\n",
" label='mixture (½ hail)', linestyle=':')\n",
"ax.set_xlabel('Doppler velocity [m/s]')\n",
"ax.set_ylabel(r'$sK_\\mathrm{dp}$ [° / km / (m s$^{-1}$)]')\n",
"ax.set_title(r'Spectral specific differential phase')\n",
"ax.set_xlim(V_MIN, V_MAX)\n",
"ax.axhline(0.0, color='gray', lw=0.6)\n",
"ax.grid(True, alpha=0.3)\n",
"ax.legend(loc='upper right')\n",
"plt.tight_layout()\n",
"plt.show()\n",
),
md("## Spectral copolar correlation coefficient $s\\rho_\\mathrm{hv}(v)$\n",
"\n",
"Per-species $s\\rho_\\mathrm{hv}$ is close to 1 wherever the species has\n",
"a single, narrowly distributed polarimetric response. Wet hail drops\n",
"to ≈ 0.81 at its fast-velocity end — the Mie-resonant mix of\n",
"water-coated oblate hailstones spans enough backscatter-phase spread\n",
"that the H–V correlation falls sharply.\n",
"\n",
"The **mixture** curve tracks hail wherever hail dominates. Where rain\n",
"contributes comparable power (v ≈ 3–5 m/s), the mixture\n",
"$s\\rho_\\mathrm{hv}$ sits *between* rain (≈ 1) and hail (≈ 0.89) —\n",
"rain is effectively diluting hail's phase spread with its own\n",
"well-correlated H–V returns. Dropping the mixture curve *below* both\n",
"components would require rain and hail to have opposite-signed\n",
"$\\delta_\\mathrm{hv}$ at the same velocity, which does not quite\n",
"happen here; the melting-layer classical $\\rho_\\mathrm{hv}$ dip\n",
"(0.85–0.9) needs that extra ingredient.\n",
"\n",
"Halving the hail concentration (gray dotted curve) drags the\n",
"mixture $s\\rho_\\mathrm{hv}$ *closer to unity* in the overlap region\n",
"— with less hail phase spread contaminating the volume, the\n",
"well-correlated rain returns dominate.\n",
"\n",
"**Compared with Lakshmi et al. (2024).** The paper reports\n",
"$s\\rho_\\mathrm{hv}$ dropping to $\\sim$0.84–0.99 near the melting\n",
"layer in the 30 Nov 2018 stratiform case (Fig. 13, altitudes 3–5 km\n",
"at 33 km range), with the lowest values coinciding with the\n",
"rain+hail mixture class in their DROPS2 hydrometeor classification\n",
"(Fig. 11). Our wet-hail-alone curve bottoms at $\\sim$0.81 at the\n",
"resonance tail and the mixture curve is pulled to $\\sim$0.92 in the\n",
"rain-hail overlap region — a quantitative match to the mid- and\n",
"low-$s\\rho_\\mathrm{hv}$ signatures Lakshmi et al. use to flag\n",
"mixed-phase volumes. More broadly, their Fig. 2 shows bulk\n",
"$\\rho_\\mathrm{hv} \\approx 1$ in pure rain below the melting layer\n",
"and a sharp drop crossing it — exactly the rain-to-mixture drop our\n",
"spectrum traces as $v$ climbs from rain-dominated (≈ 1) into the\n",
"hail-resonance tail (≈ 0.82).\n"),
code(
"fig, ax = plt.subplots(figsize=(8, 4))\n",
"ax.plot(v, r_slw.srho_hv, color='tab:blue', lw=1.2, label='SLW')\n",
"ax.plot(v, r_rain.srho_hv, color='tab:green', lw=1.2, label='rain')\n",
"ax.plot(v, r_hail.srho_hv, color='tab:red', lw=1.2, label='hail')\n",
"ax.plot(v, r_mix.srho_hv, color='black', lw=1.8,\n",
" label='mixture', linestyle='--')\n",
"ax.plot(v, r_mix_half.srho_hv, color='dimgray', lw=1.8,\n",
" label='mixture (½ hail)', linestyle=':')\n",
"ax.set_xlabel('Doppler velocity [m/s]')\n",
"ax.set_ylabel(r'$s\\rho_\\mathrm{hv}$')\n",
"ax.set_title(r'Spectral copolar correlation coefficient')\n",
"ax.set_xlim(V_MIN, V_MAX)\n",
"ax.set_ylim(0.4, 1.01)\n",
"ax.grid(True, alpha=0.3)\n",
"ax.legend(loc='lower left')\n",
"plt.tight_layout()\n",
"plt.show()\n",
),
md("## Takeaways\n",
"\n",
"* **Fall-speed separation is the whole game.** At $\\phi = 30°$, SLW,\n",
" rain, and hail occupy disjoint velocity windows and spectral\n",
" polarimetry reads each species out directly — something bulk $Z_h$\n",
" + $Z_\\mathrm{dr}$ cannot do for this mixture. Lakshmi et al. (2024)\n",
" use exactly this separation throughout their paper: ice crystals,\n",
" aggregates, graupel, and rain+hail mixtures each stake out different\n",
" Doppler windows in their RHI-scan spectra (their Figs. 7, 8, 13, 15,\n",
" 19, 20).\n",
"* **Rain $sZ_\\mathrm{dr}$ rises monotonically with $v$** because fast\n",
" drops are big and oblate — the classic size-sorting diagnostic of\n",
" Kumjian & Ryzhkov (2012) and Wang et al. (2019).\n",
"* **Wet hail's C-band resonance notch** in $sZ_\\mathrm{dr}$ and\n",
" $s\\rho_\\mathrm{hv}$ at the fast-velocity end is the fingerprint of\n",
" melting $D \\sim$ 8–10 mm hailstones — exactly the rain+hail\n",
" signature Lakshmi et al. (2024) flag in their Fig. 5 case-study\n",
" region below the melting layer.\n",
"* **Mixture $sZ_\\mathrm{dr}$ collapses toward zero** at $v \\approx$\n",
" 5 m/s, where rain's positive $Z_\\mathrm{dr}$ and hail's\n",
" resonance-driven negative $Z_\\mathrm{dr}$ carry comparable power —\n",
" a direct mixed-species signature.\n",
"* **Per-species vs mixture plots** are the practical interpretive\n",
" tool: anywhere the mixture curve diverges from the dominant\n",
" single-species curve, two (or more) populations are contributing\n",
" coherently to that velocity bin.\n",
"* **Halving the hail concentration (gray dotted curve)** re-weights\n",
" the species contributions in a physically intuitive way. $sZ_h$\n",
" drops by $\\approx$3 dB at hail-dominated velocities but is barely\n",
" touched where rain or SLW dominate. $sK_\\mathrm{dp}$ hardly moves\n",
" at all, because rain — not hail — carries almost all of the\n",
" forward-scatter differential phase in this mixture.\n",
" $sZ_\\mathrm{dr}$ and $s\\rho_\\mathrm{hv}$ shift *toward* rain\n",
" wherever rain and hail carry comparable power — with less hail\n",
" power around $v \\approx$ 5 m/s the mixture $sZ_\\mathrm{dr}$ climbs\n",
" back toward rain's positive values and $s\\rho_\\mathrm{hv}$ rises\n",
" toward unity. This is exactly the lever Lakshmi et al. (2024)\n",
" exploit when they use spectral polarimetry to infer the\n",
" *proportion* of ice-phase scatterers within a resolution volume.\n"),
]
NB13 = [
md(
"# Tutorial 13 — Wind + turbulence sensitivity of Doppler moments\n",
"\n",
"A vertically pointing radar sees a Doppler spectrum broadened by\n",
"two independent mechanisms:\n",
"\n",
"* **Turbulence** — eddies smaller than the resolution volume\n",
" rearrange drop velocities, convolving the fall-speed spectrum\n",
" with a Gaussian of width $\\sigma_t$.\n",
"* **Horizontal wind through a finite beam** — any pixel off the\n",
" boresight contributes a radial component $u_h\\sin\\theta\\cos(\\phi-\\phi_w)$\n",
" to the line-of-sight. Integrated over a Gaussian beam of one-way\n",
" HPBW $\\theta_b$ this produces a *deterministic* Gaussian\n",
" broadening of width\n",
" \n",
" $$\\sigma_\\mathrm{beam} = \\frac{|u_h|\\,\\theta_b}{2\\sqrt{2\\ln 2}}$$\n",
" \n",
" (Doviak & Zrnić 1993, §5.3). It is the Doppler equivalent of a\n",
" standard beam-filling correction.\n",
"\n",
"The two widths add in quadrature: $\\sigma^2 = \\sigma_t^2 + \\sigma_\\mathrm{beam}^2$.\n",
"Neither biases the *first* moment (both are symmetric about the\n",
"boresight); both inflate the *second* moment.\n",
"\n",
"This notebook sweeps $u_h \\in \\{0, 5, 10, 20\\}$ m/s and\n",
"$\\theta_b \\in \\{1°, 3°, 5°\\}$ at fixed $\\sigma_t^2 = 0.5$ m²/s²\n",
"($\\sigma_t \\approx 0.707$ m/s) for a W-band vertical-pointing\n",
"radar, and tabulates the observed moments against the closed-form\n",
"$\\sigma_\\mathrm{beam}$ prediction.\n",
"\n",
"**Scene assumption** — horizontal wind is constant in magnitude\n",
"over the beam (no shear, no spatial gradient). Tutorial 14 takes\n",
"up the spatially heterogeneous case.\n"),
code(
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"\n",
"from rustmatrix import Scatterer, SpectralIntegrator, spectra\n",
"from rustmatrix.psd import ExponentialPSD, PSDIntegrator\n",
"from rustmatrix.refractive import m_w_10C\n",
"from rustmatrix.tmatrix_aux import (K_w_sqr, dsr_thurai_2007,\n",
" geom_vert_back, wl_W)\n",
"\n",
"BEAMWIDTHS_DEG = (1.0, 3.0, 5.0)\n",
"SIGMA_T_SQ = 0.5\n",
"SIGMA_T = float(np.sqrt(SIGMA_T_SQ))\n",
"U_H_LIST = (0.0, 5.0, 10.0, 20.0)\n",
"V_MIN, V_MAX, N_BINS = -1.0, 14.0, 1024\n",
"FWHM = 2.0 * np.sqrt(2.0 * np.log(2.0))\n"),
md("## W-band rain scatterer (Marshall-Palmer DSD)\n"),
code(
"rain = Scatterer(wavelength=wl_W, m=m_w_10C[wl_W],\n",
" Kw_sqr=K_w_sqr[wl_W], ddelt=1e-4, ndgs=2)\n",
"integ = PSDIntegrator()\n",
"integ.D_max = 5.0\n",
"integ.num_points = 64\n",
"integ.axis_ratio_func = lambda D: 1.0 / dsr_thurai_2007(D)\n",
"integ.geometries = (geom_vert_back,)\n",
"rain.psd_integrator = integ\n",
"rain.psd_integrator.init_scatter_table(rain)\n",
"rain.psd = ExponentialPSD(N0=8000.0, Lambda=2.2, D_max=5.0)\n"),
md("## Sweep $u_h$ and $\\theta_b$\n",
"We first compute a reference spectrum with no wind and a pencil beam\n",
"(so only the turbulence Gaussian broadens the native DSD spectrum),\n",
"then quadrature-add $\\sigma_\\mathrm{beam}$ for every $(u_h, \\theta_b)$\n",
"combination.\n"),
code(
"def moments(v, sZh):\n",
" sZh = np.clip(sZh, 0.0, None)\n",
" P = sZh.sum()\n",
" if P == 0:\n",
" return np.nan, np.nan\n",
" mu = float((v * sZh).sum() / P)\n",
" var = float(((v - mu) ** 2 * sZh).sum() / P)\n",
" return mu, float(np.sqrt(max(var, 0.0)))\n",
"\n",
"def run_case(u_h, hpbw_rad):\n",
" return SpectralIntegrator(\n",
" rain,\n",
" fall_speed=spectra.brandes_et_al_2002,\n",
" turbulence=spectra.GaussianTurbulence(sigma_t=SIGMA_T),\n",
" v_min=V_MIN, v_max=V_MAX, n_bins=N_BINS,\n",
" u_h=u_h, beamwidth=hpbw_rad,\n",
" ).run()\n",
"\n",
"r_ref = run_case(0.0, 0.0)\n",
"mu_ref, sig_ref = moments(r_ref.v, r_ref.sZ_h)\n",
"print(f'reference: mu = {mu_ref:.3f} m/s, sigma = {sig_ref:.3f} m/s')\n"),
code(
"results = {}\n",
"for theta_deg in BEAMWIDTHS_DEG:\n",
" theta = np.deg2rad(theta_deg)\n",
" for u_h in U_H_LIST:\n",
" r = run_case(u_h, theta)\n",
" mu, sig = moments(r.v, r.sZ_h)\n",
" sig_beam = u_h * theta / FWHM\n",
" sig_pred = float(np.sqrt(sig_ref ** 2 + sig_beam ** 2))\n",
" results[(theta_deg, u_h)] = dict(r=r, mu=mu, sig=sig,\n",
" sig_beam=sig_beam, sig_pred=sig_pred)\n",
"\n",
"header = f'{\"theta_b\":>7} {\"u_h\":>5} {\"sigma_beam\":>10} {\"sigma_pred\":>10} {\"mu_obs\":>7} {\"sigma_obs\":>9}'\n",
"print(header)\n",
"print('-' * len(header))\n",
"for (theta_deg, u_h), r in results.items():\n",
" print(f'{theta_deg:6.1f}° {u_h:5.1f} {r[\"sig_beam\"]:10.3f} '\n",
" f'{r[\"sig_pred\"]:10.3f} {r[\"mu\"]:7.3f} {r[\"sig\"]:9.3f}')\n"),
md("## Plot the spectra and the width-scaling curve\n",
"Left — $sZ_h(v)$ for $u_h = 0$ vs $20$ m/s at each beamwidth: the\n",
"wind broadens the fall-speed peak without shifting it. Right — the\n",
"observed spectral width traced against the analytic\n",
"$\\sqrt{\\sigma_t^2 + \\sigma_\\mathrm{beam}^2}$ prediction.\n"),
code(
"colors = {1.0: 'tab:blue', 3.0: 'tab:orange', 5.0: 'tab:red'}\n",
"fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 4.5))\n",
"\n",
"for theta_deg in BEAMWIDTHS_DEG:\n",
" r_quiet = results[(theta_deg, 0.0)]['r']\n",
" r_wind = results[(theta_deg, 20.0)]['r']\n",
" c = colors[theta_deg]\n",
" ax1.semilogy(r_quiet.v, np.clip(r_quiet.sZ_h, 1e-6, None),\n",
" c=c, ls='--', alpha=0.7,\n",
" label=f'θ_b = {theta_deg:.0f}°, u_h = 0')\n",
" ax1.semilogy(r_wind.v, np.clip(r_wind.sZ_h, 1e-6, None),\n",
" c=c, lw=1.8,\n",
" label=f'θ_b = {theta_deg:.0f}°, u_h = 20 m/s')\n",
"\n",
"ax1.set_xlim(0, 10)\n",
"ax1.set_xlabel('Doppler velocity [m/s]')\n",
"ax1.set_ylabel('$sZ_h$ [mm$^6$/m$^3$ / (m/s)]')\n",
"ax1.set_title('Spectra: wind broadens, does not shift')\n",
"ax1.legend(fontsize=8)\n",
"ax1.grid(alpha=0.3)\n",
"\n",
"for theta_deg in BEAMWIDTHS_DEG:\n",
" uhs = np.array(U_H_LIST)\n",
" obs = np.array([results[(theta_deg, u)]['sig'] for u in uhs])\n",
" pred = np.array([results[(theta_deg, u)]['sig_pred'] for u in uhs])\n",
" c = colors[theta_deg]\n",
" ax2.plot(uhs, pred, c=c, lw=1.0, alpha=0.6,\n",
" label=f'θ_b = {theta_deg:.0f}° (predicted)')\n",
" ax2.plot(uhs, obs, 'o', c=c, ms=8,\n",
" label=f'θ_b = {theta_deg:.0f}° (observed)')\n",
"ax2.axhline(SIGMA_T, c='k', ls=':', alpha=0.4, label='σ_t only')\n",
"ax2.set_xlabel('$u_h$ [m/s]')\n",
"ax2.set_ylabel('Spectral width σ [m/s]')\n",
"ax2.set_title('σ scales linearly with $u_h\\\\cdot\\\\theta_b$')\n",
"ax2.legend(fontsize=8)\n",
"ax2.grid(alpha=0.3)\n",
"\n",
"fig.tight_layout();\n"),
md("## Doppler velocity is unaffected by wind\n",
"The first moment of every spectrum is within a fraction of a\n",
"velocity bin of the reference value — horizontal wind in a\n",
"symmetric beam cannot bias $\\mu$. The second moment, by contrast,\n",
"scales exactly as the Doviak–Zrnić formula predicts.\n"),
code(
"fig, ax = plt.subplots(figsize=(7, 4))\n",
"for theta_deg in BEAMWIDTHS_DEG:\n",
" uhs = np.array(U_H_LIST)\n",
" mus = np.array([results[(theta_deg, u)]['mu'] for u in uhs])\n",
" ax.plot(uhs, mus - mu_ref, 'o-', c=colors[theta_deg],\n",
" label=f'θ_b = {theta_deg:.0f}°')\n",
"ax.axhline(0, c='k', ls=':', alpha=0.5)\n",
"ax.set_xlabel('$u_h$ [m/s]')\n",
"ax.set_ylabel(r'$\\mu - \\mu_\\mathrm{ref}$ [m/s]')\n",
"ax.set_title('First-moment bias is zero within grid discretisation')\n",
"ax.legend()\n",
"ax.grid(alpha=0.3);\n"),
md("## Take-aways\n",
"* **Mean Doppler velocity is insensitive to horizontal wind** when\n",
" the beam is circularly symmetric and the scene is uniform — the\n",
" contributions from the $+\\phi$ and $-\\phi$ sides of the beam\n",
" cancel exactly.\n",
"* **Spectral width increases in quadrature** with\n",
" $\\sigma_\\mathrm{beam} = u_h \\theta_b / (2\\sqrt{2\\ln 2})$. The\n",
" observed widths match the analytic prediction to the velocity-grid\n",
" resolution.\n",
"* **Narrow beams are wind-immune.** A 1° cloud radar adds only\n",
" 0.15 m/s to $\\sigma$ even at $u_h = 20$ m/s — well below the\n",
" intrinsic DSD width. A 5° beam doubles the apparent width at the\n",
" same wind, enough to mis-attribute a retrieval.\n",
"* **Scene structure changes the story** — when the beam straddles a\n",
" reflectivity gradient or a sheared updraft, the closed form breaks\n",
" down and the moments develop real biases. That is Tutorial 14.\n"),
]
NB14 = [
md(
"# Tutorial 14 — Beam pattern × scene integration (W-band down-looking)\n",
"\n",
"A radar does not sample its boresight pixel; it samples a\n",
"pattern-weighted integral over its solid angle. When the scene is\n",
"homogeneous, the closed-form $\\sigma_\\mathrm{beam}$ used in\n",
"Tutorial 13 captures everything. When the scene has structure on\n",
"scales comparable to (or finer than) the beam footprint — convective\n",
"cells, updraft/downdraft couplets, reflectivity gradients — that\n",
"closed form fails and the beam has to be integrated explicitly.\n",
"\n",
"Two properties of the pattern matter independently:\n",
"\n",
"* **Main-lobe width.** A 1° beam at 15 km range has a ≈ 260 m\n",
" footprint; a 3° beam has ≈ 780 m. Features narrower than the\n",
" footprint are smeared.\n",
"* **Sidelobes.** A uniform circular aperture has an Airy first\n",
" sidelobe at **−17.6 dB**. A distant bright cell sitting in that\n",
" sidelobe can dominate the moments if the main lobe is pointed at a\n",
" quiet patch.\n",
"\n",
"This notebook uses the new `rustmatrix.spectra.beam` module. A\n",
"down-looking radar at 20 km altitude scans across a synthetic rain\n",
"scene: 20 dBZ Marshall-Palmer background with 500-m-wide 45 dBZ\n",
"cells spaced every 1.5 km. Each cell carries a co-located vertical\n",
"motion — three combinations are explored:\n",
"\n",
"* **uniform_updraft** — every cell is a 3 m/s updraft.\n",
"* **alternating** — adjacent cells alternate −3 / +3 m/s.\n",
"* **dipole_couplet** — each cell is an updraft/downdraft dipole\n",
" straddling the enhanced-Z peak.\n",
"\n",
"Each scene is sampled with 1° and 3° beams, in both Gaussian and\n",
"Airy patterns. Horizontal wind is zero throughout — we isolate the\n",
"*scene-structure* contribution to beam broadening here; Tutorial 13\n",
"covered the uniform-wind case.\n"),
code(
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"\n",
"from rustmatrix import Scatterer\n",
"from rustmatrix.psd import PSDIntegrator\n",
"from rustmatrix.refractive import m_w_10C\n",
"from rustmatrix.spectra import brandes_et_al_2002\n",
"from rustmatrix.spectra.beam import (AiryBeam, BeamIntegrator, GaussianBeam,\n",
" Scene, marshall_palmer_psd_factory)\n",
"from rustmatrix.tmatrix_aux import (K_w_sqr, dsr_thurai_2007,\n",
" geom_vert_back, wl_W)\n",
"\n",
"RADAR_ALTITUDE_M = 20000.0\n",
"TARGET_ALTITUDE_M = 3000.0\n",
"RANGE_M = RADAR_ALTITUDE_M - TARGET_ALTITUDE_M\n",
"RAIN_TOP_M = 5000.0\n",
"\n",
"BG_DBZ = 20.0\n",
"CELL_PEAK_DBZ = 45.0\n",
"CELL_WIDTH_M = 500.0\n",
"CELL_SIGMA_M = CELL_WIDTH_M / (2.0 * np.sqrt(2.0 * np.log(2.0)))\n",
"CELL_CENTERS_X = np.arange(-3000.0, 3001.0, 1500.0)\n",
"W_PEAK = 3.0\n",
"\n",
"SCAN_X = np.linspace(-4000.0, 4000.0, 81)\n",
"V_MIN, V_MAX, N_BINS = -5.0, 15.0, 384\n"),
md("## Beam-pattern comparison\n",
"Before running the full scene sweep, look at the two candidate\n",
"patterns at identical HPBW. The Gaussian taper has no sidelobes;\n",
"the Airy pattern has nulls and a first sidelobe at −17.6 dB.\n"),
code(
"theta = np.linspace(0, np.deg2rad(6), 600)\n",
"gb = GaussianBeam(hpbw=np.deg2rad(1.0))\n",
"ab = AiryBeam(hpbw=np.deg2rad(1.0))\n",
"\n",
"fig, ax = plt.subplots(figsize=(7, 4))\n",
"ax.plot(np.rad2deg(theta), 10 * np.log10(np.clip(gb.gain(theta), 1e-6, None)),\n",
" label='Gaussian', lw=1.5)\n",
"ax.plot(np.rad2deg(theta), 10 * np.log10(np.clip(ab.gain(theta), 1e-6, None)),\n",
" label='Airy (uniform circular aperture)', lw=1.5)\n",
"ax.axhline(-3, c='k', ls=':', alpha=0.4, label='−3 dB (HPBW)')\n",
"ax.axhline(-17.6, c='r', ls=':', alpha=0.4, label='−17.6 dB (Airy 1st sidelobe)')\n",
"ax.set_xlabel('off-axis angle θ [°]')\n",
"ax.set_ylabel('normalized gain [dB]')\n",
"ax.set_ylim(-50, 1)\n",
"ax.set_title('Gaussian vs Airy one-way power pattern (HPBW = 1°)')\n",
"ax.legend()\n",
"ax.grid(alpha=0.3);\n"),
md("## Build the scene\n",
"Each scene is a triple of callables `(Z_dBZ, w, u_h)` evaluated at\n",
"pixel positions. We bake in the rain top at 5 km (above that the\n",
"scene is empty), the cell grid, and the three vertical-motion\n",
"variants.\n"),
code(
"def make_scene(pattern):\n",
" centers = CELL_CENTERS_X\n",
" sigma = CELL_SIGMA_M\n",
" z_top = RAIN_TOP_M\n",
" Z_bg_lin = 10.0 ** (BG_DBZ / 10.0)\n",
" Z_peak_excess = 10.0 ** (CELL_PEAK_DBZ / 10.0) - Z_bg_lin\n",
"\n",
" def Z_dBZ(x, y, z):\n",
" mask = (z >= 0) & (z <= z_top)\n",
" Z_lin = np.where(mask, Z_bg_lin, 1e-10)\n",
" for xc in centers:\n",
" bump = Z_peak_excess * np.exp(-0.5 * ((x - xc) / sigma) ** 2)\n",
" Z_lin = Z_lin + bump * mask\n",
" return 10.0 * np.log10(np.maximum(Z_lin, 1e-10))\n",
"\n",
" if pattern == 'uniform_updraft':\n",
" signs = -np.ones(len(centers))\n",
" elif pattern == 'alternating':\n",
" signs = np.array([-1.0 if i % 2 == 0 else 1.0\n",
" for i in range(len(centers))])\n",
" elif pattern == 'dipole_couplet':\n",
" signs = None\n",
" else:\n",
" raise ValueError(pattern)\n",
"\n",
" def w_fn(x, y, z):\n",
" mask = (z >= 0) & (z <= z_top)\n",
" w_total = np.zeros_like(x)\n",
" for i, xc in enumerate(centers):\n",
" arg = (x - xc) / sigma\n",
" if signs is None:\n",
" w_total = w_total + W_PEAK * arg * np.exp(0.5 - 0.5 * arg ** 2)\n",
" else:\n",
" w_total = w_total + signs[i] * W_PEAK * np.exp(-0.5 * arg ** 2)\n",
" return w_total * mask\n",
"\n",
" def u_h_fn(x, y, z):\n",
" return np.zeros_like(x)\n",
"\n",
" return Scene(Z_dBZ=Z_dBZ, w=w_fn, u_h=u_h_fn, u_h_azimuth=0.0)\n"),
code(
"# Visualize each scene along the x-axis at z = TARGET_ALTITUDE_M.\n",
"x_vis = np.linspace(-4000, 4000, 801)\n",
"zero = np.zeros_like(x_vis)\n",
"z_target = np.full_like(x_vis, TARGET_ALTITUDE_M)\n",
"\n",
"fig, (axZ, axW) = plt.subplots(2, 1, figsize=(10, 5.5), sharex=True)\n",
"for name in ('uniform_updraft', 'alternating', 'dipole_couplet'):\n",
" sc = make_scene(name)\n",
" axZ.plot(x_vis, sc.Z_dBZ(x_vis, zero, z_target), label=name, lw=1.3)\n",
" axW.plot(x_vis, sc.w(x_vis, zero, z_target), label=name, lw=1.3)\n",
"axZ.set_ylabel('Z [dBZ]')\n",
"axZ.axhline(BG_DBZ, c='k', ls=':', alpha=0.3, label='bg')\n",
"axZ.legend(fontsize=8)\n",
"axZ.grid(alpha=0.3)\n",
"axW.set_xlabel('x [m]')\n",
"axW.set_ylabel('w (pos. = down) [m/s]')\n",
"axW.axhline(0, c='k', ls=':', alpha=0.3)\n",
"axW.grid(alpha=0.3);\n"),
md("## W-band rain scatterer and the scan sweep\n",
"The `BeamIntegrator` takes the scatterer, beam pattern, scene, and\n",
"a PSD factory (Z → PSD mapping — Marshall–Palmer here) and returns\n",
"a `SpectralResult` identical to the one produced by the ordinary\n",
"`SpectralIntegrator`. We run it once per scan position.\n"),
code(
"rain = Scatterer(wavelength=wl_W, m=m_w_10C[wl_W],\n",
" Kw_sqr=K_w_sqr[wl_W], ddelt=1e-4, ndgs=2)\n",
"integ = PSDIntegrator()\n",
"integ.D_max = 5.0\n",
"integ.num_points = 48\n",
"integ.axis_ratio_func = lambda D: 1.0 / dsr_thurai_2007(D)\n",
"integ.geometries = (geom_vert_back,)\n",
"rain.psd_integrator = integ\n",
"rain.psd_integrator.init_scatter_table(rain)\n",
"\n",
"psd_factory = marshall_palmer_psd_factory(N0=8000.0, D_max=5.0)\n"),
code(
"def moments(v, sZh):\n",
" sZh = np.clip(sZh, 0.0, None)\n",
" P = sZh.sum()\n",
" if P <= 0:\n",
" return -np.inf, np.nan, np.nan\n",
" dv = np.mean(np.diff(v))\n",
" Z_dBZ = 10.0 * np.log10(max(P * dv, 1e-10))\n",
" mu = float((v * sZh).sum() / P)\n",
" var = float(((v - mu) ** 2 * sZh).sum() / P)\n",
" return Z_dBZ, mu, float(np.sqrt(max(var, 0.0)))\n",
"\n",
"def sweep(scene, beam):\n",
" Zs = np.empty_like(SCAN_X)\n",
" mus = np.empty_like(SCAN_X)\n",
" sigs = np.empty_like(SCAN_X)\n",
" for i, xr in enumerate(SCAN_X):\n",
" bi = BeamIntegrator(\n",
" scatterer=rain, beam=beam, scene=scene,\n",
" psd_factory=psd_factory, fall_speed=brandes_et_al_2002,\n",
" radar_position=(xr, 0.0, RADAR_ALTITUDE_M),\n",
" boresight=(0.0, 0.0, -1.0), range_m=RANGE_M,\n",
" v_min=V_MIN, v_max=V_MAX, n_bins=N_BINS,\n",
" n_theta=16, n_phi=16,\n",
" )\n",
" r = bi.run()\n",
" Zs[i], mus[i], sigs[i] = moments(r.v, r.sZ_h)\n",
" return Zs, mus, sigs\n"),
code(
"PATTERNS = ('uniform_updraft', 'alternating', 'dipole_couplet')\n",
"BEAMS = (\n",
" ('gaussian', 1.0, GaussianBeam(hpbw=np.deg2rad(1.0))),\n",
" ('airy', 1.0, AiryBeam(hpbw=np.deg2rad(1.0))),\n",
" ('gaussian', 3.0, GaussianBeam(hpbw=np.deg2rad(3.0))),\n",
" ('airy', 3.0, AiryBeam(hpbw=np.deg2rad(3.0))),\n",
")\n",
"\n",
"results = {}\n",
"for pattern in PATTERNS:\n",
" scene = make_scene(pattern)\n",
" for kind, hpbw_deg, beam in BEAMS:\n",
" results[(pattern, kind, hpbw_deg)] = sweep(scene, beam)\n",
" print(f' {pattern:18s} done')\n"),
md("## Scan curves — moments vs radar x position\n",
"For each scene, plot $Z$, $V_R$ (Doppler velocity, first moment),\n",
"and $\\sigma$ (spectral width) as the radar sweeps across the cell\n",
"grid. Compare how narrow (1°) and wide (3°) beams resolve the\n",
"structure, and how the Airy sidelobes leak in neighbouring-cell\n",
"signal even when the main lobe is off a peak.\n"),
code(
"def plot_scan(pattern):\n",
" fig, axes = plt.subplots(3, 1, figsize=(10, 8), sharex=True)\n",
" styles = {('gaussian', 1.0): ('tab:blue', '-'),\n",
" ('airy', 1.0): ('tab:blue', '--'),\n",
" ('gaussian', 3.0): ('tab:red', '-'),\n",
" ('airy', 3.0): ('tab:red', '--')}\n",
" for (kind, hpbw_deg), (c, ls) in styles.items():\n",
" Zs, mus, sigs = results[(pattern, kind, hpbw_deg)]\n",
" lbl = f'{kind}, {hpbw_deg:.0f}°'\n",
" axes[0].plot(SCAN_X, Zs, c=c, ls=ls, lw=1.3, label=lbl)\n",
" axes[1].plot(SCAN_X, mus, c=c, ls=ls, lw=1.3)\n",
" axes[2].plot(SCAN_X, sigs, c=c, ls=ls, lw=1.3)\n",
" # mark cell centers\n",
" for xc in CELL_CENTERS_X:\n",
" for ax in axes:\n",
" ax.axvline(xc, c='k', ls=':', alpha=0.15)\n",
" axes[0].set_ylabel('Z [dBZ]')\n",
" axes[1].set_ylabel(r'$V_R$ [m/s]')\n",
" axes[2].set_ylabel(r'$\\sigma$ [m/s]')\n",
" axes[2].set_xlabel('radar x [m]')\n",
" axes[0].set_title(f'scene = {pattern}')\n",
" axes[0].legend(fontsize=8, ncol=2)\n",
" for ax in axes:\n",
" ax.grid(alpha=0.3)\n",
" fig.tight_layout()\n",
" return fig\n",
"\n",
"for p in PATTERNS:\n",
" plot_scan(p);\n"),
md("## Interpretation\n",
"* **Narrow (1°) vs wide (3°) main lobes.** The 1° beam's ~260 m\n",
" footprint is narrower than the 500 m cell, so $Z$ resolves the\n",
" full peak-to-trough swing ($\\approx$25 dB). The 3° beam averages\n",
" across neighbouring cells and the peak-to-trough contrast drops\n",
" to a few dB.\n",
"* **Gaussian vs Airy (same HPBW).** The main-lobe difference is\n",
" small — both patterns integrate a similar main-lobe footprint of\n",
" reflectivity. The *sidelobes* of the Airy pattern pick up signal\n",
" from cells up to $\\pm 2$ km away, leaking a small bias into the\n",
" inter-cell minima. The `alternating` and `dipole_couplet` scenes\n",
" show this most clearly: at the between-cell location, Airy's\n",
" $V_R$ is shifted toward the nearer cell's velocity while Gaussian\n",
" is closer to zero (the DSD-intrinsic rest frame).\n",
"* **uniform_updraft vs alternating.** Both start from the same Z\n",
" field, but in `alternating` adjacent cells pull $V_R$ in opposite\n",
" directions — when the 3° beam straddles a boundary it sees a\n",
" bimodal spectrum whose first moment can land near zero even\n",
" though both contributing populations are strongly non-still.\n",
" This is the physical origin of the \"turbulence\" signature in\n",
" spectral-width retrievals over convection.\n",
"* **dipole_couplet.** At the cell centre the beam averages equal\n",
" up- and downdraft power, giving $V_R \\approx v_t$ (DSD-only) and\n",
" an inflated $\\sigma$ — the classic bimodal-spectrum fingerprint.\n",
" The 3° beam smooths over more of the couplet, narrowing the\n",
" apparent $\\sigma$ swing.\n",
"* **Practical take.** Any moment retrieval that assumes a pencil\n",
" beam and homogeneous scene is subtracting the wrong amount of\n",
" \"beam broadening\" from the observed $\\sigma$. The\n",
" `BeamIntegrator` lets you attach a forward model to a scene and\n",
" an instrument specification, and returns spectra and moments\n",
" that the observed radar actually would see.\n"),
md("## Interactive exploration — cell spacing\n",
"The sweep above used 1.5 km spacing between cells. What happens\n",
"when that spacing shrinks below the beam footprint, or grows to\n",
"where even a wide beam resolves the cells cleanly? The slider\n",
"below varies the spacing from **100 m** (far finer than either\n",
"beam footprint) to **10 km** (fully resolved by both). At fine\n",
"spacings the `alternating` pattern's $V_R$ collapses toward zero\n",
"— adjacent up- and downdrafts average inside the pattern — while\n",
"$\\sigma$ balloons because the spectrum becomes bimodal. At coarse\n",
"spacings the 1° and 3° beams converge on the same answer because\n",
"both now resolve the scene.\n",
"\n",
"This widget drives the same `BeamIntegrator` used above, but on a\n",
"coarser scan (41 x-positions, 10×10 beam samples) so each slider\n",
"change completes in a few seconds.\n"),
code(
"import ipywidgets as widgets\n",
"\n",
"FAST_SCAN_X = np.linspace(-4000.0, 4000.0, 41)\n",
"\n",
"def _scene_at_spacing(pattern, spacing_m):\n",
" centers = np.arange(-4000.0, 4000.0 + 1e-6, spacing_m)\n",
" sigma = CELL_SIGMA_M\n",
" z_top = RAIN_TOP_M\n",
" Z_bg_lin = 10.0 ** (BG_DBZ / 10.0)\n",
" Z_peak_excess = 10.0 ** (CELL_PEAK_DBZ / 10.0) - Z_bg_lin\n",
"\n",
" def Z_dBZ(x, y, z):\n",
" mask = (z >= 0) & (z <= z_top)\n",
" Z_lin = np.where(mask, Z_bg_lin, 1e-10)\n",
" for xc in centers:\n",
" bump = Z_peak_excess * np.exp(-0.5 * ((x - xc) / sigma) ** 2)\n",
" Z_lin = Z_lin + bump * mask\n",
" return 10.0 * np.log10(np.maximum(Z_lin, 1e-10))\n",
"\n",
" if pattern == 'uniform_updraft':\n",
" signs = -np.ones(len(centers))\n",
" elif pattern == 'alternating':\n",
" signs = np.array([-1.0 if i % 2 == 0 else 1.0\n",
" for i in range(len(centers))])\n",
" elif pattern == 'dipole_couplet':\n",
" signs = None\n",
" else:\n",
" raise ValueError(pattern)\n",
"\n",
" def w_fn(x, y, z):\n",
" mask = (z >= 0) & (z <= z_top)\n",
" w_total = np.zeros_like(x)\n",
" for i, xc in enumerate(centers):\n",
" arg = (x - xc) / sigma\n",
" if signs is None:\n",
" w_total = w_total + W_PEAK * arg * np.exp(0.5 - 0.5 * arg ** 2)\n",
" else:\n",
" w_total = w_total + signs[i] * W_PEAK * np.exp(-0.5 * arg ** 2)\n",
" return w_total * mask\n",
"\n",
" def u_h_fn(x, y, z):\n",
" return np.zeros_like(x)\n",
"\n",
" return Scene(Z_dBZ=Z_dBZ, w=w_fn, u_h=u_h_fn, u_h_azimuth=0.0), centers\n",
"\n",
"def _fast_sweep(scene, beam):\n",
" mus = np.empty_like(FAST_SCAN_X)\n",
" sigs = np.empty_like(FAST_SCAN_X)\n",
" for i, xr in enumerate(FAST_SCAN_X):\n",
" bi = BeamIntegrator(\n",
" scatterer=rain, beam=beam, scene=scene,\n",
" psd_factory=psd_factory, fall_speed=brandes_et_al_2002,\n",
" radar_position=(xr, 0.0, RADAR_ALTITUDE_M),\n",
" boresight=(0.0, 0.0, -1.0), range_m=RANGE_M,\n",
" v_min=V_MIN, v_max=V_MAX, n_bins=192,\n",
" n_theta=10, n_phi=10,\n",
" )\n",
" r = bi.run()\n",
" _, mus[i], sigs[i] = moments(r.v, r.sZ_h)\n",
" return mus, sigs\n",
"\n",
"def explore_spacing(spacing_m=1500.0, pattern='alternating'):\n",
" scene, centers = _scene_at_spacing(pattern, spacing_m)\n",
" g1 = GaussianBeam(hpbw=np.deg2rad(1.0))\n",
" g3 = GaussianBeam(hpbw=np.deg2rad(3.0))\n",
" mus1, sigs1 = _fast_sweep(scene, g1)\n",
" mus3, sigs3 = _fast_sweep(scene, g3)\n",
"\n",
" x_vis = np.linspace(-4000, 4000, 801)\n",
" zero = np.zeros_like(x_vis)\n",
" z_target = np.full_like(x_vis, TARGET_ALTITUDE_M)\n",
"\n",
" fig, axes = plt.subplots(3, 1, figsize=(10, 8), sharex=True)\n",
" axes[0].plot(x_vis, scene.Z_dBZ(x_vis, zero, z_target), 'k-', lw=1.0)\n",
" axes[0].set_ylabel('Z [dBZ]')\n",
" axes[1].plot(FAST_SCAN_X, mus1, color='tab:blue', lw=1.5, label='1° Gaussian')\n",
" axes[1].plot(FAST_SCAN_X, mus3, color='tab:red', lw=1.5, label='3° Gaussian')\n",
" axes[1].set_ylabel(r'$V_R$ [m/s]')\n",
" axes[1].legend(fontsize=8, loc='best')\n",
" axes[2].plot(FAST_SCAN_X, sigs1, color='tab:blue', lw=1.5)\n",
" axes[2].plot(FAST_SCAN_X, sigs3, color='tab:red', lw=1.5)\n",
" axes[2].set_ylabel(r'$\\sigma$ [m/s]')\n",
" axes[2].set_xlabel('radar x [m]')\n",
" for xc in centers:\n",
" for ax in axes:\n",
" ax.axvline(xc, color='k', ls=':', alpha=0.12)\n",
" for ax in axes:\n",
" ax.grid(alpha=0.3)\n",
" axes[0].set_title(f'pattern = {pattern}, spacing = {spacing_m:.0f} m, '\n",
" f'{len(centers)} cells')\n",
" fig.tight_layout()\n",
" plt.show()\n",
"\n",
"spacing_slider = widgets.FloatLogSlider(\n",
" value=1500.0, base=10,\n",
" min=np.log10(100.0), max=np.log10(10000.0),\n",
" step=0.02, description='spacing [m]',\n",
" readout_format='.0f', continuous_update=False,\n",
")\n",
"pattern_dd = widgets.Dropdown(\n",
" options=['uniform_updraft', 'alternating', 'dipole_couplet'],\n",
" value='alternating', description='pattern',\n",
")\n",
"widgets.interact(explore_spacing, spacing_m=spacing_slider, pattern=pattern_dd);\n"),
]
def main() -> None:
here = Path(__file__).parent
for name, cells in (("01_sphere_mie.ipynb", NB01),
("02_raindrop_zdr.ipynb", NB02),
("03_psd_gamma_rain.ipynb", NB03),
("04_oriented_ice.ipynb", NB04),
("05_radar_band_sweep.ipynb", NB05),
("06_hd_mix.ipynb", NB06),
("07_doppler_spectrum_rain.ipynb", NB07),
("08_spectral_polarimetry_rain_ice.ipynb", NB08),
("09_zhu_2023_particle_inertia.ipynb", NB09),
("10_slw_vs_snow.ipynb", NB10),
("11_honeyager_hydrometeor_classes.ipynb", NB11),
("12_spectral_polarimetry_rain_slw_hail.ipynb", NB12),
("13_wind_turbulence_sensitivity.ipynb", NB13),
("14_beam_pattern_scene.ipynb", NB14)):
path = here / name
with open(path, "w") as f:
json.dump(notebook(cells), f, indent=1)
print(f"wrote {path}")
if __name__ == "__main__":
main()