{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Tutorial 02 \u2014 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\u2076 in the Rayleigh regime and then walks out of\n",
" it \u2014 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 \u2248 5 mm.\n",
"* **K_dp** scales with Re(f_h(0) \u2212 f_v(0)) \u2014 strictly stronger at\n",
" shorter wavelengths; X-band K_dp per drop is \u2248 2\u00d7 C-band and\n",
" \u2248 4\u00d7 S-band.\n",
"* **LDR** \u2014 linear depolarisation ratio \u2014 is set by the canting\n",
" distribution. Here we model a \u03c3 = 5\u00b0 Gaussian wobble around the\n",
" flat-lying orientation to produce realistic rain LDR in the\n",
" \u221230 to \u221225 dB range for 5+ mm drops.\n",
"\n",
"This notebook sweeps drop equivalent diameter D = 0.1\u20138 mm at S,\n",
"C, and X bands and plots all four observables. We report Z_h and\n",
"K_dp *per drop/m\u00b3* so multiplying by the drop concentration\n",
"N [m\u207b\u00b3] gives the usual bulk observables.\n"
]
},
{
"cell_type": "code",
"metadata": {},
"execution_count": null,
"outputs": [],
"source": [
"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"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Build the drop and run the sweep\n",
"\n",
"One `Scatterer` per (D, \u03bb) point, with a \u03c3 = 5\u00b0 Gaussian canting\n",
"PDF around \u03b2 = 0\u00b0 (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"
]
},
{
"cell_type": "code",
"metadata": {},
"execution_count": null,
"outputs": [],
"source": [
"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"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Plot all four observables vs. D\n",
"\n",
"Top-left Z_h tracks the D\u2076 line until it bends: X-band breaks\n",
"earliest (shortest \u03bb, smallest \u03c7 = \u03c0D/\u03bb needed), S-band last.\n",
"Top-right Z_dr rises monotonically; the C-band curve bumps above\n",
"X and S around D \u2248 5 mm \u2014 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"
]
},
{
"cell_type": "code",
"metadata": {},
"execution_count": null,
"outputs": [],
"source": [
"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\u2076 (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\u00b3]')\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| [\u00b0/km per drop/m\u00b3]')\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 \u00b0C water, \u03c3_canting = {CANTING_STD_DEG:.0f}\u00b0)')\n",
"fig.tight_layout();\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Spot values at canonical diameters\n",
"\n",
"Six diameters span the regime map: 0.5 mm (nearly spherical,\n",
"pure Rayleigh), 1\u20133 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"
]
},
{
"cell_type": "code",
"metadata": {},
"execution_count": null,
"outputs": [],
"source": [
"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 [\u00b0/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"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"name": "python",
"version": "3.11"
}
},
"nbformat": 4,
"nbformat_minor": 5
}