{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Running the Code\n",
    "\n",
    "In this tutorial you will learn: \n",
    "\n",
    "1. How to run the code! \n",
    "\n",
    "You should already be familiar with: \n",
    "\n",
    "1. How to pick condensates to run in a cloud model \n",
    "2. What is $f_{sed}$ and what number is right for my model? \n",
    "3. What are the chemical limitations in the code \n",
    "4. How to compute initial Mie scattering grid\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from bokeh.io import output_notebook \n",
    "from bokeh.plotting import show, figure\n",
    "from bokeh.palettes import Colorblind\n",
    "output_notebook()\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import astropy.units as u\n",
    "\n",
    "#cloud code\n",
    "import virga.justdoit as jdi"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simple Isothermal/Constant $K_{z}$ Example"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's take the isothermal example from the first tutorial to begin gaining intuition for using the code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pressure = np.logspace(-5,3,30) #simple isotherml PT profile (kelvin)\n",
    "temperature = np.zeros(30)+1600  \n",
    "kz = np.zeros(30)+ 1e10 #constant kz profile \n",
    "\n",
    "metallicity = 1 #atmospheric metallicity relative to Solar\n",
    "mean_molecular_weight = 2.2 # atmospheric mean molecular weight\n",
    "\n",
    "#get pyeddy recommendation for which gases to run\n",
    "recommended_gases = jdi.recommend_gas(pressure, temperature,\n",
    "                                     metallicity,mean_molecular_weight)\n",
    "\n",
    "print(recommended_gases)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Define the `Atmosphere Class` "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#set cloud species to condense\n",
    "#let's take the recommended gases at face value for now\n",
    "sum_planet = jdi.Atmosphere(['MgSiO3'],fsed=0.1,mh=metallicity,\n",
    "                 mmw = mean_molecular_weight)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Set gravity and the P/T/Kz profile"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#set the planet gravity\n",
    "sum_planet.gravity(gravity=357.00, gravity_unit=u.Unit('cm/(s**2)'))\n",
    "\n",
    "#PT \n",
    "sum_planet.ptk(df = pd.DataFrame({'pressure':pressure, 'temperature':temperature,\n",
    "                           'kz':kz}))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#you can also read in a file using pandas.read_csv() \n",
    "pd.DataFrame({'pressure':pressure, 'temperature':temperature,\n",
    "                           'kz':kz}).to_csv('test.csv')\n",
    "sum_planet.ptk(filename='test.csv',usecols = [1,2,3])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### The minimum values for `kz` \n",
    "\n",
    "There is a parameter called `kz_min` in the ptk function. The default value is `kz_min=1e5`. This number isn't necessarily a hard and fast rule, but it does avoid numerical instabilities. Another thing to note is that __this excludes negative kz values__. The code will alert you if you are tripping the `kz_min` reset. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#set the planet gravity\n",
    "sum_planet.gravity(gravity=357.00, gravity_unit=u.Unit('cm/(s**2)'))\n",
    "\n",
    "#This is what happens when you trip up the minimum KZ value\n",
    "sum_planet.ptk(df = pd.DataFrame({'pressure':pressure, 'temperature':temperature,\n",
    "                           'kz':kz*-1}))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Run the code "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#directory where mieff files are \n",
    "mieff_directory = '/data/virga/'\n",
    "sum_planet.ptk(filename='test.csv',usecols = [1,2,3])\n",
    "#start with the simplest output \n",
    "#get total opacity, single scattering, asymmetry, and individual optical depths \n",
    "df_out = sum_planet.compute(directory=mieff_directory)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#get full dictionary output \n",
    "all_out = sum_planet.compute(directory=mieff_directory, as_dict=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exploring `dict` Output\n",
    "\n",
    "There are many things to explore in the `dict` output. In the next tutorial, we will reproduce some of the most common plots used in papers to analyze cloud runs. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#see what outputs exist in the dictionary \n",
    "all_out.keys()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hopefully the names of the `dict` elements are self explanatory. If not, here is a brief description of each. `nl`=number of layers, `nw`=number of wavelengths, `ng`=number of gases.\n",
    "\n",
    "- `condensate_mmr`(nl x ng): Mean mass mixing ratio (concentration) of the condensate. See `qc` in A&M 01 Eqn. 4 and Eqn. 8.  \n",
    "- `cond_plus_gas_mmr`(nl x ng): Mean mass mixing ratio of the consensate plus the vapor. See `qt` in A&M 01 Eqn. 7. \n",
    "- `mean_particle_r`(nl x ng): Geometric mean particle radius. See `r_g` in A&M 01 Eqn. 13. \n",
    "- `droplet_eff_r`(nl x ng): Effetive (area-weighted) droplet radius. See `r_eff` in A&M 01 Eqn. 17\n",
    "- `column_density`(nl x ng): Total number concentration of particles PER LAYER. See `N` in A&M 01 Eqn. 14\n",
    "- `opd_per_layer`(nl x nw): Total extinction PER layer. This includes all gases. \n",
    "- `single_scattering`(nl x nw): Total single scattering albedo \n",
    "- `asymmetry`(nl x nw): Total asymmetry \n",
    "- `opd_by_gas`(nl x ng): Optical depth for conservative geometric scatteres separated by gas. E.g. [Fig 7 in Morley+2012](https://arxiv.org/pdf/1206.4313.pdf)  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Additionally, scalar input contains some of the original input "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_out['scalar_inputs']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}