{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0885676e-8187-4eb8-8558-1a95aeaf52d0",
   "metadata": {},
   "source": [
    "# 4. Scattering calculations\n",
    "\n",
    "This example demonstrates the basic principles of performing scattering calcualtions with ARTS. Populations of particles that scatter radiation are referred to as *scattering species*. Each scattering species defines a mapping from one or several fields of *scattering species properties* to *scattering properties* which form the input to the actual scattering calculation."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ade3f219-304f-462e-8ea6-d7e1c7587834",
   "metadata": {},
   "source": [
    "## Set up a 1D atmosphere"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "05fd23f9-9c20-4df4-b8ba-09296faa6944",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pyarts\n",
    "\n",
    "ws = pyarts.Workspace()\n",
    "ws.surface_fieldPlanet(option=\"Earth\")\n",
    "ws.surface_field[pyarts.arts.SurfaceKey(\"t\")] = 295.0\n",
    "ws.atmospheric_fieldRead(\n",
    "    toa=100e3, basename=\"planets/Earth/afgl/tropical/\", missing_is_zero=1\n",
    ")\n",
    "ws.atmospheric_fieldIGRF(time=\"2000-03-11 14:39:37\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8e88bc4b-1bf7-4575-a99a-481c4b93cbb2",
   "metadata": {},
   "source": [
    "## Add a field of scattering species properties\n",
    "\n",
    "For this example, we will use scatterers with a Henyey-Greenstein phase function to represent ice particles. We will use two *scattering species properties* to represent the scattering species *ice*: The extinction and the single-scattering albed (SSA). The ice extinction and SSA thus become atmospheric fields. Below we define scattering-species-property object that identify these fields."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fde7da44-6123-44b6-a545-66e6271a892d",
   "metadata": {},
   "source": [
    "### Ice"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "56b62cec-4f34-4e3d-9158-dd9b6973548d",
   "metadata": {},
   "outputs": [],
   "source": [
    "ice_extinction = pyarts.arts.ScatteringSpeciesProperty(\"ice\", pyarts.arts.ParticulateProperty(\"Extinction\"))\n",
    "ice_ssa = pyarts.arts.ScatteringSpeciesProperty(\"ice\", pyarts.arts.ParticulateProperty(\"SingleScatteringAlbedo\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "59bfd385-6125-4751-94da-c06f85a0f599",
   "metadata": {},
   "source": [
    "We then define a GriddedField3 representing the ice extinction and single-scattering albedo and add it to ``atm_field`` of the workspace."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f8eddf1c-5288-408d-9b4a-7047c0f69b23",
   "metadata": {},
   "outputs": [],
   "source": [
    "grids = ws.atmospheric_field[\"t\"].data.grids\n",
    "z = grids[0]\n",
    "ice_extinction_profile = np.zeros_like(z)\n",
    "ice_extinction_profile[(z > 5e3) * (z < 15e3)] = 1.0\n",
    "ice_extinction_profile = pyarts.arts.GriddedField3(data=ice_extinction_profile[..., None, None], grids=grids)\n",
    "ws.atmospheric_field[ice_extinction] = ice_extinction_profile\n",
    "\n",
    "ice_ssa_profile = np.zeros_like(z)\n",
    "ice_ssa_profile[(z > 5e3) * (z < 15e3)] = 0.5\n",
    "ice_ssa_profile = pyarts.arts.GriddedField3(data=ice_ssa_profile[..., None, None], grids=grids)\n",
    "ws.atmospheric_field[ice_ssa] = ice_ssa_profile"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9fa4cdf8-3e41-4120-a9db-c95f3973e4e5",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "1523f263-1c70-4c30-b6d4-8c5f8a2ae68b",
   "metadata": {},
   "source": [
    "### Rain"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a94b899d-f78e-4772-b46f-53378b17f276",
   "metadata": {},
   "outputs": [],
   "source": [
    "rain_extinction = pyarts.arts.ScatteringSpeciesProperty(\"rain\", pyarts.arts.ParticulateProperty(\"Extinction\"))\n",
    "rain_ssa = pyarts.arts.ScatteringSpeciesProperty(\"rain\", pyarts.arts.ParticulateProperty(\"SingleScatteringAlbedo\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bee6b305-76f3-4ecc-9a3f-bc7642384c8d",
   "metadata": {},
   "outputs": [],
   "source": [
    "rain_extinction_profile = np.zeros_like(z)\n",
    "rain_extinction_profile[z < 5e3] = 1.0\n",
    "rain_extinction_profile = pyarts.arts.GriddedField3(data=rain_extinction_profile[..., None, None], grids=grids)\n",
    "ws.atmospheric_field[rain_extinction] = rain_extinction_profile\n",
    "\n",
    "rain_ssa_profile = np.zeros_like(z)\n",
    "rain_ssa_profile[z < 5e3] = 0.5\n",
    "rain_ssa_profile = pyarts.arts.GriddedField3(data=rain_ssa_profile[..., None, None], grids=grids)\n",
    "ws.atmospheric_field[rain_ssa] = rain_ssa_profile"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "056c6b3c-7cab-4d61-bb24-0764ed1e86cc",
   "metadata": {},
   "source": [
    "## Create the scattering species"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "74c5c106-4c66-44ea-a837-79a30e1424dc",
   "metadata": {},
   "outputs": [],
   "source": [
    "from pyarts.arts import HenyeyGreensteinScatterer, ArrayOfScatteringSpecies\n",
    "\n",
    "hg_ice = HenyeyGreensteinScatterer(ice_extinction, ice_ssa, 0.5)\n",
    "hg_rain = HenyeyGreensteinScatterer(rain_extinction, rain_ssa, 0.0)\n",
    "scattering_species = ArrayOfScatteringSpecies()\n",
    "scattering_species.add(hg_ice)\n",
    "scattering_species.add(hg_rain)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d7342431-699a-40da-a183-c5ca3e03a6a1",
   "metadata": {},
   "source": [
    "## Extracting scattering properties\n",
    "\n",
    "We can now extract bulk scattering properties from the scattering species."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3093e2be-366b-4cd4-9a35-77392a285bcb",
   "metadata": {},
   "source": [
    "### Ice\n",
    "\n",
    "We can extract the ice phase function from the combind scattering species by calculating the bulk-scattering properties at an altitude above 5 km. To check the consistency of the Henyey-Greenstein scatterer, we extract the data in gridded and spectral representation and ensure that they are the same when both are converted to gridded representation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "43d853f8-62f7-4a2f-a740-9a52cbf39bb3",
   "metadata": {},
   "outputs": [],
   "source": [
    "atm_pt = ws.atmospheric_field(11e3, 0.0, 0.0)\n",
    "f_grid = np.array([89e9])\n",
    "bsp = scattering_species.get_bulk_scattering_properties_tro_spectral(atm_pt, f_grid, 32)\n",
    "pm_spectral = bsp.phase_matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6710e1fc-28e1-43e9-b074-8cbbef6b2b79",
   "metadata": {},
   "outputs": [],
   "source": [
    "za_scat_grid = pyarts.arts.ScatteringTroSpectralVector.to_gridded(pm_spectral).get_za_scat_grid()\n",
    "bsp = scattering_species.get_bulk_scattering_properties_tro_gridded(atm_pt, f_grid, za_scat_grid)\n",
    "pm_gridded = bsp.phase_matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e3da481e-1ce8-43ae-96f6-f20192748afb",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "plt.plot(pyarts.arts.ScatteringTroSpectralVector.to_gridded(pm_spectral).flatten()[::6], label=\"Converted from Spectral\")\n",
    "plt.plot(pm_gridded.flatten()[::6], ls=\"--\", label=\"Gridded\")\n",
    "plt.title(\"Henyey-Greenstein phase function\")\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e389b4b1-c997-4e50-9d46-d695753e2d1f",
   "metadata": {},
   "source": [
    "### Rain\n",
    "\n",
    "Similarly, we can extract the phase function for rain by calculating the bulk-scattering properties at a lower point in the atmosphere."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0ac9d751-328e-4b00-b611-e56cef020978",
   "metadata": {},
   "outputs": [],
   "source": [
    "atm_pt = ws.atmospheric_field(4e3, 0.0, 0.0)\n",
    "bsp = scattering_species.get_bulk_scattering_properties_tro_spectral(atm_pt, f_grid, 32)\n",
    "pm_spectral = bsp.phase_matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b345c993-a997-4925-a3d7-2ad5c1636997",
   "metadata": {},
   "outputs": [],
   "source": [
    "za_scat_grid = pyarts.arts.ScatteringTroSpectralVector.to_gridded(pm_spectral).get_za_scat_grid()\n",
    "bsp = scattering_species.get_bulk_scattering_properties_tro_gridded(atm_pt, f_grid, za_scat_grid)\n",
    "pm_gridded = bsp.phase_matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "04a485a5-cd07-4b46-abf0-7db349073c8f",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "plt.plot(pyarts.arts.ScatteringTroSpectralVector.to_gridded(pm_spectral).flatten()[::6], label=\"Converted from Spectral\")\n",
    "plt.plot(pm_gridded.flatten()[::6], ls=\"--\", label=\"Gridded\")\n",
    "plt.title(\"Henyey-Greenstein phase function\")\n",
    "plt.legend()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}