{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Spherical Integration in PICASO\n",
    "\n",
    "In this tutorial you will learn the different integration schemes that exist within `PICASO`.  "
   ]
  },
  {
   "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",
    "import picaso.opacity_factory as opa_fa\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": [
    "## Test Geometry with Brown Dwarf Sonora Example"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's try three different cases"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "wave_range = [3,5]\n",
    "\n",
    "opa = jdi.opannection(wave_range=wave_range)\n",
    "\n",
    "case1 = jdi.inputs(calculation='browndwarf')\n",
    "\n",
    "case2 = jdi.inputs(calculation='browndwarf')\n",
    "\n",
    "case3 = jdi.inputs(calculation='browndwarf')\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First case will be a full 10x10 grid, second will be the same 10x10 grid but leveraging symmetry, last will be a 1x10 grid. A 1x10 grid automaticlly leverages symmetry as we will see in the disco plots below. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "names = ['10x10', '10x10sym', '1x10']\n",
    "cases = [ case1, case2, case3]#, case5]\n",
    "\n",
    "case1.phase_angle(num_tangle=10, num_gangle=10)\n",
    "case2.phase_angle(num_tangle=10, num_gangle=10, symmetry=True)\n",
    "case3.phase_angle(num_tangle=1, num_gangle=10) #remember num_tangle=1 automatically insinuates symmetry"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Query from Sonora Profile Grid \n",
    "\n",
    "(Same as before) 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": [
    "# Note here that we do not need to provide case.star, since there is none! \n",
    "for i in cases:i.gravity(gravity=100 , gravity_unit=u.Unit('m/s**2'))\n",
    "    \n",
    "#grab sonora profile \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",
    "for i in cases:i.sonora(sonora_profile_db, Teff)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Run Spectrum"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = {}\n",
    "for i,ikey in zip(cases, names):\n",
    "    print(ikey)\n",
    "    df[ikey] = i.spectrum(opa, full_output=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Convert Units and Regrid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i,ikey in zip(cases, names):\n",
    "    print(ikey)\n",
    "    x,y = df[ikey]['wavenumber'], df[ikey]['thermal'] #units of erg/cm2/s/cm\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[ikey]['fluxnu'] = y \n",
    "    x,y = jdi.mean_regrid(1e4/x, y, R=300) #wavenumber, erg/cm2/s/Hz\n",
    "    df[ikey]['regridy'] =  y\n",
    "    df[ikey]['regridx'] = x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compare Different Integration Schemes with Sonora Grid"
   ]
  },
  {
   "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": {},
   "outputs": [],
   "source": [
    "to_compare = []\n",
    "fraction_compare = []\n",
    "for ikey in names: \n",
    "    to_compare += [df[ikey]['regridy']]\n",
    "    fraction_compare += [df[ikey]['regridy']/sony]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "to_compare += [sony]\n",
    "legend_label = names + ['sonora']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compare Accuracy of Different Number of Integration Points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "show(jpi.spectrum([x]*len(to_compare),to_compare, legend=legend_label\n",
    "                  ,plot_width=800,x_range=wave_range,y_axis_type='log'))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visualize Integration Schemes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for ikey in names:\n",
    "    asdict = df[ikey]['full_output']\n",
    "    print(ikey)\n",
    "    jpi.disco(asdict, calculation='thermal', wavelength=[4.5])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot Gauss Weights \n",
    "\n",
    "Plotting out the gauss weights you can see the difference between the purely 1D (without using chebyshev weights), 2D full grid and 2D symmetric. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = figure(height=300, x_axis_label='gangle',y_axis_label='gweight')\n",
    "for i in [fig.line, fig.circle]: i(cases[0].inputs['disco']['gangle'],\n",
    "         cases[0].inputs['disco']['gweight'],color='red',\n",
    "                                  legend_label=names[0]  )\n",
    "\n",
    "for i in [fig.line, fig.circle]: i(cases[1].inputs['disco']['gangle'],\n",
    "         cases[1].inputs['disco']['gweight'],color='blue',\n",
    "                                 legend_label=names[1] )\n",
    "\n",
    "for i in [fig.line, fig.circle]: i(cases[2].inputs['disco']['gangle'],\n",
    "         cases[2].inputs['disco']['gweight'] \n",
    "                                   ,color='green',legend_label=names[2])\n",
    "    \n",
    "\n",
    "show(fig)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Geometry with Reflectivity and Non-Zero Phase\n",
    "\n",
    "In this example, let's look at four difference cases. \n",
    "\n",
    "1) 1x10 symmetric grid with zero phase (inherently symmetric)\n",
    "\n",
    "2) 10x10 grid with zero phase but let's run without symmetry utilization\n",
    "\n",
    "3) same as 2 but flipping on symmetry \n",
    "\n",
    "4) Non-zero phase, full grid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "wave_range = [0.3,1]\n",
    "\n",
    "opa = jdi.opannection(wave_range=wave_range)\n",
    "\n",
    "case1 = jdi.inputs()\n",
    "case2 = jdi.inputs()\n",
    "case3 = jdi.inputs()\n",
    "case4 = jdi.inputs()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "names = ['1x10','10x10', '10x10sym',  '10x10_60']\n",
    "\n",
    "case1.phase_angle(phase=0, num_tangle=1, num_gangle=10) #remember num_tangle=1 automatically insinuates symmetry\n",
    "case2.phase_angle(phase=0, num_tangle=10, num_gangle=10)\n",
    "case3.phase_angle(phase=0, num_tangle=10, num_gangle=10, symmetry=True)\n",
    "case4.phase_angle(phase=60*np.pi/180, num_tangle=10, num_gangle=10)\n",
    "\n",
    "cases = [ case1, case2, case3, case4]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in cases:i.gravity(gravity=20 , gravity_unit=u.Unit('m/s**2'))\n",
    "for i in cases:i.atmosphere(filename = jdi.jupiter_pt(), delim_whitespace=True)\n",
    "for i in cases:i.approx(raman='pollack')\n",
    "for i in cases:i.star(opa, 5000,0,4.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = {}\n",
    "for i,ikey in zip(cases, names):\n",
    "    print(ikey)\n",
    "    df[ikey] = i.spectrum(opa,calculation='reflected' ,full_output=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "to_compare = []\n",
    "fraction_compare = []\n",
    "for ikey in names: \n",
    "    x,y = df[ikey]['wavenumber'], df[ikey]['albedo']\n",
    "    x,y = jdi.mean_regrid(x, y, R=100) #wavenumber, erg/cm2/s/Hz\n",
    "    df[ikey]['regridy'] =  y\n",
    "    df[ikey]['regridx'] = x\n",
    "    to_compare += [df[ikey]['regridy']]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "show(jpi.spectrum([x]*len(to_compare),to_compare, legend=names\n",
    "                  ,plot_width=800))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you zoom in, you will notice that the 1x10 is just slightly off the full 10x10 grid. This is certainly not enough to warrant doing a full 10x10=100 calculations opposed to the 10/2x1=5 calculations it too to run the 1x10 grid (you can see what RT points were ran below)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for ikey in names:\n",
    "    asdict = df[ikey]['full_output']\n",
    "    print(ikey)\n",
    "    jpi.disco(asdict, calculation='reflected', wavelength=[0.5])"
   ]
  },
  {
   "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"
  },
  "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": 4
}