{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Tutorial 08 \u2014 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\u2013931\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 \u2014 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\u00b7log\u2081\u2080(sZ_X / sZ_W) stays near 0\n",
    "at low velocities and rises to several dB at high velocities \u2014 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"
   ]
  },
  {
   "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, 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"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 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, \u03c1_ice = 0.1 g/cm\u00b3 (low-density, typical of\n",
    "  mixed-habit aggregates), axis ratio 0.6, D_max = 5 mm\n",
    "* exponential PSD, N\u2080 = 2\u00d710\u2074 m\u207b\u00b3 mm\u207b\u00b9, \u039b = 2.5 mm\u207b\u00b9\n",
    "* median-volume diameter D\u2080 = 3.67/\u039b \u2248 1.5 mm\n",
    "* IWC = (\u03c0 \u03c1_ice / \u039b\u2074) \u00b7 N\u2080 \u2248 0.16 g/m\u00b3 \u2014 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\u201320 dB \u2014 *not* the\n",
    "50 dBZ / 25 dB outlier you get if you blindly pick \u039b = 0.8 mm\u207b\u00b9 and\n",
    "D_max = 10 mm (D\u2080 \u2248 5 mm, IWC > 10 g/m\u00b3 \u2014 a reflectivity-saturating\n",
    "deep convective snow cell, not an aggregation layer).\n"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "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"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## \u03c3_b(D) at X-band vs W-band \u2014 the non-Rayleigh onset\n",
    "\n",
    "Up to about 2 mm equivalent diameter the two \u03c3_b(D) curves are\n",
    "essentially parallel D\u2076 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 \u2014 this is the imprint the dual-frequency spectrum\n",
    "will show up as a velocity-resolved sDWR.\n"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "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\u00b2]')\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"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Dual-frequency Doppler spectra\n",
    "\n",
    "Same fall-speed model (Locatelli\u2013Hobbs aggregates), same turbulence,\n",
    "same velocity grid \u2014 only the radar wavelength changes.\n"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "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 \u2014 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"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 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 \u2014 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"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "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"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python",
   "version": "3.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}