{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Getting Set up to Run the Code\n",
    "\n",
    "In this tutorial you will learn: \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",
    "\n",
    "If you are already familiar with these things, I recommend moving to the next tutorial"
   ]
  },
  {
   "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",
    "#here is virga\n",
    "import virga.justdoit as jdi"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will be pretty thorough in this first tutorial. In the very first step, we will set up our calculation run. This includes: \n",
    "\n",
    "1. Defining what gas condensates we want to include in the calculation. \n",
    "\n",
    "2. Defining the sedimentation efficiency $f_{sed}$. \n",
    "\n",
    "3. Defining metallicity and atmospheric mean molecular weight\n",
    "\n",
    "3. Defining grid of particle radii for the Mie code\n",
    "\n",
    "Let's briefly go through these one-by-one."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## How to Pick Gas Condensates \n",
    "\n",
    "We highly recommend looking at your climate profile first to give an informed decision to the cloud model. First check out what current condensates are available"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#you can print the available condensate species in the code\n",
    "jdi.available()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then you can let `virga` give you a recommendation with some visualization aid"
   ]
  },
  {
   "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)+1300  \n",
    "metallicity = 1 #atmospheric metallicity relative to Solar\n",
    "mean_molecular_weight = 2.2 # atmospheric mean molecular weight\n",
    "\n",
    "#get virga recommendation for which gases to run\n",
    "recommended = jdi.recommend_gas(pressure, temperature, metallicity,mean_molecular_weight, \n",
    "                #Turn on plotting \n",
    "                 plot=True)\n",
    "#print the results\n",
    "print(recommended)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### How to Compute Temperature Condensation Curves"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`jdi.recommend` is comparing your pressure-temperature curve against __condensation temperature curves__ for each species. If you are interested in computing your own you can do so!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cond_t_fig = figure(y_axis_type='log',y_range=[1e2,1e-3],height=400, width=400,\n",
    "                  y_axis_label='Pressure(bars)',x_axis_label='Temperature (K)',\n",
    "                   title='Condensation Temperatures')\n",
    "metallicity = 1 \n",
    "mean_molecular_weight = 2.2\n",
    "for gas_name, c in zip(['Cr','MgSiO3','MnS'],Colorblind[8]): #case sensitive names\n",
    "    #grab p,t from eddysed\n",
    "    p,t = jdi.condensation_t(gas_name, metallicity, mean_molecular_weight)\n",
    "    #plot it up\n",
    "    cond_t_fig.line(t,p, legend_label=gas_name,color=c,line_width=3)\n",
    "\n",
    "show(cond_t_fig)    \n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## How to pick $f_{sed}$\n",
    "\n",
    "This parameter will set the vertical extent of our cloud. The general equation is: \n",
    "\n",
    "`condensate concentration (z)` $\\propto$ `concentration at cloud base` $exp(-f_{sed}z/H)$\n",
    "\n",
    "Therefore, high values of $f_{sed}$ deplete the cloud quickly, low values of $f_{sed}$ create large vertically thick clouds. __For context, it is thought that Hot Jupiters have fractional $f_{sed}$ values around 0.1, while Jupiter clouds can be parameterized with much larger values of 3-6.__\n",
    "\n",
    "Here is a little schematic to show the general behavior of $f_{sed}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#km, arbitrary planet scale height \n",
    "H = 27 \n",
    "\n",
    "#define the cloud base \n",
    "cloud_base = 0.8 #bars  \n",
    "\n",
    "#and some arbitrary height\n",
    "z = np.linspace(0,300,100)#km \n",
    "pressure_grid=100*np.exp(-z/H)#bars \n",
    "\n",
    "#setup figure\n",
    "toy_model = figure(x_axis_type='log',x_range=[1e-3,1.1],y_range=[1,1e-2],height=400, width=400,\n",
    "                  y_axis_label='Pressure(bars)',x_axis_label='Arbitrary Fake Cloud Concentration Units',\n",
    "                  title='F_sed Toy Model')\n",
    "\n",
    "#loop through some fseds and plot them up \n",
    "for fsed,c in zip([10,6,3,1,0.1],['red','orange','green','blue','purple']):\n",
    "    cloud=np.zeros(100)+1e-13 #start with empty array \n",
    "    \n",
    "    #exponential equation I mentioned above \n",
    "    fake_cloud = np.exp(-fsed*z[pressure_grid<cloud_base]/H) #fill with exponential fake cloud \n",
    "    cloud[pressure_grid<cloud_base] = fake_cloud/np.max(fake_cloud) #fill cloud with fake cloud and normalize\n",
    "    toy_model.line(cloud,pressure_grid,color=c,legend_label='F_sed='+str(fsed))\n",
    "    toy_model.legend.location='bottom_left'\n",
    "show(toy_model)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Chemical Limitations of the Code (e.g. M/H)\n",
    "\n",
    "Although we are trying to make `virga` as general as possible. There are limitations to the metallicity and mean molecular weights that users enter. \n",
    "\n",
    "#### Retrieving pre-set gas properties \n",
    "\n",
    "__There are two values that are hard-coded for H2 atmospheres__\n",
    "\n",
    "1. The diameter of the atmospheric molecule (We assume d = 2.827e-8 for H2) \n",
    "2. The depth of the Lennard-Jones potential well for the atmosphere (we use 59.7 K, the H2 value)\n",
    "\n",
    "You may ask... *Seems REAL easy to make these an input.. What gives!?* \n",
    "\n",
    "The reason we fix this is because the chemical concentrations (taken from Morley et al 2012) cannot be trivially scaled and applied to non-H2 very high metallicity atmospheres. \n",
    "\n",
    "To see what gas properties are currently in the code, you can explore `gas_properties`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "atmo_mean_molecular_weight = 2.2 \n",
    "metallicity = 1 #non-log, relative to solar \n",
    "\n",
    "#retrieve : \n",
    "#condensate molecular weight, gas mass mixing ratio, density of the condensate\n",
    "gas_mw, gas_mmr, rho_p = jdi.gas_properties.H2O(atmo_mean_molecular_weight, metallicity )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Some of our functions have a metallicity dependence (such as ZnS)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#of course there are other gases which do contain metallicity dependencies\n",
    "#e.g. ZnS \n",
    "gas_mw, gas_mmr, rho_p = jdi.gas_properties.ZnS(atmo_mean_molecular_weight, 100*metallicity )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Creating the Mie scattering Database \n",
    "\n",
    "The wavelength dependence of the code relies on computing the scattering and absorbing coefficients from Mie theory. \n",
    "\n",
    "Our python code uses [PyMieScatt](https://pymiescatt.readthedocs.io/en/latest/), which has [extensive documentation](https://pymiescatt.readthedocs.io/en/latest/forward.html) on the how the calculations are done assuming that all particles are spherical and homogeneous. \n",
    "\n",
    "What you need to know is that before running `virga` you must first create a database of these parameters given a particles complex refractive index ($m = n+ik$), which are basically just a condensate species optical properties. These optical properties are highly sought after properties, which are measured in the lab. \n",
    "\n",
    "### Complex refractive indices\n",
    "\n",
    "$n$ and $k$ are a function of wavelength. Mie scattering and absorbing coefficients are a function of wavelength and particle size. \n",
    "\n",
    "You can learn a lot about a particle itself by exploring the complex refractive indices. Let's look at both.\n",
    "\n",
    "`virga` comes with a folder with files that end in `.refrind`. We've fixed the grid of these refractive indices. However, note that the wavelength grid needs to be fine enough to resolve any spectral features you are interested in."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#directory where you have put `refrind`\n",
    "refrind_dir = '/data/virga'\n",
    "wave, n, k  = jdi.get_refrind('Na2S', refrind_dir)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lets reproduce Figure 2 in [Morley et al 2012](https://arxiv.org/pdf/1206.4313.pdf).\n",
    "\n",
    "The single __most useful thing to remember about refractive indices is:__ imaginary part (also called the extinction coefficient) relates to the absorbing power of the species. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from bokeh.layouts import row\n",
    "k_fig = figure(height = 200,width=300 ,y_axis_label = 'Imaginary Part',x_axis_type='log',y_axis_type='log',x_axis_label='Wave (um)') \n",
    "n_fig = figure(height = 200, width=300,y_axis_label = 'Real Part',x_axis_type='log',x_axis_label='Wave (um)') \n",
    "n_fig.line(wave,n)\n",
    "k_fig.line(wave,k)\n",
    "show(row([n_fig, k_fig]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Pre-set wavelength Grid\n",
    "\n",
    "Note that we have supplied $n$ and $k$ on a pre-set wavelength grid. Some people have inquired about __increasing the resolution of the grid__. Doing so requires more robust laboratory data. Ask your local theorist what might be available for your problem. \n",
    "\n",
    "### Running the Mie Database Code\n",
    "\n",
    "The function below will create a set of files called e.g. `ZnS.Mieff` in whatever directory you specify. So beware of file overwrites. \n",
    "\n",
    "__Just for example, we will run a single gas with a 10 radius grid. BUT IRL USE AT LEAST THE DEFAULT nradii=40 and for all the gases__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "out_dir = os.getcwd()\n",
    "\n",
    "#JUST FOR THE EXAMPLE, let's run one species on a small radii grid \n",
    "gas_name = 'H2O'\n",
    "qext, qscat, g_qscat, radii,wave = jdi.calc_mie_db([gas_name], \n",
    "                                refrind_dir, out_dir, rmin = 1e-5, nradii = 10) \n",
    "\n",
    "#Your IRL run, which will take some time to compute\n",
    "#qext, qscat, cos_qscat, radii = ed.calc_mie_db(ed.avaialable(), refrind_dir, out_dir) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Analyzing Mie Parameters\n",
    "\n",
    "For reflected light calculations, the most common/useful plots to look at are \n",
    "\n",
    "- __Asymmetry__ (`g_qscat/qscat`) : will show how forward/backscattering the gas is \n",
    "- __Qscat / Qext__ (`qscat/qext`) : will show you how bright the condensate is expected to be"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from bokeh.palettes import viridis\n",
    "from bokeh.models import LogColorMapper, ColorBar, LogTicker\n",
    "\n",
    "colors = viridis(len(radii))\n",
    "qfig = figure(width=300, height=300,\n",
    "              x_axis_type='log',y_axis_label ='Asymmetry',\n",
    "              x_axis_label='Wavelength(um)',title=gas_name)\n",
    "\n",
    "wfig = figure(width=400, height=300,\n",
    "              x_axis_type='log',y_axis_label ='Qscat/Qext',\n",
    "              x_axis_label='Wavelength(um)',title=gas_name)\n",
    "\n",
    "for i in range(len(radii)):\n",
    "    qfig.line(wave, g_qscat[:,i,0]/qscat[:,i,0], color=colors[i])\n",
    "    wfig.line(wave, qscat[:,i,0]/qext[:,i,0], color=colors[i])\n",
    "    \n",
    "#color bar\n",
    "color_mapper = LogColorMapper(palette=\"Viridis256\", low=min(radii*1e4), high=max(radii*1e4))\n",
    "color_bar = ColorBar(color_mapper=color_mapper, ticker=LogTicker(),\n",
    "                     label_standoff=12, border_line_color=None, location=(0,0),title=\"Radius (um)\")\n",
    "wfig.add_layout(color_bar, 'right')\n",
    "\n",
    "show(row(qfig, wfig))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can also quickly grab your files after you've computed them. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qext, qscat, g_qscat, nwave,radii,wave = jdi.get_mie(gas_name,out_dir)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qext.shape"
   ]
  },
  {
   "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
}