{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Tutorial 04 \u2014 Oriented ice crystals at W-band\n",
    "\n",
    "Ice crystals fall with a preferred orientation and a spread of\n",
    "canting angles around it. The orientation PDF captures the spread\n",
    "and feeds directly into the dual-pol observables. This notebook\n",
    "compares the three orientation-averaging strategies rustmatrix\n",
    "ships with: none, fixed-quadrature, and adaptive.\n"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "import time\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from rustmatrix import Scatterer, orientation as rs_orient, radar\n",
    "from rustmatrix.tmatrix_aux import geom_horiz_back, wl_W\n",
    "from rustmatrix.refractive import mi\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## One crystal, three averaging schemes\n"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "D_eq_mm = 0.5\n",
    "ice_m = mi(wl_W, 0.9)\n",
    "pdf = rs_orient.gaussian_pdf(std=20.0, mean=90.0)\n",
    "base = dict(radius=D_eq_mm/2, wavelength=wl_W, m=ice_m,\n",
    "            axis_ratio=0.5, ddelt=1e-4, ndgs=2)\n",
    "\n",
    "schemes = [\n",
    "    ('single', rs_orient.orient_single, {'or_pdf': pdf}),\n",
    "    ('fixed 6\u00d712', rs_orient.orient_averaged_fixed,\n",
    "     {'or_pdf': pdf, 'n_alpha': 6, 'n_beta': 12}),\n",
    "    ('adaptive', rs_orient.orient_averaged_adaptive, {'or_pdf': pdf}),\n",
    "]\n",
    "rows = []\n",
    "for name, orient, extra in schemes:\n",
    "    s = Scatterer(**base, **extra)\n",
    "    s.orient = orient\n",
    "    s.set_geometry(geom_horiz_back)\n",
    "    t0 = time.perf_counter()\n",
    "    zdr = radar.Zdr(s)\n",
    "    elapsed = time.perf_counter() - t0\n",
    "    rows.append((name, 10*np.log10(zdr), elapsed))\n",
    "    print(f'{name:<12} Z_dr = {10*np.log10(zdr):+.4f} dB   ({elapsed*1000:.1f} ms)')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Canting-angle spread sensitivity\n",
    "\n",
    "Sweep the Gaussian PDF's `std` parameter to see how rapidly Z_dr\n",
    "collapses as the orientation spread widens.\n"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "stds = np.array([1, 5, 10, 15, 20, 30, 45, 60])\n",
    "zdrs = np.empty_like(stds, dtype=float)\n",
    "for i, sd in enumerate(stds):\n",
    "    s = Scatterer(**base)\n",
    "    s.orient = rs_orient.orient_averaged_fixed\n",
    "    s.or_pdf = rs_orient.gaussian_pdf(std=float(sd), mean=90.0)\n",
    "    s.n_alpha, s.n_beta = 6, 12\n",
    "    s.set_geometry(geom_horiz_back)\n",
    "    zdrs[i] = 10 * np.log10(radar.Zdr(s))\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(7, 4))\n",
    "ax.plot(stds, zdrs, 'C0-o')\n",
    "ax.set_xlabel('canting-angle \u03c3 [deg]')\n",
    "ax.set_ylabel('Z_dr [dB]')\n",
    "ax.set_title('W-band prolate ice column, oriented')\n",
    "ax.grid(True, alpha=0.3)\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
}