{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Approximations for Spherical Harmonics Radiative Transfer in Reflected Light\n",
    "\n",
    "In [Rooney et al 2023](add-link) we rigorously derive the spherical harmonics method for reflected light and benchmark the 4-term method (SH4) against [Toon et al. 1989](https://ui.adsabs.harvard.edu/abs/1989JGR....9416287T/abstract) and two independent methods. Here, we provide the code to reproduce the analysis that compares Toon89 with the higher fidelity 4-term spherical harmonics method for reflected light calculations. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "#picaso\n",
    "from picaso import justdoit as jdi \n",
    "from picaso import justplotit as jpi\n",
    "import picaso.test as ptest\n",
    "\n",
    "jpi.output_notebook()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Benchmark Jupiter-like case profile\n",
    "\n",
    "Here we use the same profile explored in the two-stream radiative transfer tutorial. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "opa = jdi.opannection(wave_range=[0.3,1])#, resample=100)\n",
    "case1 = jdi.inputs()\n",
    "#phase \n",
    "case1.phase_angle(0) #radians\n",
    "\n",
    "#gravity \n",
    "case1.gravity(gravity = 25, gravity_unit=jdi.u.Unit('m/(s**2)'))\n",
    "\n",
    "#star\n",
    "case1.star(opa, 6000,0.0122,4.437) #kelvin, log metal, log cgs\n",
    "\n",
    "#atmosphere\n",
    "case1.atmosphere(filename= jdi.jupiter_pt(), delim_whitespace=True)\n",
    "\n",
    "#set model clouds\n",
    "case1.clouds( filename= jdi.jupiter_cld(), delim_whitespace=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setup comparison with two-stream radiative transfer "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#we'll use these labels to keep track of the cases we have created\n",
    "labels = []\n",
    "albs = []\n",
    "#run\n",
    "\n",
    "multi_phase = 'N=2' #two legendre polynomials (toon default)\n",
    "single_phase='TTHG_ray' #two term HG with rayleigh (toon default)\n",
    "DE = True #delta eddington correction, (toon default)\n",
    "raman='none' #lets turn raman off, for clarity in this benchmark case\n",
    "\n",
    "#set all approximatinos\n",
    "case1.approx(toon_coefficients='quadrature', \n",
    "             multi_phase=multi_phase,\n",
    "             single_phase=single_phase,\n",
    "             delta_eddington=DE, raman=raman)\n",
    "df = case1.spectrum(opa)\n",
    "wno, alb, fpfs = df['wavenumber'] , df['albedo'] , df['fpfs_reflected'] \n",
    "wno, alb = jdi.mean_regrid(wno, alb, R=150)\n",
    "\n",
    "labels+=['Toon89']\n",
    "albs+=[alb]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setting up Spherical Harmonics Approximations to compare with Toon\n",
    "\n",
    "Spherical harmonics allows us to remain consistent with scattering functions throughout the methodology. In Toon when calculation the two stream solution for multiple layers, the phase functions are hard-coded set to be one term HG. However, when implementing the source function technique to derive the outgoing intensity we introduce a two-term HG for direct-scattering beam in attempt to capture the back-scattering radiation observed on Neptune. These details are further discussed in Cahoy et al 2010, Batalha et al. 2019, Feng et al. 2017. \n",
    "\n",
    "In spherical harmonics, we can simply consider a two-term HG phase function throughout the calculation. In order to better compare with Toon however (for historical purposes), we made sure that our SH routines could also match that Toon approach. Therefore, below you see we have **two** different single scattering forms. Note that this is only for comparison with Toon89. The default for SH is to use a TTHG phase function with Rayleigh throughout the radiative transfer calculation. \n",
    "\n",
    "Further details can be found in Rooney et al. 2023. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# NOT SH default; enforcing for toon comparison which cannot handle \n",
    "# pre-processed two term phase functions\n",
    "w_single_form='OTHG'; w_single_rayleigh='off'\n",
    "w_multi_form='OTHG'; w_multi_rayleigh='on'\n",
    "\n",
    "#defaults in SH technique\n",
    "psingle_form='TTHG'; psingle_rayleigh='on'\n",
    "#2 or 4 stream (SH routines handles both though 4 is the default)\n",
    "stream=4\n",
    "\n",
    "case1.approx(rt_method='SH', stream=stream, \n",
    "              w_single_form=w_single_form, w_single_rayleigh=w_single_rayleigh,\n",
    "              w_multi_form=w_multi_form, w_multi_rayleigh=w_multi_rayleigh,\n",
    "             psingle_form=psingle_form, psingle_rayleigh=psingle_rayleigh,\n",
    "             delta_eddington=DE, raman=raman)\n",
    "# case1.approx(rt_method='SH', stream=4, raman='none')\n",
    "df3 = case1.spectrum(opa)\n",
    "wno3, alb3, fpfs3 = df3['wavenumber'] , df3['albedo'] , df3['fpfs_reflected'] \n",
    "wno3, alb3 = jdi.mean_regrid(wno3, alb3, R=150)\n",
    "labels+=['SH4 (OTHG multi)']\n",
    "albs+=[alb3]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we relax the constraint for OTHG to see how the multile scattering approximation in Toon impacts the spectrum.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "w_multi_form='TTHG'; w_multi_rayleigh='on'\n",
    "\n",
    "case1.approx(rt_method='SH', stream=4, \n",
    "              w_single_form=w_single_form, w_single_rayleigh=w_single_rayleigh,\n",
    "              w_multi_form=w_multi_form, w_multi_rayleigh=w_multi_rayleigh,\n",
    "             psingle_form=psingle_form, psingle_rayleigh=psingle_rayleigh,\n",
    "             delta_eddington=DE, raman='none')\n",
    "# case1.approx(rt_method='SH', stream=4, raman='none')\n",
    "df5 = case1.spectrum(opa)\n",
    "wno5, alb5, fpfs5 = df5['wavenumber'] , df5['albedo'] , df5['fpfs_reflected'] \n",
    "wno5, alb5 = jdi.mean_regrid(wno5, alb5, R=150)\n",
    "labels+=['SH4 (TTHG multi)']\n",
    "albs+=[alb5]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Reproducing Figure 5a from Rooney et al. 2023"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "jpi.show(jpi.spectrum([wno]*3, albs, labels, width=700))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Reproducing Figure 5b from Rooney et al. 2023\n",
    "\n",
    "Here we reproduce the same procedure from above, but with a different cloud setup."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#set model clouds, note these are lists since you can specify multiple cloud layers\n",
    "case1.clouds( g0=[0.9], w0=[0.8], opd=[0.5], p = [0.0], dp=[1.0])  # Slab cloud from 1.0 bar up to 0.1 bar\n",
    "labels = []\n",
    "albs = []"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "multi_phase = 'N=2'\n",
    "single_phase='TTHG_ray'\n",
    "DE = True\n",
    "\n",
    "case1.approx(multi_phase=multi_phase,single_phase=single_phase,\n",
    "             delta_eddington=DE, raman='none')\n",
    "df = case1.spectrum(opa)\n",
    "wno, alb, fpfs = df['wavenumber'] , df['albedo'] , df['fpfs_reflected'] \n",
    "wno, alb = jdi.mean_regrid(wno, alb, R=150)\n",
    "labels+=['Toon89']\n",
    "albs+=[alb]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#run SH4\n",
    "w_single_form='OTHG'; w_single_rayleigh='off' \n",
    "w_multi_form='OTHG'; w_multi_rayleigh='on'\n",
    "psingle_form='TTHG'; psingle_rayleigh='on'\n",
    "\n",
    "case1.approx(rt_method='SH', stream=4, \n",
    "              w_single_form=w_single_form, w_single_rayleigh=w_single_rayleigh,\n",
    "              w_multi_form=w_multi_form, w_multi_rayleigh=w_multi_rayleigh,\n",
    "             psingle_form=psingle_form, psingle_rayleigh=psingle_rayleigh,\n",
    "             delta_eddington=DE, raman='none')\n",
    "# case1.approx(rt_method='SH', stream=4, delta_eddington=DE, raman='none')\n",
    "df2 = case1.spectrum(opa)\n",
    "wno2, alb2, fpfs2 = df2['wavenumber'] , df2['albedo'] , df2['fpfs_reflected'] \n",
    "wno2, alb2 = jdi.mean_regrid(wno2, alb2, R=150)\n",
    "labels+=['SH4 (OTHG single)']\n",
    "albs+=[alb2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#run SH2\n",
    "w_multi_form='TTHG'; w_multi_rayleigh='on'\n",
    "\n",
    "case1.approx(rt_method='SH', stream=4, \n",
    "              w_single_form=w_single_form, w_single_rayleigh=w_single_rayleigh,\n",
    "              w_multi_form=w_multi_form, w_multi_rayleigh=w_multi_rayleigh,\n",
    "             psingle_form=psingle_form, psingle_rayleigh=psingle_rayleigh,\n",
    "             delta_eddington=DE, raman='none')\n",
    "df3 = case1.spectrum(opa)\n",
    "wno3, alb3, fpfs3 = df3['wavenumber'] , df3['albedo'] , df3['fpfs_reflected'] \n",
    "wno3, alb3 = jdi.mean_regrid(wno3, alb3, R=150)\n",
    "labels+=['SH4 (TTHG multi)']\n",
    "albs+=[alb3]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "jpi.show(jpi.spectrum([wno]*3, albs,labels, width=700))"
   ]
  }
 ],
 "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.8.12"
  },
  "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": 2
}