{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Adding Cloud Opacity, Scattering, Assymetry\n",
    "\n",
    "\n",
    "Clouds play a critical role in shaping your reflected light spectrum. \n",
    "\n",
    "Until it's release, we will guide you through the file format of the cloud input."
   ]
  },
  {
   "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",
    "\n",
    "jpi.output_notebook()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Start with the same inputs as before"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "opa = jdi.opannection(wave_range=[0.3,1])\n",
    "case1 = jdi.inputs()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#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)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## User Defined Cloud Input\n",
    "\n",
    "Just as with the Jupiter PT, we will use a Jupiter cloud profile as a demo. This should also be included in your base_cases reference file directory."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print( jdi.jupiter_cld())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#set model clouds\n",
    "case1.clouds( filename= jdi.jupiter_cld(), delim_whitespace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#run \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)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "jpi.show(jpi.spectrum(wno, alb, plot_width=500))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Huge** difference from our cloud free cases. High altitude water clouds have increased the albedo dramatically. Let's make sure we understand what is going on with the input file"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### File format for user defined input\n",
    "\n",
    "We define our cloud properties by : \n",
    "\n",
    "1) **g0** : Asymmetry factor of the cloud (value between 0 - 1)\n",
    "\n",
    "2) **w0** : Single scattering albedo (value between 0 - 1)\n",
    "\n",
    "3) **opd** : Total extinction (tau) of each layer (**not** cumulative) \n",
    "\n",
    "\n",
    "#### Cloud file specifications from a model (such as `virga`)\n",
    "\n",
    "1) Pressure and temperature are specified at the **boundaries** of each level. Cloud parameters are specified for each slab. Therefore, your cloud file should have **(number of PT points -1) layers**. \n",
    "\n",
    "2) Cloud files are pressure AND wavelength dependent. The file format could look something like this: \n",
    "\n",
    "    nlayer nwave pressure wavenumber w0 g0 opd\n",
    "    \n",
    "    1      1     #  #  #  #  # \n",
    "    \n",
    "    1      2     #  #  #  #  # \n",
    "    \n",
    "    1      3     #  #  #  #  # \n",
    "    \n",
    "    .      .     .  .  .  .  .\n",
    "    \n",
    "    1      nwno  .  .  .\n",
    "    \n",
    "    2      1     .  .  .\n",
    "    \n",
    "    2      2     .  .  .\n",
    "    \n",
    "    2      3     .  .  .\n",
    "    \n",
    "    .      .     .  .  .\n",
    "    \n",
    "    .      .     .  .  .\n",
    "    \n",
    "    .      .     .  .  .\n",
    "    \n",
    "    npt-1  nwno  .  .  .  .  .\n",
    "\n",
    "3)  It's assumed that these pressure grids correspond to your atmospheric grid. You will see errors through if #cld points != #atmo points - 1 \n",
    "\n",
    "4) You only have to speciy wavenumber grid if you aren't using `eddysed`\n",
    "\n",
    "Let's look at the premade example to make sure we understand. Since this is from `eddysed`, and since i know my pressure grid is the same as my previous input, I don't have to have columns for pressure and temperature"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Let's look at the premade example\n",
    "cld = pd.read_csv(jdi.jupiter_cld(), delim_whitespace=True)\n",
    "cld.head()\n",
    "#note, it is not required that lvl wv and sigma be in there"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plotting Cloud Input\n",
    "\n",
    "It's always good to make sure you intuitievly understand what is going on in your input before you create your spectrum. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nwno = 196 #this is the default number for A&M cloud code (see below if your wave grid is different)\n",
    "nlayer = 60 #one less than the number of PT points in your input \n",
    "jpi.show(jpi.plot_cld_input( nwno, nlayer,filename=jdi.jupiter_cld(),delim_whitespace=True))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The dark regions of the cloud optical depth plot specify where the bulk of the cloud opacity is. The light regions of the single scattering plot, indicate regions of wavenumber where your atmosphere is highly scattering."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cloud Input for Using Single Parameter Box Model\n",
    "\n",
    "Sometimes we are looking for something simpler than a full cloud model. In these cases we can speciy our cloud layer with `g0`, `w0`, `opd` and then `p` (the pressure level in log10 bars) and `dp` (the cloud thickness also in log10 bars)\n",
    "\n",
    "$\\mathrm{cloud\\_base(bars)} = 10^{p}$\n",
    "\n",
    "$\\mathrm{cloud\\_top(bars)} = 10^{p-dp}$\n",
    "\n",
    "**`dp` should never be negative!**\n"
   ]
  },
  {
   "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"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plotting simple cloud profile"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nwno = 196 #this is the default number for A&M cloud code (see below if your wave grid is different)\n",
    "nlayer = 60 #one less than the number of PT points in your input \n",
    "jpi.show(jpi.plot_cld_input(nwno, nlayer,df=case1.inputs['clouds']['profile']))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "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",
    "jpi.show(jpi.spectrum(wno, alb, plot_width=500))"
   ]
  }
 ],
 "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
}