{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0925ec42",
   "metadata": {},
   "source": [
    "# Computing scattering and extinction efficiencies\n",
    "\n",
    "Different physical assumptions and numerical methods can be used to model how light propagates through a solid material. The `xdust.scatteringmodel` module provides several interchangeable scattering models:\n",
    "\n",
    "- `xdust.scatteringmodel.Mie()` -- the Bohren & Huffman (1983) Mie-scattering algorithm, sped up with vectorized computations. This can be demanding on RAM depending on the wavelength resolution and number of grain radii used, so some care is needed at high resolution.\n",
    "- `xdust.scatteringmodel.RGscattering()` -- the Rayleigh-Gans approximation, relevant for grains much larger than the wavelength of light and relatively transparent to it (i.e. $|m - 1| \\ll 1$, where $m$ is the complex index of refraction).\n",
    "- `xdust.scatteringmodel.PAH` -- reads and interpolates tables of extinction properties for polycyclic aromatic hydrocarbons (PAHs) from Li & Draine (2001).\n",
    "\n",
    "By the end of this tutorial, you will be able to:\n",
    "\n",
    "- run a scattering calculation with `RGscattering` and `Mie`\n",
    "- understand the shape and meaning of the resulting efficiency arrays (`qsca`, `qabs`, `qext`, `diff`)\n",
    "- compute and verify a differential scattering cross-section as a function of angle\n",
    "- avoid exceeding available memory when using the `Mie` model\n",
    "\n",
    "To compute extinction, a scattering model needs a `xdust.graindist.composition.Composition` object (which supplies the optical constants of the compound) and an array of grain radii. If the grain radii have no attached units, they are assumed to be in microns.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "466c8134",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import astropy.units as u\n",
    "import astropy.constants as c\n",
    "\n",
    "import xdust"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "af05b09f",
   "metadata": {},
   "source": [
    "## The Rayleigh-Gans model\n",
    "\n",
    "This example computes Rayleigh-Gans scattering using the Drude approximation for the index of refraction.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "75fa33dd",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Energy grid for the calculation\n",
    "energy_grid = np.logspace(-1, 1, 30) * u.keV\n",
    "\n",
    "# A single grain size\n",
    "grain_radius = np.array([1.0]) * u.micron\n",
    "\n",
    "# Angular grid, used later for the differential cross-section\n",
    "theta_grid = np.logspace(-5., np.log10(np.pi), 1000) * u.rad"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "55d3c9f4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# A Composition object using the Drude approximation\n",
    "drude = xdust.graindist.composition.CmDrude()\n",
    "\n",
    "# The Rayleigh-Gans scattering model\n",
    "rg_model = xdust.scatteringmodel.RGscattering()\n",
    "\n",
    "rg_model.calculate(energy_grid, grain_radius, drude)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "43a7fa48",
   "metadata": {},
   "source": [
    "The results of the calculation -- various attenuation efficiencies -- are stored as attributes on the `ScatteringModel` object. Efficiency is the physical cross-section divided by the geometric cross-section; for a spherical grain,\n",
    "\n",
    "$$ Q = \\frac{\\sigma}{\\pi a^2} $$\n",
    "\n",
    "where $\\sigma$ is the cross-section for the physical interaction (e.g. scattering) and $a$ is the grain radius. Efficiencies are unitless.\n",
    "\n",
    "- `ScatteringModel.qsca` -- scattering efficiency\n",
    "- `ScatteringModel.qabs` -- absorption efficiency\n",
    "- `ScatteringModel.qext` -- extinction efficiency (`qsca` + `qabs`)\n",
    "- `ScatteringModel.diff` -- differential scattering cross-section, divided by the geometric cross-section\n",
    "\n",
    "`qsca`, `qabs`, and `qext` are 2-D arrays with shape `(NE, NA)`, where `NE` is the length of the input energy/wavelength grid and `NA` is the length of the input grain-size grid.\n",
    "\n",
    "`diff` is a 3-D array with shape `(NE, NA, NTH)`, where `NTH` is the length of the angular grid passed via the `theta` keyword of `ScatteringModel.calculate`.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b5410a4e",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(energy_grid, rg_model.qsca)\n",
    "plt.loglog()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d4179e45",
   "metadata": {},
   "source": [
    "The Rayleigh-Gans-plus-Drude approximation gives a scattering cross-section that decays as $E^{-2}$.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ba0d0daf",
   "metadata": {},
   "source": [
    "## Computing a differential scattering cross-section\n",
    "\n",
    "To compute the differential cross-section as a function of angle, pass an array of angles via the `theta` keyword (with or without units -- unitless values are assumed to be radians). By default `theta` is `0`; you only need to supply it when you want scattering as a function of angle.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4edc622b",
   "metadata": {},
   "outputs": [],
   "source": [
    "rg_with_angles = xdust.scatteringmodel.RGscattering()\n",
    "rg_with_angles.calculate(energy_grid, grain_radius, drude, theta=theta_grid)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a53bb52d",
   "metadata": {},
   "source": [
    "It helps to check the shape of the resulting arrays before working with them.\n",
    "\n",
    "`qsca` is, by definition, integrated over all scattering angles, so it depends only on energy/wavelength and grain radius -- its shape is `(len(energy_grid), len(grain_radius))`.\n",
    "\n",
    "`diff` also depends on angle, so its shape is `(len(energy_grid), len(grain_radius), len(theta_grid))`.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4fe23ef4",
   "metadata": {},
   "outputs": [],
   "source": [
    "rg_with_angles.qsca.shape\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "19e0e278",
   "metadata": {},
   "outputs": [],
   "source": [
    "rg_with_angles.diff.shape\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6e78a818",
   "metadata": {},
   "source": [
    "With only one grain radius in the grid, we always index the second dimension with `0`. Below, we compare the differential scattering cross-section at the lowest energy (index `0`) and the highest energy (index `-1`, the last value in `energy_grid`).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7f8eee91",
   "metadata": {},
   "outputs": [],
   "source": [
    "angle_unit = 'arcmin'\n",
    "plt.plot(theta_grid.to(angle_unit), rg_with_angles.diff[0, 0, :], 'b-', lw=2, label=energy_grid[0])\n",
    "plt.plot(theta_grid.to(angle_unit), rg_with_angles.diff[-1, 0, :], 'k--', lw=2, label=energy_grid[-1])\n",
    "plt.title(\"Grain radius: {}\".format(grain_radius[0]))\n",
    "plt.xlabel(r'$\\theta$ ({})'.format(angle_unit), size=12)\n",
    "plt.ylabel(r'$dQ_{\\rm sca}/d\\Omega$ (ster$^{-1}$)', size=12)\n",
    "plt.loglog()\n",
    "plt.legend()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dee82f21",
   "metadata": {},
   "source": [
    "As a check, the differential cross-section should integrate over solid angle to the same scattering efficiency computed above. This is generally accurate to within a few percent.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c434ab1f",
   "metadata": {},
   "outputs": [],
   "source": [
    "i = -1  # index of the energy value to check\n",
    "integrand = np.trapezoid(\n",
    "    rg_with_angles.diff[i, 0, :] * 2.0 * np.pi * np.sin(theta_grid.to('rad').value),\n",
    "    theta_grid.to('rad').value,\n",
    ")\n",
    "\n",
    "print(integrand / rg_with_angles.qsca[i, 0])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4be70b26",
   "metadata": {},
   "source": [
    "## The Mie model\n",
    "\n",
    "This example computes a Mie scattering model for silicate dust grains at a single wavelength.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e1344f22",
   "metadata": {},
   "outputs": [],
   "source": [
    "# A visible wavelength\n",
    "wavelength_V = 4500. * u.angstrom\n",
    "\n",
    "# A silicate composition\n",
    "silicate = xdust.graindist.composition.CmSilicate()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "62c26e31",
   "metadata": {},
   "outputs": [],
   "source": [
    "mie_model = xdust.scatteringmodel.Mie()\n",
    "mie_model.calculate(wavelength_V, grain_radius, silicate, theta=theta_grid)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a2bd82dc",
   "metadata": {},
   "source": [
    "With a single wavelength value, the shape of `diff` is worth checking again.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8c13870f",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.shape(mie_model.diff)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e29443ad",
   "metadata": {},
   "outputs": [],
   "source": [
    "angle_unit = 'rad'\n",
    "plt.plot(theta_grid.to(angle_unit), mie_model.diff[0, 0, :])\n",
    "plt.semilogy()\n",
    "plt.xlabel(r'$\\theta$ ({})'.format(angle_unit), size=12)\n",
    "plt.ylabel(r'$dQ_{\\rm sca}/d\\Omega$ (ster$^{-1}$)', size=12)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "139dd86b",
   "metadata": {},
   "source": [
    "As before, verify that the differential cross-section integrates to the scattering efficiency.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "771490f4",
   "metadata": {},
   "outputs": [],
   "source": [
    "check = np.trapezoid(\n",
    "    mie_model.diff * 2.0 * np.pi * np.sin(theta_grid.to('rad').value),\n",
    "    theta_grid.to('rad').value,\n",
    ")\n",
    "print(check / mie_model.qsca)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d5071f4d",
   "metadata": {},
   "source": [
    "### Note: the Mie model can use a lot of memory\n",
    "\n",
    "The `Mie` model uses multi-dimensional array operations instead of for-loops, which is fast but can easily exceed the RAM available on a typical laptop or desktop and crash your system. To prevent this, `Mie.calculate` estimates the memory the computation would require *before* running it, and prints a warning instead of performing the calculation if that estimate exceeds a limit (8 GB by default, adjustable with the `memlim` keyword).\n",
    "\n",
    "The cell below is a deliberately unrealistic example -- an energy grid reaching 5 GeV, far outside any regime where dust scattering is physically relevant -- included only to demonstrate what the warning looks like. **This is expected behavior, not an error:** the calculation is intentionally skipped, and the resulting `qsca`, `qext`, `qabs`, `gsca`, and `diff` attributes are all set to `0.0` rather than the real arrays.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cc2db165",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Deliberately oversized grid, chosen only to trigger the memory-limit warning\n",
    "oversized_energy_grid = np.linspace(1000., 5000., 2) * u.keV\n",
    "oversized_radius_grid = np.linspace(0.1, 0.5, 20) * u.micron\n",
    "\n",
    "mie_oversized = xdust.scatteringmodel.Mie()\n",
    "mie_oversized.calculate(oversized_energy_grid, oversized_radius_grid, silicate, theta=theta_grid)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6e22a59f",
   "metadata": {},
   "source": [
    "If you have more RAM available and want to proceed anyway, raise `memlim` (in GB) to a value above the printed estimate:\n",
    "\n",
    "```python\n",
    "mie_oversized = xdust.scatteringmodel.Mie()\n",
    "mie_oversized.calculate(oversized_energy_grid, oversized_radius_grid, silicate, theta=theta_grid, memlim=8.3)\n",
    "```\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "72b0d95d",
   "metadata": {},
   "source": [
    "## A Mie scattering calculation with all three dimensions\n",
    "\n",
    "The example below uses a moderate grid so that energy/wavelength, grain radius, and scattering angle are all resolved with more than one value, which makes it easier to see how `diff` varies along each dimension.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "45426b6a",
   "metadata": {},
   "outputs": [],
   "source": [
    "wavelength_grid = np.linspace(1000., 5000., 5) * u.angstrom\n",
    "radius_grid = np.linspace(0.1, 0.5, 20) * u.micron\n",
    "\n",
    "mie_grid = xdust.scatteringmodel.Mie()\n",
    "mie_grid.calculate(wavelength_grid, radius_grid, silicate, theta=theta_grid)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a0543a84",
   "metadata": {},
   "source": [
    "Compare the differential scattering cross-section at two different wavelengths (same grain radius):\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a8da46af",
   "metadata": {},
   "outputs": [],
   "source": [
    "angle_unit = 'rad'\n",
    "plt.plot(theta_grid.to(angle_unit), mie_grid.diff[0, 0, :], label=wavelength_grid[0])\n",
    "plt.plot(theta_grid.to(angle_unit), mie_grid.diff[-1, 0, :], label=wavelength_grid[-1])\n",
    "plt.semilogy()\n",
    "\n",
    "plt.title(\"Grain radius: {}\".format(radius_grid[0]), size=12)\n",
    "plt.xlabel(r'$\\theta$ ({})'.format(angle_unit), size=12)\n",
    "plt.ylabel(r'$dQ_{\\rm sca}/d\\Omega$ (ster$^{-1}$)', size=12)\n",
    "plt.legend()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "465b1efa",
   "metadata": {},
   "source": [
    "...or at two different grain sizes (same wavelength):\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b5fbecfa",
   "metadata": {},
   "outputs": [],
   "source": [
    "angle_unit = 'rad'\n",
    "plt.plot(theta_grid.to(angle_unit), mie_grid.diff[0, 0, :], label=radius_grid[0])\n",
    "plt.plot(theta_grid.to(angle_unit), mie_grid.diff[0, -1, :], label=radius_grid[-1])\n",
    "plt.semilogy()\n",
    "\n",
    "plt.title(\"Wavelength: {}\".format(wavelength_grid[0]), size=12)\n",
    "plt.xlabel(r'$\\theta$ ({})'.format(angle_unit), size=12)\n",
    "plt.ylabel(r'$dQ_{\\rm sca}/d\\Omega$ (ster$^{-1}$)', size=12)\n",
    "plt.legend()\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python",
   "pygments_lexer": "ipython3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}