{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Tutorial 01 \u2014 A dielectric sphere at X-band, checked against Mie\n",
    "\n",
    "A T-matrix solver for nonspherical particles must reduce to classical\n",
    "Mie theory in the `axis_ratio = 1` limit. That reduction is the\n",
    "simplest end-to-end check we can do. This tutorial builds a 1 mm\n",
    "water sphere at X-band, computes its scattering and extinction\n",
    "cross-sections via the T-matrix path, and compares against the\n",
    "closed-form Mie expressions shipped with rustmatrix.\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, mie_qsca, mie_qext, scatter\n",
    "from rustmatrix.tmatrix_aux import geom_horiz_back, wl_X\n",
    "from rustmatrix.refractive import m_w_10C\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Build a sphere scatterer at X-band\n"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "radius_mm = 1.0\n",
    "wavelength_mm = wl_X\n",
    "m = m_w_10C[wl_X]\n",
    "\n",
    "s = Scatterer(radius=radius_mm, wavelength=wavelength_mm, m=m,\n",
    "              axis_ratio=1.0, ddelt=1e-4, ndgs=2)\n",
    "s.set_geometry(geom_horiz_back)\n",
    "S, Z = s.get_SZ()\n",
    "S\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cross-sections \u2014 T-matrix vs Mie\n"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "size_param = 2.0 * np.pi * radius_mm / wavelength_mm\n",
    "sigma_sca_tm = scatter.sca_xsect(s, h_pol=True)\n",
    "sigma_ext_tm = scatter.ext_xsect(s, h_pol=True)\n",
    "\n",
    "geom = np.pi * radius_mm ** 2\n",
    "sigma_sca_mie = mie_qsca(size_param, m.real, m.imag) * geom\n",
    "sigma_ext_mie = mie_qext(size_param, m.real, m.imag) * geom\n",
    "\n",
    "print(f'rel err \u03c3_sca = {abs(sigma_sca_tm-sigma_sca_mie)/sigma_sca_mie:.2e}')\n",
    "print(f'rel err \u03c3_ext = {abs(sigma_ext_tm-sigma_ext_mie)/sigma_ext_mie:.2e}')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Q_sca and Q_ext across a size-parameter sweep\n",
    "\n",
    "Plotting both curves side by side makes the Mie ripple pattern\n",
    "visible and confirms that the T-matrix follows it point-for-point.\n"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "xs = np.linspace(0.1, 10.0, 40)\n",
    "q_sca_mie = np.array([mie_qsca(x, m.real, m.imag) for x in xs])\n",
    "q_ext_mie = np.array([mie_qext(x, m.real, m.imag) for x in xs])\n",
    "q_sca_tm = np.empty_like(xs)\n",
    "q_ext_tm = np.empty_like(xs)\n",
    "for i, x in enumerate(xs):\n",
    "    r = x * wavelength_mm / (2.0 * np.pi)\n",
    "    si = Scatterer(radius=r, wavelength=wavelength_mm, m=m,\n",
    "                   axis_ratio=1.0, ddelt=1e-4, ndgs=2)\n",
    "    si.set_geometry(geom_horiz_back)\n",
    "    g = np.pi * r ** 2\n",
    "    q_sca_tm[i] = scatter.sca_xsect(si) / g\n",
    "    q_ext_tm[i] = scatter.ext_xsect(si) / g\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(7, 4))\n",
    "ax.plot(xs, q_sca_mie, 'k-', label='Q_sca (Mie)')\n",
    "ax.plot(xs, q_ext_mie, 'k--', label='Q_ext (Mie)')\n",
    "ax.plot(xs, q_sca_tm, 'C1o', markersize=4, label='Q_sca (T-matrix)')\n",
    "ax.plot(xs, q_ext_tm, 'C2s', markersize=4, label='Q_ext (T-matrix)')\n",
    "ax.set_xlabel('size parameter  x = 2\u03c0r/\u03bb')\n",
    "ax.set_ylabel('efficiency')\n",
    "ax.legend()\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
}