{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Tutorial 03 \u2014 Gamma-PSD rain at C-band\n",
    "\n",
    "Radar volumes sample thousands of drops. The observed Z_h, Z_dr,\n",
    "K_dp, and specific attenuation A_i are PSD-weighted integrals of\n",
    "the single-drop quantities. This notebook tabulates S(D) and Z(D)\n",
    "once, then sweeps the integrated observables across a range of\n",
    "normalised gamma PSDs parameterised by the median volume diameter\n",
    "D0.\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, radar, psd as rs_psd\n",
    "from rustmatrix.tmatrix_aux import (dsr_thurai_2007, geom_horiz_back,\n",
    "                                      geom_horiz_forw, K_w_sqr, wl_C)\n",
    "from rustmatrix.refractive import m_w_10C\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Tabulate once, evaluate many\n"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "s = Scatterer(wavelength=wl_C, m=m_w_10C[wl_C],\n",
    "              Kw_sqr=K_w_sqr[wl_C], ddelt=1e-4, ndgs=2)\n",
    "integ = rs_psd.PSDIntegrator()\n",
    "integ.D_max = 8.0\n",
    "integ.num_points = 64\n",
    "integ.axis_ratio_func = lambda D: 1.0 / dsr_thurai_2007(D)\n",
    "integ.geometries = (geom_horiz_back, geom_horiz_forw)\n",
    "s.psd_integrator = integ\n",
    "s.psd_integrator.init_scatter_table(s)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Sweep median diameter D0\n"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "D0s = np.linspace(0.5, 3.0, 12)\n",
    "Zh = np.empty_like(D0s)\n",
    "Zdr = np.empty_like(D0s)\n",
    "Kdp = np.empty_like(D0s)\n",
    "Ai = np.empty_like(D0s)\n",
    "\n",
    "for i, D0 in enumerate(D0s):\n",
    "    s.psd = rs_psd.GammaPSD(D0=D0, Nw=8e3, mu=4)\n",
    "    s.set_geometry(geom_horiz_back)\n",
    "    Zh[i] = 10 * np.log10(radar.refl(s))\n",
    "    Zdr[i] = 10 * np.log10(radar.Zdr(s))\n",
    "    s.set_geometry(geom_horiz_forw)\n",
    "    Kdp[i] = radar.Kdp(s)\n",
    "    Ai[i] = radar.Ai(s)\n",
    "\n",
    "fig, axes = plt.subplots(2, 2, figsize=(9, 6), sharex=True)\n",
    "axes[0, 0].plot(D0s, Zh, 'C0-o'); axes[0, 0].set_ylabel('Z_h [dBZ]')\n",
    "axes[0, 1].plot(D0s, Zdr, 'C1-o'); axes[0, 1].set_ylabel('Z_dr [dB]')\n",
    "axes[1, 0].plot(D0s, Kdp, 'C2-o'); axes[1, 0].set_ylabel('K_dp [\u00b0/km]')\n",
    "axes[1, 1].plot(D0s, Ai, 'C3-o'); axes[1, 1].set_ylabel('A_i [dB/km]')\n",
    "for ax in axes.flat:\n",
    "    ax.set_xlabel('D0 [mm]')\n",
    "    ax.grid(True, alpha=0.3)\n",
    "fig.suptitle('C-band gamma-PSD rain observables (Nw=8e3, mu=4)')\n",
    "fig.tight_layout();\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python",
   "version": "3.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}