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