{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Tutorial 10 \u2014 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\u201310234 (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 (\u2272200 \u00b5m), spherical, cold refractive index of\n",
    "  water, falls at a few cm/s. Rayleigh at both X- and W-band \u21d2 very\n",
    "  low Z_h, Z_dr \u2248 0 dB, no polarimetric signature.\n",
    "* **Snow**: millimetre-scale oblate low-density aggregates, falls at\n",
    "  ~0.5\u20131 m/s. Moderately non-Rayleigh at W-band \u21d2 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 \u03c3_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"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "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"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Build a supercooled-liquid-water scatterer\n",
    "\n",
    "Small spherical drops (D in 10\u2013200 \u00b5m), refractive index at 0\u00b0C,\n",
    "log-normal-ish gamma PSD centred around 30 \u00b5m. We use mm units\n",
    "throughout to match `rustmatrix` conventions.\n"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "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 \u00b5m), Nw = 1e11 m\u207b\u00b3 mm\u207b\u00b9, \u03bc = 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"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Build a snow scatterer (low-density oblate aggregates)\n"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "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"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## \u03c3_b(D) \u2014 Rayleigh vs. Mie across the two populations\n"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "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\u00b2]')\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"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 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 \u2014 the D\u2076 Rayleigh penalty at 30-\u00b5m scale kills\n",
    "the signal. Z_dr is zero for SLW (spheres) but positive for snow\n",
    "(oblate). Dual-wavelength ratio DWR = Z_X \u2212 Z_W is near zero for SLW\n",
    "and positive for snow (non-Rayleigh at W-band).\n"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "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"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 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 \u2248 0 with a narrow spread,\n",
    "snow sits near 0.5\u20131.0 m/s with a broader tail. Plotted in dBZ space\n",
    "the two modes are clearly separable even with realistic noise.\n"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# Use W-band for the spectrum \u2014 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"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 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"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python",
   "version": "3.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}