{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Approximations for Toon89 Two-Stream Radiative Transfer in Reflected Light\n",
    "\n",
    "Like any code, there are several approximations that go into computing intesity from various phase functions. In reflected light models, some of these approximations drastically change the output spectra. \n",
    "\n",
    "In this notebook you will: \n",
    "\n",
    "- learn how to use the `approx` method to access different ways of computing radiative transfer \n",
    "- focusing on Toon et al 1989 methodology, learn how each approximation affects reflected light spectrum \n",
    "- how to run run a benchmark test against reference data within the PICASO code (leverages data from [Dlugach & Yanovitskij (1974)](https://ui.adsabs.harvard.edu/abs/1974Icar...22...66D/abstract) to replicates the study of [Batalha et al. (2019)](https://ui.adsabs.harvard.edu/abs/2019ApJ...878...70B/abstract) )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "#picaso\n",
    "from picaso import justdoit as jdi \n",
    "from picaso import justplotit as jpi \n",
    "from bokeh.plotting import show, figure\n",
    "from bokeh.layouts import column\n",
    "from bokeh.palettes import Colorblind8\n",
    "from bokeh.io import output_notebook \n",
    "import astropy.units as u\n",
    "output_notebook()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Using the `approx` key"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "opa = jdi.opannection(wave_range=[0.3,1])\n",
    "cloud_free = jdi.inputs()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that all the `approx` keys have predefined inputs. These are our recommendations for how to run the code. But, users should always be weary of these and test their sensitivity to your results."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## How to see what radiative transfer scattering options exist?\n",
    "\n",
    "All of the approximations have to do with how scattering is handled. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('Options for Direct Scattring Phase: ', jdi.single_phase_options())\n",
    "\n",
    "print('Options for Multiple Scattring Phase: ', jdi.multi_phase_options())\n",
    "\n",
    "print('Options for Raman Scattering: ', jdi.raman_options())\n",
    "\n",
    "print('Options for Toon89 Coefficients: ', jdi.toon_phase_coefficients())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Set inputs normally"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cloud_free.phase_angle(0) #phase in radians\n",
    "\n",
    "cloud_free.gravity(gravity=25, gravity_unit=u.Unit('m/(s**2)' ))\n",
    "\n",
    "#set star \n",
    "cloud_free.star(opa, 6000, 0.0122, 4.437)\n",
    "\n",
    "#set atmosphere comp and temp \n",
    "cloud_free.atmosphere(filename=jdi.jupiter_pt(),delim_whitespace=True)\n",
    "\n",
    "#make a copy to have a separate cloud input dict\n",
    "from copy import deepcopy\n",
    "cloudy=deepcopy(cloud_free)\n",
    "cloudy.clouds( filename=jdi.jupiter_cld(),delim_whitespace=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Toon Coefficients \n",
    "\n",
    "Though the default is Quadrature, Table 1 of [Toon et al. 1989]() gives two options for scattering phase functions for solar, reflected light."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#let's make two different figures for this\n",
    "fig_cloudy = figure(x_range=[0.3,1],y_range=[0,1.2], width=500, height=300\n",
    "                    ,x_axis_label='Wavelength [μm]', y_axis_label='Geometric Albedo',)\n",
    "fig_no_cloud = figure(x_range=[0.3,1],y_range=[0,1.2], width=500, height=300\n",
    "                      ,x_axis_label='Wavelength [μm]', y_axis_label='Geometric Albedo',)\n",
    "\n",
    "#define our options\n",
    "options = jdi.toon_phase_coefficients()\n",
    "colors = Colorblind8[0:len(options)]\n",
    "\n",
    "#loop through all approximations \n",
    "for approx, c in zip(options, colors):\n",
    "    #set approximations\n",
    "    cloud_free.approx(toon_coefficients = approx)    \n",
    "    cloudy.approx(toon_coefficients = approx)  \n",
    "    df = cloud_free.spectrum(opa)\n",
    "    wno_nc, alb_nc = jdi.mean_regrid(df['wavenumber'] , df['albedo'],R=150)\n",
    "    df = cloudy.spectrum(opa)\n",
    "    wno_c, alb_c = jdi.mean_regrid(df['wavenumber'] , df['albedo'],R=150)\n",
    "    fig_no_cloud.line(1e4/wno_nc, alb_nc, legend_label=approx, color=c, line_width=3)\n",
    "    fig_cloudy.line(1e4/wno_c, alb_c,  color=c, line_width=3)\n",
    "jpi.plot_format(fig_cloudy)\n",
    "jpi.plot_format(fig_no_cloud)\n",
    "show(column(fig_no_cloud,fig_cloudy ))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Direct Scattering Approximation\n",
    "\n",
    "The [derivation documentation](https://natashabatalha.github.io/picaso_dev#slide02) has a full description of these direct scattering approximations. Briefly I'll describe them here: \n",
    "\n",
    "At the center of each is the [One Term Henyey-Greenstein Phase Function (OTHG)](http://adsabs.harvard.edu/abs/1941ApJ....93...70H) and the [Two Term HG Phase Function (TTHG)](http://adsabs.harvard.edu/abs/1965ApJ...142.1563I). \n",
    "\n",
    "We also know that planet atmospheres have high degrees of Rayleigh scattering. [Cahoy+2010](http://adsabs.harvard.edu/abs/2010ApJ...724..189C) developed a methodology for incorporating Rayleigh into the direct scattering component. \n",
    "\n",
    "A more rrobust way of dealing with Rayleigh is to directly fold it's phase function into the TTHG phase function (TTHG_Ray). \n",
    "\n",
    "We'll run each case with and without a cloud so you can see what happens in both regimes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#let's make two different figures for this\n",
    "fig_cloudy = figure(x_range=[0.3,1],y_range=[0,1.2], width=500, height=300\n",
    "                    ,x_axis_label='Wavelength [μm]', y_axis_label='Geometric Albedo',)\n",
    "fig_no_cloud = figure(x_range=[0.3,1],y_range=[0,1.2], width=500, height=300\n",
    "                      ,x_axis_label='Wavelength [μm]', y_axis_label='Geometric Albedo',)\n",
    "\n",
    "#define our options\n",
    "options = jdi.single_phase_options()\n",
    "colors = Colorblind8[0:len(options)]\n",
    "\n",
    "#loop through all approximations \n",
    "for approx, c in zip(options, colors):\n",
    "    #set approximations\n",
    "    cloud_free.approx(single_phase = approx,raman='pollack')\n",
    "    cloudy.approx(single_phase = approx,raman='pollack')\n",
    "    df = cloud_free.spectrum(opa)\n",
    "    wno_nc, alb_nc = jdi.mean_regrid(df['wavenumber'] , df['albedo'],R=150)\n",
    "    df = cloudy.spectrum(opa)\n",
    "    wno_c, alb_c = jdi.mean_regrid(df['wavenumber'] , df['albedo'],R=150)\n",
    "    fig_no_cloud.line(1e4/wno_nc, alb_nc, legend_label=approx, color=c, line_width=3)\n",
    "    fig_cloudy.line(1e4/wno_c, alb_c,  color=c, line_width=3)\n",
    "jpi.plot_format(fig_cloudy)\n",
    "jpi.plot_format(fig_no_cloud)\n",
    "show(column(fig_no_cloud,fig_cloudy ))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Multiple Scattering Approximations\n",
    "\n",
    "Again, [derivation documentation](https://natashabatalha.github.io/picaso_dev#slide03) has a full description of these multiple scattering approximations. \n",
    "\n",
    "To complete the multiple scattering integration over all _diffuse angles_, we have to use some mathematical tricks. [Legendre Polynomials](http://mathworld.wolfram.com/LegendrePolynomial.html) are often used to complete this integration to varying degrees. For Solar System/Exoplanet papers, we often stop the expansion at either `N=1` or `N=2`. Our standard will be to run with `N=2`, but below we show how to run each."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = figure(x_range=[0.3,1],y_range=[0,1.2], width=500, height=300,\n",
    "             x_axis_label='Wavelength [μm]', y_axis_label='Geometric Albedo',)\n",
    "\n",
    "options = jdi.multi_phase_options()\n",
    "colors = Colorblind8[0:len(options)*2]\n",
    "\n",
    "for approx, c1,c2 in zip(options, colors[0:2], colors[2:]):\n",
    "    cloud_free.approx(multi_phase= approx)\n",
    "    cloudy.approx(multi_phase = approx)\n",
    "    df = cloud_free.spectrum(opa)\n",
    "    wno_nc, alb_nc =jdi.mean_regrid( df['wavenumber'] , df['albedo'],R=150)\n",
    "    df = cloudy.spectrum(opa)\n",
    "    wno_c, alb_c =jdi.mean_regrid( df['wavenumber'] , df['albedo'],R=150)\n",
    "    fig.line(1e4/wno_nc, alb_nc, color=c1, line_width=3)\n",
    "    fig.line(1e4/wno_c, alb_c,  color=c2, line_width=3)\n",
    "jpi.plot_format(fig)\n",
    "show(fig)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Raman Scattering Approximations\n",
    "\n",
    "We all know the importance of Rayleigh scattering in planetary atmospheres. Raman scattering also has important implications for our spectra (these features have been observed in Solar System planets). In particular, at short wavelengths, Raman scattering imprints molecular features from the star on the planetary spectrum. \n",
    "\n",
    "The most complete analysis of all Raman approximations is in [Sromosvky+2005](http://adsabs.harvard.edu/abs/2005Icar..173..254S). From these, we use the _Pollack Approximation_ that was used in [Cahoy+2010](http://adsabs.harvard.edu/abs/2010ApJ...724..189C) and others. \n",
    "\n",
    "We include the original Pollack methodology, but also include a modified version with [Oklopcic et al 2018](http://iopscience.iop.org/article/10.3847/0004-637X/832/1/30/meta) cross sections and updated methodology to include effects of stellar spectrum.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig_cloudy = figure(x_range=[0.3,1],y_range=[0,1.2], width=500, height=300\n",
    "                    ,x_axis_label='Wavelength [μm]', y_axis_label='Geometric Albedo',)\n",
    "fig_no_cloud = figure(x_range=[0.3,1],y_range=[0,1.2], width=500, height=300\n",
    "                      ,x_axis_label='Wavelength [μm]', y_axis_label='Geometric Albedo',)\n",
    "\n",
    "options = jdi.raman_options()\n",
    "colors = Colorblind8[0:len(options)]\n",
    "\n",
    "for approx, c in zip(options, colors):\n",
    "    cloud_free.approx(raman = approx)\n",
    "    cloud_free.star(opa, 6000, 0.0122, 4.437)\n",
    "    cloudy.approx(raman = approx)\n",
    "    cloudy.star(opa, 6000, 0.0122, 4.437)\n",
    "    df = cloud_free.spectrum(opa)\n",
    "    wno_nc, alb_nc = jdi.mean_regrid(df['wavenumber'] , df['albedo'],R=150)\n",
    "    df = cloudy.spectrum(opa)\n",
    "    wno_c, alb_c = jdi.mean_regrid(df['wavenumber'] , df['albedo'],R=150)\n",
    "    fig_no_cloud.line(1e4/wno_nc, alb_nc, legend_label=approx, color=c, line_width=3)\n",
    "    fig_cloudy.line(1e4/wno_c, alb_c,  color=c, line_width=3)\n",
    "\n",
    "jpi.plot_format(fig_cloudy)\n",
    "jpi.plot_format(fig_no_cloud)\n",
    "show(column(fig_no_cloud,fig_cloudy ))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Effect of Stellar Spectrum\n",
    "\n",
    "With the updated Raman scattering approximation, you will notice imprints of the stellar spectrum in the planet reflected light spectrum. Take a look below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig_cloudy = figure(x_range=[0.3,1],y_range=[0,1.2], width=500, height=300\n",
    "                    ,x_axis_label='Wavelength [μm]', y_axis_label='Geometric Albedo',)\n",
    "fig_no_cloud = figure(x_range=[0.3,1],y_range=[0,1.2], width=500, height=300\n",
    "                      ,x_axis_label='Wavelength [μm]', y_axis_label='Geometric Albedo',)\n",
    "\n",
    "#lets play around with different stellar spectra Teff\n",
    "stellar_teff = [6000,4000,3500]\n",
    "\n",
    "colors = Colorblind8[0:len(options)]\n",
    "\n",
    "cloud_free.approx(raman = 'oklopcic')\n",
    "cloudy.approx(raman = 'oklopcic')\n",
    "\n",
    "for approx, c in zip(stellar_teff, colors):\n",
    "    cloud_free.star(opa, approx, 0.0122, 4.437) \n",
    "    cloudy.star(opa, approx ,0.0122, 4.437) \n",
    "    df = cloud_free.spectrum(opa)\n",
    "    wno_nc, alb_nc = jdi.mean_regrid(df['wavenumber'] , df['albedo'],R=150)\n",
    "    df = cloudy.spectrum(opa)\n",
    "    wno_c, alb_c = jdi.mean_regrid(df['wavenumber'] , df['albedo'],R=150)\n",
    "    fig_no_cloud.line(1e4/wno_nc, alb_nc, legend_label=str(approx), color=c, line_width=3)\n",
    "    fig_cloudy.line(1e4/wno_c, alb_c,  color=c, line_width=3)\n",
    "\n",
    "cloud_free.approx(raman = 'pollack')\n",
    "cloudy.approx(raman = 'pollack')\n",
    "\n",
    "df = cloud_free.spectrum(opa)\n",
    "wno_nc, alb_nc = jdi.mean_regrid(df['wavenumber'] , df['albedo'],R=150)\n",
    "df = cloudy.spectrum(opa)\n",
    "wno_c, alb_c =jdi.mean_regrid( df['wavenumber'] , df['albedo'],R=150)\n",
    "fig_no_cloud.line(1e4/wno_nc, alb_nc, legend_label='Pollack', color='black', line_width=2, line_dash='dashed')\n",
    "fig_cloudy.line(1e4/wno_c, alb_c,  color='black', line_width=2, line_dash='dashed')\n",
    "\n",
    "\n",
    "jpi.plot_format(fig_cloudy)\n",
    "jpi.plot_format(fig_no_cloud)\n",
    "show(column(fig_no_cloud,fig_cloudy ))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Benchmark Toon89 w/ Dlugach & Yanovitskij (1974) \n",
    "\n",
    "This replicates Figure 9 from Batalha et al. 2019, which is a benchmark study with [Dlugach & Yanovitskij (1974)](https://ui.adsabs.harvard.edu/abs/1974Icar...22...66D/abstract)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import picaso.test as ptest"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Dlugach, Toon89= ptest.dlugach_test(delta_eddington=True,toon_coefficients='quadrature',opd=0.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "FigToon = jpi.rt_heatmap(Toon89, \n",
    "                        cmap_kwargs={'palette':jpi.pals.viridis(11), \n",
    "                                    'low':0,'high':0.8}, \n",
    "                        figure_kwargs={'title':'Toon89'})\n",
    "FigDlugach = jpi.rt_heatmap(Dlugach, \n",
    "                        cmap_kwargs={'palette':jpi.pals.viridis(11), \n",
    "                                    'low':0,'high':0.8} , \n",
    "                        figure_kwargs={'title':'Dlugach'})\n",
    "FigDiff = jpi.rt_heatmap((Toon89-Dlugach)/Dlugach, \n",
    "                        cmap_kwargs={'palette':jpi.pals.RdGy11, \n",
    "                                    'low':-0.4,'high':0.4} , \n",
    "                        figure_kwargs={'title':'% Diff'})\n",
    "jpi.show(jpi.row([FigToon, FigDlugach,FigDiff ]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A second way to visualize the albedo from using Toon, and Dlugach"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f = jpi.figure(x_range=[0.7,1], y_axis_label='Albedo',x_axis_label='Single Scattering Albedo',\n",
    "              height=300)\n",
    "for i,c in enumerate(Toon89.index): \n",
    "    f.line(Toon89.columns.astype(float), Toon89.loc[c,:], \n",
    "           color=jpi.pals.Spectral7[i],line_width=3)\n",
    "    f.line(Dlugach.columns.astype(float), Dlugach.loc[c,:], line_dash='dashed',\n",
    "           color=jpi.pals.Spectral7[i],line_width=3)\n",
    "jpi.show(f)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.11.5"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": true,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}