{
 "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
}