{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Tutorial 13 \u2014 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** \u2014 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** \u2014 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\u0107 1993, \u00a75.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\u00b0, 3\u00b0, 5\u00b0\\}$ at fixed $\\sigma_t^2 = 0.5$ m\u00b2/s\u00b2\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** \u2014 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"
]
},
{
"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, 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"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## W-band rain scatterer (Marshall-Palmer DSD)\n"
]
},
{
"cell_type": "code",
"metadata": {},
"execution_count": null,
"outputs": [],
"source": [
"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"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 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"
]
},
{
"cell_type": "code",
"metadata": {},
"execution_count": null,
"outputs": [],
"source": [
"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"
]
},
{
"cell_type": "code",
"metadata": {},
"execution_count": null,
"outputs": [],
"source": [
"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}\u00b0 {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"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Plot the spectra and the width-scaling curve\n",
"Left \u2014 $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 \u2014 the\n",
"observed spectral width traced against the analytic\n",
"$\\sqrt{\\sigma_t^2 + \\sigma_\\mathrm{beam}^2}$ prediction.\n"
]
},
{
"cell_type": "code",
"metadata": {},
"execution_count": null,
"outputs": [],
"source": [
"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'\u03b8_b = {theta_deg:.0f}\u00b0, 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'\u03b8_b = {theta_deg:.0f}\u00b0, 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'\u03b8_b = {theta_deg:.0f}\u00b0 (predicted)')\n",
" ax2.plot(uhs, obs, 'o', c=c, ms=8,\n",
" label=f'\u03b8_b = {theta_deg:.0f}\u00b0 (observed)')\n",
"ax2.axhline(SIGMA_T, c='k', ls=':', alpha=0.4, label='\u03c3_t only')\n",
"ax2.set_xlabel('$u_h$ [m/s]')\n",
"ax2.set_ylabel('Spectral width \u03c3 [m/s]')\n",
"ax2.set_title('\u03c3 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"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 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 \u2014 horizontal wind in a\n",
"symmetric beam cannot bias $\\mu$. The second moment, by contrast,\n",
"scales exactly as the Doviak\u2013Zrni\u0107 formula predicts.\n"
]
},
{
"cell_type": "code",
"metadata": {},
"execution_count": null,
"outputs": [],
"source": [
"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'\u03b8_b = {theta_deg:.0f}\u00b0')\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"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Take-aways\n",
"* **Mean Doppler velocity is insensitive to horizontal wind** when\n",
" the beam is circularly symmetric and the scene is uniform \u2014 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\u00b0 cloud radar adds only\n",
" 0.15 m/s to $\\sigma$ even at $u_h = 20$ m/s \u2014 well below the\n",
" intrinsic DSD width. A 5\u00b0 beam doubles the apparent width at the\n",
" same wind, enough to mis-attribute a retrieval.\n",
"* **Scene structure changes the story** \u2014 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"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"name": "python",
"version": "3.11"
}
},
"nbformat": 4,
"nbformat_minor": 5
}