{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# **Reflected Light Phase Curves**\n",
    "\n",
    "**Citation: Hamill et al. (2024) : Reflected Light Phase Curves with PICASO: A Kepler-7b Case Study**\n",
    "\n",
    "Computing a reflected light phase curve is similar to computing 3D reflected light spectra, but now we will compute reflected spectra for an entire grid of phase angles.\n",
    "\n",
    "From the previous tutorials you should understand:\n",
    "\n",
    "1. How to convert GCM input to PICASO's required Xarray\n",
    "2. How to post-process output to append to GCM ouput\n",
    "3. How to run a thermal phase curve and analyze the output\n",
    "\n",
    "In this tutorial, you will learn:\n",
    "\n",
    "1. How to compute a reflected light phase curve\n",
    "2. How to analyze the output"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "\n",
    "from astropy import units as u\n",
    "\n",
    "import xarray as xr\n",
    "\n",
    "from picaso import justdoit as jdi\n",
    "from picaso import justplotit as jpi\n",
    "jpi.output_notebook()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setting Up Reflected Light Phase Curve\n",
    "\n",
    "Reflected light calculations consider only scattering from the dayside hemisphere. The dayside hemisphere will continually come in or out of view depending on the phase angle.\n",
    "\n",
    "Let's load in opacity for the visible to near-infrared and the 3D hot Jupiter atmosperic model:\n",
    "\n",
    "We are resampling our opacities here to make phase curve calculations faster! Use resampling with caution in your own phase curve calculations!!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "opacity = jdi.opannection(wave_range=[0.50,0.90], resample=100)  # this is the wavelength band of Kepler"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gcm_out =  jdi.HJ_pt_3d(as_xarray=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Add chemistry. \n",
    "\n",
    "We can use user-input chemistry for this example.\n",
    "\n",
    "To use post-processed chemical abundances, use the same methods in [Modeling a Phase Curve pt 2 (Robbins-Blanch et al. 2022) : Run Thermal Phase Curve w/ Post-Processed Chemical Abundances.](https://natashabatalha.github.io/picaso/notebooks/9f_PhaseCurves-wChemEq.html)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create coords\n",
    "lon = gcm_out.coords.get('lon').values\n",
    "lat = gcm_out.coords.get('lat').values\n",
    "pres = gcm_out.coords.get('pressure').values\n",
    "\n",
    "fake_chem_H2O = np.random.rand(len(lon), len(lat),len(pres))*0.1+0.1 # create fake data\n",
    "fake_chem_H2 = 1-fake_chem_H2O # create data\n",
    "\n",
    "# put data into a dataset\n",
    "ds_chem = jdi.xr.Dataset(\n",
    "    data_vars=dict(\n",
    "        H2O=([\"lon\", \"lat\",\"pressure\"], fake_chem_H2O,{'units': 'v/v'}),\n",
    "        H2=([\"lon\", \"lat\",\"pressure\"], fake_chem_H2,{'units': 'v/v'}),\n",
    "\n",
    "    ),\n",
    "    coords=dict(\n",
    "        lon=([\"lon\"], lon,{'units': 'degrees'}), #required\n",
    "        lat=([\"lat\"], lat,{'units': 'degrees'}), #required\n",
    "        pressure=([\"pressure\"], pres,{'units': 'bar'})#required*\n",
    "    ),\n",
    "    attrs=dict(description=\"coords with vectors\"),\n",
    ")\n",
    "all_gcm = gcm_out.update(ds_chem)\n",
    "\n",
    "all_gcm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_gcm['temperature'].isel(pressure=34).plot(x='lon',y='lat')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Set up phase curve grid \n",
    "\n",
    "Define the number of phases, and the minimum and maximum phase.\n",
    "\n",
    "For reflected phase curves, the minimum spatial resolution accepted by PICASO is 6 x 6 (num_gangle=6, num_tangle=6)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "case_3d = jdi.inputs()\n",
    "\n",
    "n_phases = 8\n",
    "min_phase = 0\n",
    "max_phase = 2*np.pi\n",
    "phase_grid = np.linspace(min_phase,max_phase,n_phases)#between 0 - 2pi\n",
    "#send params to phase angle routine\n",
    "#case_3d.phase_angle(phase_grid=phase_grid,\n",
    "#                    num_gangle=6, num_tangle=6,calculation='reflected')\n",
    "case_3d.phase_curve_geometry('reflected', phase_grid=phase_grid, num_gangle=6, num_tangle=6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A few notes: \n",
    "\n",
    "- Phase = 0 is the secondary eclipse (maximum reflection)\n",
    "- Phase = 2*pi (~6.28) is the primary eclipse (no reflection).\n",
    "\n",
    "Things to check: \n",
    "- Ensure that (0 degrees longitude, 0 degrees latitude) represents your substellar point on the planet. If it does not, you will need to change your reference point in your GCM.\n",
    "- For the reflected case, zero_point must be set to 'secondary eclipse'.\n",
    "- For a tidally locked planet, set shift to 0 at all phase angles."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "case_3d.atmosphere_4d(all_gcm, shift = np.zeros(n_phases), zero_point='secondary_eclipse',\n",
    "                                     plot=True,verbose=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### `Atmosphere_4d`\n",
    "\n",
    "Atmosphere_4d function works by adding 'phase' as another coordinate in the GCM. 'lat2d' and 'lon2d' are also created and store the specific lats/lons for which the dayside hemisphere is visible to the observer at each individual phase."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "case_3d.inputs['atmosphere']['profile']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compute the reflected light phase curve.\n",
    "\n",
    "Note: Planetary radius, orbital distance, and stellar radius is needed to compute reflected phase curves in units of planet/star flux ratio (Fp/Fs)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "case_3d.gravity(radius=1,radius_unit=jdi.u.Unit('R_jup'),\n",
    "                mass=1, mass_unit=jdi.u.Unit('M_jup')) #any astropy units available\n",
    "case_3d.star(opacity,5000,0,4.0, radius=1, radius_unit=jdi.u.Unit('R_sun'), semi_major=0.06, semi_major_unit=u.Unit('au'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Phase curve computation is below. Computation time will vary widely depending on number of phases and spatial resolution used. We recommend a high performance computer for any phase curves greater than 10x10 spatial resolution. Please see the Appendix of Hamill et al. (2024) for example computation times."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "case_3d.inputs['atmosphere']['profile'].isel(phase=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "allout = case_3d.phase_curve(opacity, n_cpu = 3,#jdi.cpu_count(),\n",
    "                             full_output=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot the phase curve!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "to_plot = 'fpfs_reflected' \n",
    "#collapse = [0.40 ]#micron, could select one or more specific wavelengths instead of the mean \n",
    "collapse='np.mean'\n",
    "\n",
    "phases, all_curves, all_ws, fig=jpi.phase_curve(allout, to_plot,\n",
    "                                                collapse=collapse, R=10,reorder_output=True)\n",
    "\n",
    "jpi.show(fig)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Phase-resolved spectra used to create the phase curve.\n",
    "\n",
    "\n",
    "Since chemistry is randomized for this tutorial and there are no clouds, the spectra is mostly symmetrical across 0 degrees phase angle."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#same old same old\n",
    "wno =[];fpfs=[];legend=[]\n",
    "for iphase in allout.keys():\n",
    "    w,f = jdi.mean_regrid(allout[iphase]['wavenumber'],\n",
    "                               allout[iphase]['fpfs_reflected'], R=100)\n",
    "    wno+=[w]\n",
    "    fpfs+=[f*1e6]\n",
    "    legend +=[str(int(iphase*180/np.pi))]\n",
    "jpi.show(jpi.spectrum(wno, fpfs, plot_width=500,legend=legend,\n",
    "                     palette=jpi.pals.viridis(n_phases)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Reflected phase curve with `Virga` clouds \n",
    "\n",
    "In reflected light brightness can be heavily influenced by cloud layers. Let's add a user-input cloud layer and compare ouputs**\n",
    "\n",
    "To use Virga clouds, please use the same steps found in [Modeling a 3D spectrum (Adams et al. 2022) : Run Cloudy 3D spectra](https://natashabatalha.github.io/picaso/notebooks/9c_PostProcess3Dinput-Clouds.html#Post-Process-Clouds:-virga-cloud-model) or [Modeling a Phase Curve pt 2 : Run Thermal Phase Curves w/ Post-Processed Virga Models](https://natashabatalha.github.io/picaso/notebooks/9f_PhaseCurves-wChemEq.html#Run-Thermal-Phase-Curve-w/-Post-Processed-virga-Models)\n",
    "\n",
    "If you use Virga clouds for reflected phase curves, we highly advise against regridding beforehand. This means your Virga clouds may take a while to compute, but it prevents issues with the re-mapping used in `clouds_4d`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create coords\n",
    "lon = gcm_out.coords.get('lon').values\n",
    "lat = gcm_out.coords.get('lat').values\n",
    "pres = gcm_out.coords.get('pressure').values[:-1]\n",
    "wno_grid = np.linspace(33333,1e4,10) #cloud properties are defined on a wavenumber grid\n",
    "\n",
    "#create box-band cloud model\n",
    "fake_opd = np.zeros((len(lon), len(lat),len(pres), len(wno_grid))) # create fake data\n",
    "where_lat = np.where(((lat>-50) & (lat<50)))#creating a grey cloud band\n",
    "where_lon = np.where(((lon>-90) & (lon<0)))#creating a grey cloud band\n",
    "where_pres = np.where(((pres<0.01) & (pres>.00001)))#creating a grey cloud band, 10 mbar to 0.01 mbar\n",
    "for il in where_lat[0]:\n",
    "    for ilon in where_lon[0]:\n",
    "        for ip in where_pres[0]:\n",
    "            fake_opd[ilon,il,ip,:]=10 #optical depth of 10 (>>1)\n",
    "\n",
    "#make up asymmetry and single scattering properties\n",
    "fake_asymmetry_g0 = 0.5 + np.zeros((len(lon), len(lat),len(pres), len(wno_grid)))\n",
    "fake_ssa_w0 = 0.9 + np.zeros((len(lon), len(lat),len(pres), len(wno_grid)))\n",
    "\n",
    "# put data into a dataset\n",
    "ds_cld= xr.Dataset(\n",
    "    data_vars=dict(\n",
    "        opd=([\"lon\", \"lat\",\"pressure\",\"wno\"], fake_opd,{'units': 'depth per layer'}),\n",
    "        g0=([\"lon\", \"lat\",\"pressure\",\"wno\"], fake_asymmetry_g0,{'units': 'none'}),\n",
    "        w0=([\"lon\", \"lat\",\"pressure\",\"wno\"], fake_ssa_w0,{'units': 'none'}),\n",
    "    ),\n",
    "    coords=dict(\n",
    "        lon=([\"lon\"], lon,{'units': 'degrees'}),#required\n",
    "        lat=([\"lat\"], lat,{'units': 'degrees'}),#required\n",
    "        pressure=([\"pressure\"], pres,{'units': 'bar'}),#required\n",
    "        wno=([\"wno\"], wno_grid,{'units': 'cm^(-1)'})#required for clouds\n",
    "    ),\n",
    "    attrs=dict(description=\"coords with vectors\"),\n",
    ")\n",
    "\n",
    "ds_cld"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Sanity Check \n",
    "\n",
    "Plot of our optically thick cloud layer that lies on the western dayside. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds_cld['opd'].isel(pressure=where_pres[0][0],wno=0).plot(x='lon',y='lat')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Run clouds_4d \n",
    "\n",
    "Set up phase curve calculation same as before. The only difference is we are now running clouds_4d as well."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "case_3d_clouds = jdi.inputs()\n",
    "\n",
    "n_phases = 8\n",
    "min_phase = 0\n",
    "max_phase = 2*np.pi\n",
    "phase_grid = np.linspace(min_phase,max_phase,n_phases)#between 0 - 2pi\n",
    "#send params to phase angle routine\n",
    "#case_3d.phase_angle(phase_grid=phase_grid,\n",
    "#                    num_gangle=6, num_tangle=6,calculation='reflected')\n",
    "\n",
    "case_3d_clouds.phase_curve_geometry('reflected', phase_grid=phase_grid, num_gangle=6, num_tangle=6)\n",
    "case_3d_clouds.inputs['atmosphere']['profile']\n",
    "\n",
    "case_3d_clouds.atmosphere_4d(gcm_out, shift = np.zeros(n_phases), zero_point='secondary_eclipse',\n",
    "                                     plot=True,verbose=False)\n",
    "\n",
    "# run clouds_4d (which must always run after atmosphere_4d)\n",
    "case_3d_clouds.clouds_4d(ds_cld, iz_plot=34, plot=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## **Calculate cloudy phase curve**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "case_3d_clouds.gravity(radius=1,radius_unit=jdi.u.Unit('R_jup'),\n",
    "                mass=1, mass_unit=jdi.u.Unit('M_jup')) #any astropy units available\n",
    "case_3d_clouds.star(opacity,5000,0,4.0, radius=1, radius_unit=jdi.u.Unit('R_sun'), semi_major=0.06, semi_major_unit=u.Unit('au'))\n",
    "\n",
    "allout_clouds = case_3d_clouds.phase_curve(opacity, n_cpu = 3,#jdi.cpu_count(),\n",
    "                             full_output=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "to_plot = 'fpfs_reflected' #or thermal\n",
    "#collapse = [0.40 ]#micron\n",
    "collapse='np.mean'\n",
    "#phases, all_curves, all_ws, fig=jpi.phase_curve(allout, to_plot, collapse=collapse, R=R)\n",
    "phases, all_curves, all_ws, fig=jpi.phase_curve(allout, to_plot, collapse=collapse, R=10,plot_width=1500, reorder_output=True)\n",
    "phases_clouds, all_curves_clouds, all_ws_clouds, fig_clouds=jpi.phase_curve(allout_clouds, to_plot, collapse=collapse, R=10,plot_width=1500, reorder_output=True)\n",
    "#jpi.show(fig)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compare cloudy and cloud free\n",
    "\n",
    "Below is a comparison between our cloudless and cloudy phase curves. The cloudy phase curve is much brighter. We also see a shift in the phase maximum due to our inhomogeneous cloud layer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "plt.figure(figsize=(15, 7.5))\n",
    "\n",
    "plt.plot(phases,all_curves * 1e6, label=\"Cloudless\")\n",
    "plt.plot(phases_clouds,all_curves_clouds * 1e6, label=\"Cloudy\")\n",
    "\n",
    "plt.xlabel('Orbital Phase')\n",
    "plt.ylabel('F$_{P}$/F$_{S}$ (ppm)')\n",
    "plt.title('Reflected Light Phase Curves')\n",
    "plt.legend(loc='upper left')\n",
    "plt.rcParams['font.size'] = 15\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Cloudy phase spectra\n",
    "\n",
    "The spectra for phase angles < 180 is flat due to clouds."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#same old same old\n",
    "wno =[];fpfs=[];legend=[]\n",
    "for iphase in allout_clouds.keys():\n",
    "    w,f = jdi.mean_regrid(allout_clouds[iphase]['wavenumber'],\n",
    "                               allout_clouds[iphase]['fpfs_reflected'], R=100)\n",
    "    wno+=[w]\n",
    "    fpfs+=[f*1e6]\n",
    "    legend +=[str(int(iphase*180/np.pi))]\n",
    "jpi.show(jpi.spectrum(wno, fpfs, plot_width=500,legend=legend,\n",
    "                     palette=jpi.pals.viridis(n_phases)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To further analyze your output, you can use all of the same tools found in [Modeling a Phase Curve pt 1 (Robbins-Blanch et al. 2022) > Phase Curve Plotting Tools: Phase Snaps](https://natashabatalha.github.io/picaso/notebooks/9e_PhaseCurves.html#Phase-Curve-Plotting-Tools:-Phase-Snaps)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "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.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}