{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# How clouds affect albedo spectroscopy\n",
    "\n",
    "`PICASO` is used to compute reflected light and thermal emission of exoplanets. It is not a `virga` dependency though. This tutorial is mostly used as a module to help people gain an intuition for how various cloud outputs will affect reflected light observations.\n",
    "\n",
    "What you will learn: \n",
    "\n",
    "1. How different cloud models affect reflected light spectra\n",
    "\n",
    "What you should already know: \n",
    "\n",
    "1. How to run `PICASO` via [these tutorials](https://natashabatalha.github.io/picaso/notebooks/1_GetStarted.html). \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "#main programs\n",
    "from picaso import justdoit as pj\n",
    "from virga import justdoit as vj\n",
    "#plot tools\n",
    "from picaso import justplotit as picplt\n",
    "from virga import justplotit as cldplt\n",
    "\n",
    "\n",
    "import astropy.units as u\n",
    "import pandas as pd\n",
    "from bokeh.plotting import show, figure\n",
    "from bokeh.io import output_notebook \n",
    "output_notebook()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's set up a basic spectrum plot "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "opacity = pj.opannection(wave_range=[0.3,1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sum_planet = pj.inputs()\n",
    "sum_planet.phase_angle(0) #radians\n",
    "sum_planet.gravity(gravity=25, gravity_unit=u.Unit('m/(s**2)')) #any astropy units available\n",
    "sum_planet.star(opacity, 5000,0,4.0) #opacity db, pysynphot database, temp, metallicity, logg\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's do simple cloud run with the `jupiter_pt` profile and constant $K_z$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_atmo = pd.read_csv(pj.jupiter_pt(), delim_whitespace=True)\n",
    "#you will have to add kz to the picaso profile\n",
    "df_atmo['kz'] = [1e9]*df_atmo.shape[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#business as usual\n",
    "sum_planet.atmosphere(df=df_atmo)\n",
    "\n",
    "#let's get the cloud free spectrum for reference\n",
    "cloud_free = sum_planet.spectrum(opacity)\n",
    "\n",
    "x_cld_free, y_cld_free = pj.mean_regrid(cloud_free['wavenumber'], cloud_free['albedo'], R=150)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Submit `virga` run within `picaso`\n",
    "\n",
    "You will be able to run `virga` runs directly through `picaso` so you don't have to go through all the steps twice. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "metallicity = 1 #atmospheric metallicity relative to Solar\n",
    "mean_molecular_weight = 2.2 # atmospheric mean molecular weight\n",
    "directory ='/data/virga/'\n",
    "\n",
    "#we can get the same full output from the virga run\n",
    "cld_out = sum_planet.virga(['H2O'],directory, fsed=1,mh=metallicity,\n",
    "                 mmw = mean_molecular_weight)\n",
    "out = sum_planet.spectrum(opacity, full_output=True)\n",
    "\n",
    "x_cldy, y_cldy = pj.mean_regrid(out['wavenumber'], out['albedo'], R=150)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analyzing Cloudy versus Cloud-free\n",
    "\n",
    "Let's take a moment to understand the difference between the cloudy and the cloud free case "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "show(picplt.spectrum([x_cld_free, x_cldy],\n",
    "                     [y_cld_free, y_cldy],plot_width=500, plot_height=300,\n",
    "                  legend=['Cloud Free','Cloudy']))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Brighter or Darker?\n",
    "\n",
    "Our spectra clearly got much brighter towards 1 micron. Let's explore why this is by looking at the optical properties and the photon attenuation plot. \n",
    "\n",
    "#### Photon Attenuation Depth\n",
    "\n",
    "This will tell us at what pressure we attain $\\tau \\sim 1$ optical depth for the molecular opacity, cloud opacity, and rayleigh opacity. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "show(picplt.photon_attenuation(out['full_output'],\n",
    "                            plot_width=500, plot_height=300))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If there were no cloud present, in this case Rayleigh scattering would dominate at the bluer wavelengths and gas opacity at the redder. But we see that the cloud reaches a significant optical depth well above the other two, so it will certainly influence the computed spectrum. Therefore, we know that our spectrum is going to be heavily influenced by this particular cloud model. **But now we need to know what the optical properties of the cloud are at those pressures** \n",
    "\n",
    "Let's ask ourselves the following questions: \n",
    "\n",
    "1. What are the particle radii at 0.1 bar where the cloud opacity is sitting? \n",
    "2. What are the Mie properties (asymmetry and single scattering) or the main condensate at that particle radii?\n",
    "\n",
    "#### 1) Look at particle radii\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, dndr = cldplt.radii(cld_out,at_pressure=0.1)\n",
    "show(fig)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Eyeballing the plot, it looks like our mean particle radius is about 30 microns. Let's pull our Mie parameters and see what the single scattering and asymmetry look like at that point. \n",
    "\n",
    "#### 2) Look at Mie parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#grab your mie parameters\n",
    "qext, qscat, g_qscat, nwave,radii,wave = vj.get_mie('H2O',directory)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from bokeh.layouts import row,column\n",
    "ind = cldplt.find_nearest_1d(radii,30e-4) #remember the radii are in cm\n",
    "\n",
    "qfig = figure(width=300, height=300,\n",
    "              x_axis_type='log',y_axis_label ='Asymmetry',\n",
    "              x_axis_label='Wavelength(um)')\n",
    "\n",
    "wfig = figure(width=300, height=300,\n",
    "              x_axis_type='log',y_axis_label ='Qscat/Qext',\n",
    "              x_axis_label='Wavelength(um)')\n",
    "\n",
    "qfig.line(wave[:,0], g_qscat[:,ind]/qscat[:,ind])\n",
    "wfig.line(wave[:,0], qscat[:,ind]/qext[:,ind])\n",
    "\n",
    "show(row(qfig, wfig))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This plot covers a large wavelength range. Zoom in on 0.3,1 micron, where we are computing our spectrum. **What asymmetry and single scattering values do you observe?**. \n",
    "\n",
    "Single scattering is essentially 1!! Water clouds are highly scattering! The asymmetry is around 0.8, which is quite forward scattering. However, the **high single scattering** creates a very bright cloud deck."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## How varying cloud species affects reflectance spectrum\n",
    "\n",
    "Let's ramp up the temperature a bit and move away from H2O clouds to Na2S and ZnS. **WARNING**: We are about to artificially increase the temperature. Don't try this at home kids. IRL I would have to change all my chemistry as well. But here I am just making a point about how cloud composition/temp affects the reflectance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hot_atmo = df_atmo\n",
    "hot_atmo['temperature'] = hot_atmo['temperature'] + 600\n",
    "\n",
    "#remember we can use recommend_gas function to look at what the condensation curves look like\n",
    "recommended = vj.recommend_gas(hot_atmo['pressure'], hot_atmo['temperature'], metallicity,mean_molecular_weight,\n",
    "                              plot=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Run `picaso` and `virga` with different cloud species"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#business as usual\n",
    "sum_planet.atmosphere(df=hot_atmo)\n",
    "\n",
    "#make sure clouds are turned off\n",
    "sum_planet.clouds_reset()\n",
    "\n",
    "#let's get the cloud free spectrum for reference\n",
    "cloud_free = sum_planet.spectrum(opacity)\n",
    "x_cld_free, y_cld_free = pj.mean_regrid(cloud_free['wavenumber'], cloud_free['albedo'], R=150)\n",
    "\n",
    "#now the cloudy runs\n",
    "cld_out = sum_planet.virga(['Na2S','ZnS'],directory, fsed=1,mh=metallicity,\n",
    "                 mmw = mean_molecular_weight)\n",
    "\n",
    "out = sum_planet.spectrum(opacity, full_output=True)\n",
    "x_cld, y_cld = pj.mean_regrid(out['wavenumber'], out['albedo'], R=150)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compare Na2S/ZnS against H2O clouds\n",
    "\n",
    "Look at the huge difference here. Our Na2S/ZnS did not have the effect the water clouds did. Let's look at photon attenuation, particle radii and Mie parameters to see if we can deduce why. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "w = [x_cld_free, x_cld]\n",
    "a = [y_cld_free, y_cld]\n",
    "show(picplt.spectrum(w,a,plot_width=500, plot_height=300,\n",
    "                  legend=['Cloud Free','Cloudy']))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "#### Is the cloud opacity at a different pressure than the water cloud? "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "show(picplt.photon_attenuation(out['full_output'],\n",
    "                            plot_width=500, plot_height=300))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Yes! The cloud opacity is quite a bit lower than our H2O run. \n",
    "\n",
    "#### Are the particle radii different?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, dndr = cldplt.radii(cld_out,at_pressure=0.5)\n",
    "show(fig)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Not drastically different. We have particles that are approximately 10 microns instead of 30 microns. This is probably not a huge driver.\n",
    "\n",
    "#### Are the optical properties different? Run this script for Na2S and ZnS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#grab your mie parameters\n",
    "gas_name = 'Na2S' #ZnS\n",
    "qext, qscat, g_qscat, nwave,radii,wave = vj.get_mie(gas_name,directory)\n",
    "ind = cldplt.find_nearest_1d(radii,10e-4) #remember the radii are in cm\n",
    "\n",
    "qfig = figure(width=300, height=300,\n",
    "              x_axis_type='log',y_axis_label ='Asymmetry',\n",
    "              x_axis_label='Wavelength(um)')\n",
    "\n",
    "wfig = figure(width=300, height=300,\n",
    "              x_axis_type='log',y_axis_label ='Qscat/Qext',\n",
    "              x_axis_label='Wavelength(um)')\n",
    "\n",
    "qfig.line(wave[:,0], g_qscat[:,ind]/qscat[:,ind])\n",
    "wfig.line(wave[:,0], qscat[:,ind]/qext[:,ind])\n",
    "\n",
    "show(row(qfig, wfig))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Yes!! What did you notice about the optical properties between 0.3-1 micron? While these condensates are significantly less forward scattering, the single scattering albedos are much lower than one. \n",
    "\n",
    "**In conclusion! Your reflectance spectra are going to drastically change depending on what condensate species you are running!!!** For this case the higher pressure and lower single scattering albedos were the culprit in decreasing the overall brightness of the spectrum."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## How sedimentation efficiency affects albedo spectrum?\n",
    "\n",
    "Let's go back to your normal cool case with H2O clouds. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_atmo = pd.read_csv(pj.jupiter_pt(), delim_whitespace=True)\n",
    "df_atmo['kz'] = [1e9]*df_atmo.shape[0]\n",
    "\n",
    "sum_planet.atmosphere(df = df_atmo)\n",
    "\n",
    "all_fseds = [0.1, 1, 6, 10]\n",
    "w = []\n",
    "a = []\n",
    "all_outs,cld_outs=[],[]\n",
    "for fs in all_fseds:\n",
    "    cld = sum_planet.virga(['H2O'],directory, fsed=fs,mh=metallicity,\n",
    "                 mmw = mean_molecular_weight)\n",
    "    cld_outs += [cld]\n",
    "    out = sum_planet.spectrum(opacity,full_output=True)\n",
    "    x,y = pj.mean_regrid(out['wavenumber'], out['albedo'], R=150)\n",
    "    w += [x]\n",
    "    a += [y]\n",
    "    all_outs += [out['full_output']]\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "show(picplt.spectrum(w,a,plot_width=500, plot_height=300,\n",
    "                     legend=['fs= '+str(i) for i in all_fseds]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### For H2O, why does lower $f_{sed}$ lead to higher albedo ?\n",
    "\n",
    "We can use the same exercises that we have used for the previous exercise in order to assess what is going on. \n",
    "\n",
    "For the sake of less code, let's just compare the $f_{sed}$ 0.1 and 6 case. \n",
    "\n",
    "### How $f_{sed}$ affects cloud optical depth"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "show(column(picplt.photon_attenuation(all_outs[0],title='fs=0.1',\n",
    "                            plot_width=500, plot_height=300),\n",
    "   picplt.photon_attenuation(all_outs[2],title='fs=6',\n",
    "                            plot_width=500, plot_height=300)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Zoom out so you can see where the cloud opacity is at low fsed! The cloud opacity at low fsed is about an order of magnitude lower pressure than the high fsed case. In the high fsed case, at longer wavelengths, the molecular opacity starts competing with the cloud opacity. \n",
    "\n",
    "### How $f_{sed}$ affects particle size "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#low fsed case\n",
    "fig, dndr = cldplt.radii(cld_outs[0],at_pressure=1e-3)\n",
    "show(fig)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#high fsed case\n",
    "fig, dndr = cldplt.radii(cld_outs[2],at_pressure=1)\n",
    "show(fig)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Mean particle sizes where the cloud is most optically thick have changed significantly (almost two orders of magnitude). Therefore, even though we have the same cloud species, the Mie parameters will be different!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#grab your mie parameters\n",
    "gas_name = 'H2O'\n",
    "qext, qscat, g_qscat, nwave,radii,wave = vj.get_mie(gas_name,directory)\n",
    "\n",
    "ind_low = cldplt.find_nearest_1d(radii,1e-4) #look at the plots to get these numberss\n",
    "ind_high = cldplt.find_nearest_1d(radii,50e-4)#remember radii is in cm\n",
    "\n",
    "qfig = figure(width=300, height=300,\n",
    "              x_axis_type='log',y_axis_label ='Asymmetry',\n",
    "              x_axis_label='Wavelength(um)')\n",
    "\n",
    "wfig = figure(width=300, height=300,\n",
    "              x_axis_type='log',y_axis_label ='Qscat/Qext',\n",
    "              x_axis_label='Wavelength(um)')\n",
    "\n",
    "for ind,l,c in zip([ind_low,ind_high],['fsed=0.1','fsed=6'],['red','blue']):\n",
    "    qfig.line(wave[:,0], g_qscat[:,ind]/qscat[:,ind],legend_label=l,color=c)\n",
    "    wfig.line(wave[:,0], qscat[:,ind]/qext[:,ind],legend_label=l,color=c)\n",
    "\n",
    "show(row(qfig, wfig))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## How $K_{z}$ affects albedo spectrum?\n",
    "\n",
    "Lastly, we can do the same thing with K_zz, the mixing parameter. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_atmo = pd.read_csv(pj.jupiter_pt(), delim_whitespace=True)\n",
    "\n",
    "all_kzz = [1e6, 1e8, 1e10]\n",
    "w = []\n",
    "a = []\n",
    "all_outs,cld_outs=[],[]\n",
    "for kz in all_kzz:\n",
    "    df_atmo['kz'] = [kz]*df_atmo.shape[0]\n",
    "    sum_planet.atmosphere(df = df_atmo)\n",
    "    cld = sum_planet.virga(['H2O'],directory, fsed=2,mh=metallicity,\n",
    "                 mmw = mean_molecular_weight)\n",
    "    cld_outs += [cld]\n",
    "    out = sum_planet.spectrum(opacity,full_output=True)\n",
    "    x,y = pj.mean_regrid(out['wavenumber'], out['albedo'], R=150)\n",
    "    w += [x]\n",
    "    a += [y]\n",
    "    all_outs += [out['full_output']]\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By looking at the spectra, you might suspect that increasing K_zz simply decreases the optical thickness of the cloud. Let's take a look at what else is changing in the cloud model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f = picplt.spectrum(w,a,plot_width=500, plot_height=300,\n",
    "                     legend=['Kzz= '+str(i) for i in all_kzz])\n",
    "f.legend.location = 'bottom_left'\n",
    "show(f)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From the lowest to the highest case, things don't immediately look very different. Close inspection of the photon attenuation plot reveals that the $\\tau\\sim$1 level of the highest mixing case is indeeed a bit lower than the lower mixing mixing case. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "show(column(picplt.photon_attenuation(all_outs[0],title='kz=1e6',\n",
    "                            plot_width=500, plot_height=300),\n",
    "   picplt.photon_attenuation(all_outs[2],title='kz=1e10',\n",
    "                            plot_width=500, plot_height=300)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lastly, inspecting the particle size distributions is where we see the largest differences. From the lowest, to highest K_zz cases we ran, we see an increase in 4 orders of magnitude in particle size! "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, dndr = cldplt.radii(cld_outs[0],at_pressure=1e-2)\n",
    "print('kz=1e6')\n",
    "show(fig)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, dndr = cldplt.radii(cld_outs[2],at_pressure=1e-1)\n",
    "print('kz=1e10')\n",
    "show(fig)"
   ]
  },
  {
   "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
}