{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Brown Dwarf Models\n",
    "\n",
    "Modeling brown dwarfs is very similar to modeling thermal emission for exoplanets. The only difference is there is no stellar input!\n",
    "\n",
    "In this tutorial you will learn:\n",
    "\n",
    "1. How to turn that feature off \n",
    "2. Query a profile from the [Sonora Grid](https://zenodo.org/record/1309035#.Xo5GbZNKjGJ). Note, this is note necessary -- just convenient! \n",
    "3. Create a Brown Dwarf Spectrum "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import astropy.units as u\n",
    "\n",
    "#picaso\n",
    "from picaso import justdoit as jdi \n",
    "from picaso import justplotit as jpi\n",
    "#plotting\n",
    "from bokeh.io import output_notebook\n",
    "\n",
    "output_notebook()\n",
    "from bokeh.plotting import show,figure"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Start with the same inputs as before"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "wave_range = [3,5]\n",
    "\n",
    "opa = jdi.opannection(wave_range=wave_range)\n",
    "\n",
    "bd = jdi.inputs(calculation='browndwarf')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Note here that we do not need to provide case.star, since there is none! \n",
    "bd.gravity(gravity=100 , gravity_unit=u.Unit('m/s**2'))\n",
    "\n",
    "#this is the integration setup that was used to compute the Sonora grid \n",
    "#take a look at Spherical Integration Tutorial to get a look at what these parameters do\n",
    "bd.phase_angle(0)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Download and Query from Sonora Profile Grid \n",
    "\n",
    "Download the profile files that are located in the [profile.tar file](https://zenodo.org/record/1309035#.Xo5GbZNKjGJ)\n",
    "\n",
    "Once you untar the file you can set the file path below. You do not need to unzip each profile. `picaso` will do that upon read in. `picaso` will find the nearest neighbor and attach it to your class. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#point to where you untared your sonora profiles\n",
    "sonora_profile_db = '/data/sonora_profile/'\n",
    "Teff = 900 \n",
    "#this function will grab the gravity you have input above and find the nearest neighbor with the \n",
    "#note sonora chemistry grid is on the same opacity grid as our opacities (1060). \n",
    "bd.sonora(sonora_profile_db, Teff)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note if you have added anything to your atmosphere profile (`bd.inputs['atmosphere']['profile']`), `sonora` function will overwrite it! Therefore, **make sure that you add new column fields after running sonora.** "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Run Spectrum"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = bd.spectrum(opa, full_output=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Convert to $F_\\nu$ Units and Regrid\n",
    "\n",
    "`PICASO` outputs the raw flux as: \n",
    "\n",
    "$$ F_\\lambda ( \\frac{erc}{cm^2 * sec * cm}) $$\n",
    "\n",
    "Typical fluxes shown in several Brown Dwarf papers are: \n",
    "\n",
    "$$ F_\\nu ( \\frac{erc}{cm^2 * sec * Hz}) $$\n",
    "\n",
    "Below is a little example of how to convert units. \n",
    "\n",
    "**NOTE**: Some people like to plot out `Eddington Flux`, $H_\\nu$. This gets confusing as the units appear to be erg/cm2/s/Hz but you will notice a factor of four difference: \n",
    "\n",
    "$$H_\\nu = \\frac{F_\\nu}{4}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x,y = df['wavenumber'], df['thermal'] #units of erg/cm2/s/cm\n",
    "\n",
    "xmicron = 1e4/x\n",
    "\n",
    "flamy = y*1e-8 #per anstrom instead of per cm\n",
    "sp = jdi.psyn.ArraySpectrum(xmicron, flamy, \n",
    "                            waveunits='um', \n",
    "                            fluxunits='FLAM')   \n",
    "sp.convert(\"um\")\n",
    "sp.convert('Fnu') #erg/cm2/s/Hz\n",
    "\n",
    "x = sp.wave #micron\n",
    "y= sp.flux #erg/cm2/s/Hz\n",
    "df['fluxnu'] = y \n",
    "x,y = jdi.mean_regrid(1e4/x, y, R=300) #wavenumber, erg/cm2/s/Hz\n",
    "df['regridy'] =  y\n",
    "df['regridx'] = x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compare with Sonora Grid\n",
    "\n",
    "The corresponding spectra are also available at the same link above. `PICASO` doesn't provide query functions\n",
    "for this. So if you want to compare, you will have to read in the files as is done below"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "son = pd.read_csv('sp_t900g100nc_m0.0',delim_whitespace=True, \n",
    "                 skiprows=3,header=None,names=['w','f'])\n",
    "sonx, sony =  jdi.mean_regrid(1e4/son['w'], son['f'], newx=x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "show(jpi.spectrum([x]*2,[df['regridy'], sony], legend=['PICASO', 'Sonora']\n",
    "                  ,plot_width=800,x_range=wave_range,y_axis_type='log'))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Thermal Contribution Function \n",
    "\n",
    "Contribution functions give us an understanding of where pressures flux is being emitted (e.g. [Lothringer et al 2018](https://iopscience.iop.org/article/10.3847/1538-4357/aadd9e#apjaadd9es3-3-1) Figure 12) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax, CF = jpi.thermal_contribution(df['full_output'], norm=jpi.colors.LogNorm(vmin=1e7, vmax=1e11))"
   ]
  },
  {
   "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.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
}