{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Using the Baum Ice Crystal Database with HR"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import sasktran as sk"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {
    "raw_mimetype": "text/restructuredtext"
   },
   "source": [
    "SASKTRAN supports calculations using the Baum ice crystal database.  The core of this functionality is contained within the :py:class:`BaumIceCrystal <sasktran.BaumIceCrystal>` optical property.  In principle including ice crystals is as simple as creating a climatology to go along with the optical property and adding it to the engine, however complications involving the radiative transfer of ice clouds usually means that the HR engine must be configured in a specific way to obtain accurate calculations.\n",
    "\n",
    "The :py:class:`SpeciesBaumIceCloud <sasktran.SpeciesBaumIceCloud>` species combines the :py:class:`BaumIceCrystal <sasktran.BaumIceCrystal>` optical property with a climatology that creates a Gaussian shaped cloud with a specified size and optical thickness."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "baum_species = sk.SpeciesBaumIceCloud(particlesize_microns=20, \n",
    "                                      cloud_top_m=16500,\n",
    "                                      cloud_width_fwhm_m=500,\n",
    "                                      vertical_optical_depth=0.01,\n",
    "                                      vertical_optical_depth_wavel_nm=750)"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {
    "raw_mimetype": "text/restructuredtext"
   },
   "source": [
    "Which creates a cloud species that consists of 20 micron sized particles, a cloud top altitude of 16.5 km, a full width at half maximum in altitude of 500 m, and has a vertical optical depth of 0.01 at 750 nm.  Note that limb path lengths can be ~100 times that of vertical path lengths, so even a vertically thin cloud may appear thick in the limb geometry.  Shown below is the calculated Gaussian number density."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "altitudes = np.arange(10000, 20000, 10)\n",
    "\n",
    "cloud_dens = baum_species.climatology.get_parameter('icecloud', latitude=0, longitude=0, mjd=54372, altitudes=altitudes)\n",
    "\n",
    "plt.plot(cloud_dens, altitudes)\n",
    "plt.ylabel('Altitude [m]')\n",
    "plt.xlabel('Number Density [/cm3]')"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {
    "raw_mimetype": "text/restructuredtext"
   },
   "source": [
    "We can now set up the HR engine and add the species to it as we would any other species.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tanalts_km = np.arange(10, 50, 1)\n",
    "\n",
    "# First recreate our geometry and atmosphere classes\n",
    "geometry = sk.VerticalImage()\n",
    "geometry.from_sza_saa(sza=60, saa=60, lat=0, lon=0, tanalts_km=tanalts_km, mjd=54372, locallook=0,\n",
    "                      satalt_km=600, refalt_km=20)\n",
    "\n",
    "atmosphere = sk.Atmosphere()\n",
    "\n",
    "atmosphere['ozone'] = sk.Species(sk.O3OSIRISRes(), sk.Labow())\n",
    "atmosphere['air'] = sk.Species(sk.Rayleigh(), sk.MSIS90())\n",
    "atmosphere['cloud'] = baum_species\n",
    "\n",
    "# And now make the engine\n",
    "engine = sk.EngineHR(geometry=geometry, atmosphere=atmosphere)\n",
    "\n",
    "# Choose some wavelengths to do the calculation at\n",
    "engine.wavelengths = [340, 800]"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {
    "raw_mimetype": "text/restructuredtext"
   },
   "source": [
    "But before we run the calculation we should modify some settings in the engine to account for the cloud accurately."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setting the model parameters"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {
    "raw_mimetype": "text/restructuredtext"
   },
   "source": [
    "The default vertical resolution of SASKTRAN-HR is ~1000 m, however clouds typically have fine structures and the cloud included above is calculated at a 10 m resolution by default.  We could set the grid spacing of HR to be 10 m, however this would result in a very slow calculation.  It also is not necessary because the cloud is localized over a few km, we do not need 10 m resolution throughout the entire atmosphere.  There are also several settings related to the adaptive integration of rays that should be used when including clouds.  For convenience the :py:meth:`configure_for_cloud <sasktran.EngineHR.configure_for_cloud>` method will set most of these things."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "engine.configure_for_cloud(cloud_altitudes=baum_species.altitude_grid, grid_spacing_m=1000)\n",
    "\n",
    "rad = engine.calculate_radiance()\n",
    "plt.plot(rad.T, tanalts_km)\n",
    "plt.ylabel('Altitude [km]')\n",
    "plt.xlabel('Radiance')"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Raw Cell Format",
  "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.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}