Interactives/Exoplanet_Sim.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exoplanet System Simulation Interactive\n",
    "\n",
    "This interactive application simulates the orbit of of a single exoplanet around a star, taking into account\n",
    "\n",
    "- the masses of the planet and the star, and\n",
    "- the shape of and orientation of the orbit in space\n",
    "\n",
    "to accurately model the orbit as well as the observable radial velocites and light curve of the exoplanetary system as observed from Earth.  This interactively crudely makes two estimates:\n",
    "\n",
    "- The temperature of the planet is estimated by assuming it reflect light like Earth does and that its temperature doesn't change during its orbit (so it uses the average distance from the star for temperature calculations). \n",
    "- The radius of the planet is estimated by fitting the 636 exoplanets with both estimated masses and radii (obtained from the NASA Exoplanet Archive on 15 August 2018).  It assumes no physics, and so is probably a very rough estimate.\n",
    "\n",
    "**Note**: This code takes an approach originally seen in the *TwoStars* code by Carroll and Ostlie.  That code was used for modelling binary stars and has been extensively changed to handle exoplanets and optimized for the Python programming language.\n",
    "\n",
    "## Model versus Reality \n",
    "\n",
    "A reminder that the model shown on the left hand side of the exoplanet orbiting the star is supposed to be just that, a model.  Remember we rarely actually visually see the exoplanets in system like this as they are in reality billions of times less luminous than the star.  What we actually observe is a radial velocity or light curve, this model is an explanation for those observations. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Developed by Juan Cabanela starting June 19, 2018 and proceeding through the Summer of 2018\n",
    "#\n",
    "# This simulation is a modified version of the Binary Star Simulation Interactive reworked \n",
    "# to handle exoplanetary systems."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.display import display, HTML\n",
    "import numpy as np\n",
    "import ipywidgets as widgets\n",
    "import traitlets\n",
    "import pythreejs as p3j\n",
    "import bqplot as bq\n",
    "import tempNcolor as tc\n",
    "import starlib as star\n",
    "import number_formatting as nf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## FUNCTIONS ##\n",
    "\n",
    "def property_update(change=None):\n",
    "    '''\n",
    "    This function updates the stellar properties (and orbital properties) and is meant to be used\n",
    "    when there are changes to stellar properties as controled by various ipywidgets on screen.\n",
    "    '''\n",
    "    global bsm\n",
    "\n",
    "    # Turn off continuous updating to allow quick changing of several variables\n",
    "    # NOTE: Masses automatically update in the binary star model.\n",
    "    bsm.continuous_update = False  \n",
    "    (starradius.value, startemp.value, hexcolor1) = star.ConfigStar(starmass_slider.value)\n",
    "    (exoradius.value, exotemp.value, hexcolor2) = ConfigPlanet(exomass_slider.value, \n",
    "                                                                      semimajor_slider.value*AU2RSun, \n",
    "                                                                      startemp.value, \n",
    "                                                                      starradius.value)\n",
    "    # Copy radius and temperature info to binary star model\n",
    "    bsm.rad1 = starradius.value\n",
    "    bsm.rad2 = exoradius_in_RSun.value\n",
    "    bsm.temp1 = startemp.value\n",
    "    bsm.temp2 = exotemp.value\n",
    "    \n",
    "    # Turn on continuous updating again\n",
    "    bsm.continuous_update = True \n",
    "\n",
    "    # Force a binary star model update (which should update view as well)\n",
    "    bsm.force_update()\n",
    "    \n",
    "    # If there is a collision, stop the phase changes, otherwise allow access to changing phase\n",
    "    if (bsm.collision):\n",
    "        phase_title.value = phase_title_collision\n",
    "        phase_play.step = 0\n",
    "        phase_play.disabled = True\n",
    "        phase_slider.disabled = True\n",
    "    else:\n",
    "        # Make sure orbital phase can be adjusted\n",
    "        phase_title.value = phase_title_default\n",
    "        phase_play.step = 1\n",
    "        phase_play.disabled = False\n",
    "        phase_slider.disabled = False\n",
    "\n",
    "    # Set luminosities\n",
    "    L1_output.value = str(nf.SigFig((startemp.value/star.Te_Sun)**4 * starradius.value**2, 3))\n",
    "    L2_output.value = str(nf.SigFig((exotemp.value/star.Te_Sun)**4 * exoradius.value**2, 3))\n",
    "    \n",
    "    # Update the on screen label to indicate unreal scaling if necessary\n",
    "    if (exosystem_view.multiplier > 1):\n",
    "        sys_title.value = sys_title_scaling\n",
    "    else:\n",
    "        sys_title.value=sys_title_default\n",
    "        \n",
    "    # Revise orbital readouts\n",
    "    P_output.value = \"{0:.1f}\".format(bsm.P)\n",
    "    ap_output.value = \"{0:.1f}\".format(bsm.ap*AU2RSun)\n",
    "    aa_output.value = \"{0:.1f}\".format(bsm.aa*AU2RSun)\n",
    "    \n",
    "    # Perform the inclination update to force an update of the display\n",
    "    inclination_update()\n",
    "\n",
    "    \n",
    "def inclination_update(change=None):\n",
    "    '''\n",
    "    This function updates the star system's inclination to the plane of the sky and updates the graph shown\n",
    "    '''\n",
    "    global bsm\n",
    "    \n",
    "    # Retrieve updated radial velocity or light curves and then update figure\n",
    "    # (in OO version of this code, the bsm should have updated automatically following the inclination angle change)\n",
    "    if (fig_selector.value == rv_val):\n",
    "        update_radvel_curve(bsm.radvel_info, bsm.orbit_info)\n",
    "    else:\n",
    "        update_light_curve(bsm.lc_info)\n",
    "    \n",
    "    \n",
    "def position_update(change=None):\n",
    "    '''\n",
    "    This function updates the two stars' phase line on whatever graph is shown.  The viewer changes are handled\n",
    "    automatically by linking the slider to the t_idx value in the exosystem_view object.\n",
    "    '''\n",
    "    # Update the phase line in radial velocity or light curve\n",
    "    if (fig_selector.value == rv_val):\n",
    "        rv_phase_line.x = [bsm.radvel_info['phase'][phase_slider.value], bsm.radvel_info['phase'][phase_slider.value]]\n",
    "    else:\n",
    "        lc_phase_line.x = [bsm.lc_info['phase'][phase_slider.value], bsm.lc_info['phase'][phase_slider.value]]\n",
    "\n",
    "        \n",
    "def scene_update(change=None):\n",
    "    '''\n",
    "    This function updates the notificaiton of the scaling of the objects and the grid/orbit drawing if necessary.\n",
    "    '''\n",
    "    if (exosystem_view.multiplier > 1):\n",
    "        sys_title.value = sys_title_scaling\n",
    "    else:\n",
    "        sys_title.value=sys_title_default\n",
    "    \n",
    "    # Reset the grid and orbit\n",
    "    exosystem_view.reset_grid_and_orbit()\n",
    "\n",
    "def zoom_update(change=None):\n",
    "    '''\n",
    "    This function handles updates to the level of magnification.\n",
    "    '''\n",
    "    exosystem_view._starcam.zoom=zoom_slider.value\n",
    "        \n",
    "def graph_update(change=None):\n",
    "    '''\n",
    "    This function switches which graph to plot\n",
    "    '''\n",
    "    global bsm, graph_fig\n",
    "    \n",
    "    # Check if graph exists, if it doesn't, create it.\n",
    "    try:\n",
    "        graph_fig\n",
    "    except NameError:\n",
    "        graph_fig = None\n",
    "    \n",
    "    if (fig_selector.value == rv_val):\n",
    "        bsm.wipe_lc_info()\n",
    "        bsm.set_radvel_info()\n",
    "        new_fig = create_radvel_curve(bsm.radvel_info, bsm.orbit_info)\n",
    "    else:\n",
    "        bsm.wipe_radvel_info()\n",
    "        bsm.set_lc_info()\n",
    "        new_fig = create_light_curve(bsm.lc_info)\n",
    "    \n",
    "    if (graph_fig is None):\n",
    "        graph_fig = new_fig\n",
    "    else:\n",
    "        graph_fig.marks = new_fig.marks\n",
    "        graph_fig.axes = new_fig.axes\n",
    "        graph_fig.title = new_fig.title\n",
    "        graph_fig.layout = new_fig.layout\n",
    "\n",
    "\n",
    "def create_light_curve(light_curve):\n",
    "    '''\n",
    "    Initialize the entire light curve \n",
    "    '''\n",
    "    global lc_line, lc_phase_line\n",
    "    \n",
    "    # Determine deepest possible eclipse (90 degree inclination)\n",
    "    edgeon_lc = star.LightCurveInfo(bsm.orbit_info, 90, bsm.rad1, bsm.rad2, bsm.temp1, bsm.temp2)\n",
    "    min_flux = np.floor(np.min(edgeon_lc['F_norm'])*100)/100   # Round minimum to nearest percent\n",
    "    \n",
    "    # Set scales\n",
    "    sc_x = bq.LinearScale(min=0, max=1)\n",
    "    sc_y = bq.LinearScale(min=min_flux, max=1)\n",
    "\n",
    "    # Build the light curve\n",
    "    lc_line = bq.Lines(x=light_curve['phase'], y=light_curve['F_norm'], scales={'x': sc_x, 'y': sc_y}, \n",
    "                      colors=['Black'])\n",
    "    \n",
    "    # Indicate the current phase\n",
    "    x_phase = [light_curve['phase'][phase_slider.value], light_curve['phase'][phase_slider.value]]\n",
    "    y_phase = [min_flux, 1] \n",
    "    lc_phase_line = bq.Lines(x=x_phase, y=y_phase, scales={'x': sc_x, 'y': sc_y}, \n",
    "                      colors=['Red'])\n",
    "    \n",
    "    # Setup axes and return figure\n",
    "    ax_x = bq.Axis(scale=sc_x, label='Phase')\n",
    "    ax_y = bq.Axis(scale=sc_y, orientation='vertical', label='Fraction of Maximum Flux')\n",
    "    return bq.Figure(marks=[lc_line, lc_phase_line], axes=[ax_x, ax_y], title='Light Curve',\n",
    "                     layout=widgets.Layout(width=graph_width, height=graph_height, margin='5px 5px 5px 5px'))\n",
    "\n",
    "\n",
    "def update_light_curve(lc_info):\n",
    "    '''\n",
    "    Update the light curve in event of inclination change\n",
    "    '''\n",
    "    global lc_line\n",
    "    lc_line.x=lc_info['phase']\n",
    "    lc_line.y=lc_info['F_norm']\n",
    "\n",
    "    \n",
    "def create_radvel_curve(radvel_info, orbit_info):\n",
    "    '''\n",
    "    Initialize the radial velocity curve\n",
    "    '''\n",
    "    global star_line, rv_phase_line\n",
    "    \n",
    "    # Set up the scales\n",
    "    sc_x = bq.LinearScale()\n",
    "    sc_y = bq.LinearScale()\n",
    "\n",
    "    # Indicate the current phase\n",
    "    x_phase = [radvel_info['phase'][phase_slider.value], radvel_info['phase'][phase_slider.value]]\n",
    "    \n",
    "    # Scale to min/max velocity at all inclinations (assume no systemic radial velocity)\n",
    "    maxval = -np.min(orbit_info['vx1']*1000)  # use m/s\n",
    "    minval = -np.max(orbit_info['vx1']*1000)  # use m/s\n",
    "    y_phase = [minval, maxval] \n",
    "    rv_phase_line = bq.Lines(x=x_phase, y=y_phase, scales={'x': sc_x, 'y': sc_y}, \n",
    "                      colors=['Red'])\n",
    "    \n",
    "    # Draw the radial velocity curves (in m/s)\n",
    "    star_line = bq.Lines(x=radvel_info['phase'], y=radvel_info['v1r']*1000, scales={'x': sc_x, 'y': sc_y},\n",
    "                         colors=['DarkOrange'], labels=['Star'], display_legend=True)\n",
    "    \n",
    "    # Setup axes and return (initially invisible) figure\n",
    "    ax_x = bq.Axis(scale=sc_x, label='Phase')\n",
    "    ax_y = bq.Axis(scale=sc_y, orientation='vertical', label='Radial velocity (m/s)')\n",
    "    ax_y.label_offset = '3.5em'\n",
    "    return bq.Figure(marks=[star_line, rv_phase_line], axes=[ax_x, ax_y], title='Radial Velocity Curve',\n",
    "                     layout=widgets.Layout(width=graph_width, height=graph_height, margin='5px 5px 5px 5px'))\n",
    "\n",
    "\n",
    "def update_radvel_curve(radvel_info, orbit_info):\n",
    "    '''\n",
    "    Update the radial velocity curve (only stellar curve would be visible)\n",
    "    '''\n",
    "    global star_line, rv_phase_line\n",
    "    \n",
    "    star_line.x=radvel_info['phase']\n",
    "    star_line.y=radvel_info['v1r']*1000   # use m/s\n",
    "    \n",
    "    # Rescale phase line limits if orbit changes\n",
    "    maxval = -np.min(orbit_info['vx1']*1000)  # use m/s\n",
    "    minval = -np.max(orbit_info['vx1']*1000)  # use m/s\n",
    "    rv_phase_line.y = [minval, maxval] \n",
    "\n",
    "\n",
    "def ConfigPlanet(M_exo, semimajor, T_star, R_star, albedo = 0.3):\n",
    "    \"\"\"\n",
    "    Configure planet parameters.\n",
    "    \n",
    "    Requires R_star and semimajor be in the SAME units.\n",
    "    Requires M_exo be in Jupiter masses.\n",
    "    \n",
    "    Radius is based on a fit of exoplanet mass to radius, it is very crude.\n",
    "    \n",
    "    For temperature we use a crude non-greenhouse model that assumes \n",
    "    average distance from the star (semimajor axis) sets the temperature.\n",
    "    \"\"\"\n",
    "    # Estimate the radius of the planet in Jupiter Radii given mass in Jupiter Masses\n",
    "    # Limit the masses to range from 0.01 to 20 Jupiter masses to stay within fit limits.\n",
    "    MvsR_coeff = [0.003216146262957223, 0.02517307913211346, 0.042740981066343345,\n",
    "                  -0.10205588511119912, -0.29830452706869, 0.3442309265319224, \n",
    "                  0.04882412681603183]\n",
    "    logM_vs_logR_fit = np.poly1d(MvsR_coeff)\n",
    "    logR_est = logM_vs_logR_fit(np.log10(M_exo))\n",
    "    R_est = pow(10, logR_est)\n",
    "    \n",
    "    # Set color of planet to be 'cool star' like\n",
    "    planet_colortemp = 4000\n",
    "    hexcolor = tc.rgb2hex(tc.temp2rgb(planet_colortemp))[0]\n",
    "    \n",
    "    # Compute temperature\n",
    "    # Works as long as R_star and semimajor axis are in the same units\n",
    "    T_exo = pow((1 - albedo), 0.25)*pow(R_star/(2*semimajor), 0.5)*T_star\n",
    "    \n",
    "    return (R_est, T_exo, hexcolor)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## INTERACTIVE/DISPLAY WIDGETS ##\n",
    "\n",
    "# Conversion constants\n",
    "R_Jupiter = 6.9911e7  # Radius of Jupiter in km\n",
    "M_Jupiter = 1.8982e27   # Mass of Jupiter in kg\n",
    "RJup_in_RSun = R_Jupiter/star.R_Sun\n",
    "MJup_in_MSun = M_Jupiter/star.M_Sun\n",
    "RSun2AU = star.R_Sun/star.AU\n",
    "AU2RSun = 1/RSun2AU\n",
    "\n",
    "# Define constants\n",
    "min_mass = 0.1   # Minimum stellar mass in solar masses\n",
    "max_mass = 5    # Maximum stellar mass in solar masses\n",
    "mass_step = 0.1  # Step size for stellar mass slider in solar masses\n",
    "init_mass = 1    # Initial mass of star in solar masses\n",
    "\n",
    "min_exomass = 0.1   # Minimum planetary mass in Jupiter masses\n",
    "max_exomass = 5    # Maximum planetary mass in Jupiter masses\n",
    "exomass_step = 0.1  # Step size for mass sliders in Jupiter masses\n",
    "init_exomass = 1    # Initial mass of planet in Jupiter masses\n",
    "\n",
    "min_radius = 0.2   # Minimum stellar radius in solar radii\n",
    "max_radius = 100   # Maximum stellar radius in solar radii\n",
    "radius_step = 0.1  # Step size for stellar radius slider in solar radii\n",
    "init_radius = 1    # Initial radius of star in solar radii\n",
    "\n",
    "min_exoradius = 0.05   # Minimum planetary radius in Jupiter radii\n",
    "max_exoradius = 20   # Maximum planetary radius in Jupiter radii\n",
    "exoradius_step = 0.05  # Step size for planetary radius slider in Jupiter radii\n",
    "init_exoradius = 1    # Initial radius of planet in Jupiter radii\n",
    "\n",
    "min_temp = 3100    # Minimum stellar temperature in K\n",
    "max_temp = 40000   # Maximum stellar temperature in K\n",
    "temp_step = 10     # Step size for temperature slider in K\n",
    "init_temp = int(star.Te_Sun/10)*10 # Initial temperature of star in K\n",
    "\n",
    "min_a = 0.02      # Minimum semimajor axis of system in AU\n",
    "max_a = 5         # Maximum semimajor axis of system in AU\n",
    "step_a = 0.02     # Step size for separation slider in AU\n",
    "init_a = 0.3      # Start off with the two objects this close together\n",
    "\n",
    "init_incl = 0.0   # Initial inclination value (non-zero value avoids odd orientation of FOV at start)\n",
    "init_phi = 0      # Initial semimajor axis phase angle\n",
    "init_ecc = 0.2    # Initial orbital eccentricity\n",
    "\n",
    "view_factor = 1.25  # How many times the maximum distance to place the viewer\n",
    "N = 3000          # Number of time steps to use for orbit\n",
    "Na = 100          # Number of annuli to break up stars into for computing eclipse fraction\n",
    "Ntheta = 360      # Number of angular steps to break up stars into for computing eclipse fraction\n",
    "\n",
    "#\n",
    "# Define some widths to use throughout for layout of controls\n",
    "#\n",
    "\n",
    "# Set simulation size\n",
    "view_width = 350\n",
    "view_height = view_width\n",
    "\n",
    "# Initialize slider sizes\n",
    "EntireWidth = '1000px'\n",
    "SimWidth = '{0:.0f}px'.format(view_width)\n",
    "ControlColWidth = '450px'\n",
    "slider_width = '300px'\n",
    "slider_minwidth = '250px'\n",
    "readout_width = '70px'\n",
    "lum_width = '120px'\n",
    "inform_width = '200px'\n",
    "graph_width = '450px'\n",
    "graph_height = '300px'\n",
    "\n",
    "##\n",
    "## Create control for selecting Figure\n",
    "##\n",
    "rv_val = 'Radial Velocity Curve'\n",
    "lc_val = 'Light Curve'\n",
    "fig_selector = widgets.ToggleButtons(options=[rv_val, lc_val],\n",
    "                                    value=rv_val,\n",
    "                                    #description='Plot to Display:',\n",
    "                                    #style = {'description_width': 'initial'},\n",
    "                                    disabled=False,\n",
    "                                    orientation='horizontal',\n",
    "                                    layout=widgets.Layout(width=ControlColWidth,\n",
    "                                                          height='70px', \n",
    "                                                          overflow='visible')\n",
    "                                   )\n",
    "\n",
    "\n",
    "\n",
    "##\n",
    "##Create controls for stellar parameters \n",
    "##\n",
    "\n",
    "## Mass\n",
    "\n",
    "starmass_slider = widgets.FloatSlider(\n",
    "    value=init_mass,\n",
    "    min=min_mass,\n",
    "    max=max_mass+(mass_step/2),\n",
    "    step=mass_step,\n",
    "    description=\"Star mass\",\n",
    "    style = {'description_width': 'initial'},\n",
    "    disabled=False,\n",
    "    continuous_update=False,\n",
    "    orientation='horizontal',\n",
    "    readout=False,\n",
    "    readout_format='.1f',\n",
    "    layout=widgets.Layout(width=slider_width, min_width=slider_minwidth, \n",
    "                          overflow='visible')\n",
    ")\n",
    "\n",
    "exomass_slider = widgets.FloatSlider(\n",
    "    value=init_mass,\n",
    "    min=min_exomass,\n",
    "    max=max_exomass+(exomass_step/2),\n",
    "    step=exomass_step,\n",
    "    description=\"Planet mass\",\n",
    "    style = {'description_width': 'initial'},\n",
    "    disabled=False,\n",
    "    continuous_update=False,\n",
    "    orientation='horizontal',\n",
    "    readout=False,\n",
    "    readout_format='.1f',\n",
    "    layout=widgets.Layout(width=slider_width, min_width=slider_minwidth, \n",
    "                          overflow='visible')\n",
    ")\n",
    "\n",
    "# Define text boxes for readout\n",
    "starmass_readout = widgets.BoundedFloatText(min=starmass_slider.min, max=starmass_slider.max, \n",
    "                                         value=starmass_slider.value, \n",
    "                                         layout=widgets.Layout(width=readout_width, \n",
    "                                                               overflow='visible'))\n",
    "exomass_readout = widgets.BoundedFloatText(min=exomass_slider.min, max=exomass_slider.max, \n",
    "                                         value=exomass_slider.value, \n",
    "                                         layout=widgets.Layout(width=readout_width,\n",
    "                                                               overflow='visible'))\n",
    "\n",
    "# Link slider and textboxes\n",
    "widgets.jslink((starmass_readout, 'value'), (starmass_slider, 'value'))\n",
    "widgets.jslink((exomass_readout, 'value'), (exomass_slider, 'value'))\n",
    "\n",
    "# Create the individual controls for stellar masses\n",
    "solar_mass = widgets.HTML('M<sub>&#x2609;</sub>')\n",
    "jupiter_mass = widgets.HTML('M<sub>Jup</sub>')\n",
    "starmass_cntl = widgets.HBox([starmass_slider, starmass_readout, solar_mass], \n",
    "                          layout=widgets.Layout(width=ControlColWidth, \n",
    "                                                overflow='visible'))\n",
    "exomass_cntl = widgets.HBox([exomass_slider, exomass_readout, jupiter_mass], \n",
    "                          layout=widgets.Layout(width=ControlColWidth, \n",
    "                                                overflow='visible'))\n",
    "\n",
    "# Define planetary parameters in solar units as globals\n",
    "exomass_in_MSun = widgets.FloatText(value=exomass_slider.value*MJup_in_MSun)\n",
    "def ConvertMass(change):\n",
    "    global exomass_in_MSun\n",
    "    exomass_in_MSun.value = exomass_slider.value*MJup_in_MSun\n",
    "\n",
    "# Automatically update the values in solar units when values for planet change\n",
    "exomass_slider.observe(ConvertMass, names=['value'])\n",
    "\n",
    "\n",
    "## Radius\n",
    "\n",
    "# Define text boxes for readout\n",
    "starradius = widgets.FloatText(value=init_radius,\n",
    "                               disabled=True,\n",
    "                               layout=widgets.Layout(width=readout_width,\n",
    "                                                     overflow='visible'))\n",
    "exoradius = widgets.FloatText(value=init_exoradius,\n",
    "                              disabled=True,\n",
    "                              layout=widgets.Layout(width=readout_width,\n",
    "                                                    overflow='visible'))\n",
    "\n",
    "# Create the individual controls for stellar masses\n",
    "radius_label = widgets.HTML('Radius: ')\n",
    "solar_radius = widgets.HTML('R<sub>&#x2609;</sub>')\n",
    "jupiter_radius = widgets.HTML('R<sub>Jup</sub>')\n",
    "\n",
    "# Define planetary parameters in solar units as globals\n",
    "exoradius_in_RSun = widgets.FloatText(value=exoradius.value*RJup_in_RSun)\n",
    "\n",
    "def ConvertRadius(change):\n",
    "    global exoradius_in_RSun\n",
    "    exoradius_in_RSun.value = exoradius.value*RJup_in_RSun\n",
    "\n",
    "# Automatically update the values in solar units when values for planet change\n",
    "exoradius.observe(ConvertRadius, names=['value'])\n",
    "\n",
    "## Temperature\n",
    "\n",
    "startemp = widgets.IntText(value=init_temp,\n",
    "                           disabled=True,\n",
    "                           layout=widgets.Layout(width=readout_width,\n",
    "                                                 overflow='visible'))\n",
    "exotemp = widgets.IntText(value=init_temp,\n",
    "                          disabled=True,\n",
    "                          layout=widgets.Layout(width=readout_width,\n",
    "                                                overflow='visible'))\n",
    "\n",
    "# Create the individual controls for stellar masses\n",
    "temp_label = widgets.HTML('Temp.: ')\n",
    "Kelvin = widgets.HTML('K')\n",
    "\n",
    "## Luminosity\n",
    "lum_label = widgets.HTML('Lum.:')\n",
    "solar_lum = widgets.HTML('L<sub>&#x2609;</sub>')\n",
    "\n",
    "L1_output = widgets.Text(value = str(0),\n",
    "                         style = {'description_width': 'initial'},\n",
    "                         disabled = True, \n",
    "                         layout=widgets.Layout(width=readout_width, \n",
    "                                               overflow='visible'))\n",
    "\n",
    "L2_output = widgets.Text(value = str(0),\n",
    "                         style = {'description_width': 'initial'},\n",
    "                         disabled = True, \n",
    "                         layout=widgets.Layout(width=readout_width, \n",
    "                                               overflow='visible'))\n",
    "\n",
    "## Information strips on star and planet\n",
    "gap = widgets.HTML('&nbsp;', layout=widgets.Layout(width='10px', overflow='visible'))\n",
    "star_info = widgets.HBox([radius_label, starradius, solar_radius, gap,\n",
    "                          temp_label, startemp, Kelvin, gap,\n",
    "                          lum_label, L1_output, solar_lum],\n",
    "                          layout=widgets.Layout(width=ControlColWidth, \n",
    "                                                overflow='visible'))\n",
    "exo_info = widgets.HBox([radius_label, exoradius, jupiter_radius, gap,\n",
    "                         temp_label, exotemp, Kelvin, gap,\n",
    "                         lum_label, L2_output, solar_lum],\n",
    "                        layout=widgets.Layout(width=ControlColWidth,\n",
    "                                              overflow='visible'))\n",
    "\n",
    "##\n",
    "##Create controls for system properties\n",
    "##\n",
    "\n",
    "# These sliders change entire orbital model and should NOT be continously updated\n",
    "semimajor_slider = widgets.FloatSlider(value=init_a, \n",
    "                                       min=min_a, \n",
    "                                       max=max_a,\n",
    "                                       step=step_a,\n",
    "                                       description=\"Semimajor Axis\",\n",
    "                                       style = {'description_width': 'initial'},\n",
    "                                       disabled=False,\n",
    "                                       continuous_update=False,\n",
    "                                       orientation='horizontal',\n",
    "                                       readout=False,\n",
    "                                       readout_format='.0f',\n",
    "                                       layout=widgets.Layout(width=slider_width, min_width=slider_minwidth,\n",
    "                                                             overflow='visible') )\n",
    "\n",
    "ecc_slider = widgets.FloatSlider(value=init_ecc,\n",
    "                                 min=0,\n",
    "                                 max=0.8,\n",
    "                                 step=0.02,\n",
    "                                     description=\"Eccentricity\",\n",
    "                                 style = {'description_width': 'initial'},\n",
    "                                 disabled=False,\n",
    "                                 continuous_update=False,\n",
    "                                     orientation='horizontal',\n",
    "                                     readout=False,\n",
    "                                 readout_format='.2f',\n",
    "                                 layout=widgets.Layout(width=slider_width, min_width=slider_minwidth, \n",
    "                                                       overflow='visible') )\n",
    "\n",
    "phi_slider = widgets.FloatSlider(value=init_phi,\n",
    "                                 min=0,\n",
    "                                 max=180,\n",
    "                                 step=1,\n",
    "                                 description=\"Major Axis Longitude\",\n",
    "                                 style = {'description_width': 'initial'},\n",
    "                                 disabled=False,\n",
    "                                 continuous_update=False,\n",
    "                                 orientation='horizontal',\n",
    "                                     readout=False,\n",
    "                                 readout_format='.0f',\n",
    "                                 layout=widgets.Layout(width=slider_width, min_width=slider_minwidth,\n",
    "                                                       overflow='visible') )\n",
    "\n",
    "incl_slider = widgets.FloatSlider(value=init_incl,\n",
    "                                  min=0,\n",
    "                                  max=90,\n",
    "                                  step=1,\n",
    "                                  description=\"Inclination\",\n",
    "                                  style = {'description_width': 'initial'},\n",
    "                                  disabled=False,\n",
    "                                  continuous_update=False,\n",
    "                                  orientation='horizontal',\n",
    "                                  readout=False,\n",
    "                                  readout_format='.0f',\n",
    "                                  layout=widgets.Layout(width=slider_width, min_width=slider_minwidth, \n",
    "                                                        overflow='visible') )\n",
    "\n",
    "# Define text boxes for readout\n",
    "semimajor_readout = widgets.BoundedFloatText(min=semimajor_slider.min, max=semimajor_slider.max, \n",
    "                                             value=semimajor_slider.value, \n",
    "                                             layout=widgets.Layout(width=readout_width, \n",
    "                                                                   overflowy='visible'))\n",
    "ecc_readout = widgets.BoundedFloatText(min=ecc_slider.min, max=ecc_slider.max, \n",
    "                                       value=ecc_slider.value, \n",
    "                                       layout=widgets.Layout(width=readout_width, \n",
    "                                                             overflow='visible'))\n",
    "incl_readout = widgets.BoundedFloatText(min=incl_slider.min, max=incl_slider.max, \n",
    "                                        value=incl_slider.value, \n",
    "                                        layout=widgets.Layout(width=readout_width, \n",
    "                                                              overflow='visible'))\n",
    "phi_readout = widgets.BoundedFloatText(min=phi_slider.min, max=phi_slider.max, \n",
    "                                       value=phi_slider.value, \n",
    "                                       layout=widgets.Layout(width=readout_width, \n",
    "                                                             overflow='visible'))\n",
    "# Link slider and textboxes\n",
    "widgets.jslink((semimajor_readout, 'value'), (semimajor_slider, 'value'))\n",
    "widgets.jslink((ecc_readout, 'value'), (ecc_slider, 'value'))\n",
    "widgets.jslink((incl_readout, 'value'), (incl_slider, 'value'))\n",
    "widgets.jslink((phi_readout, 'value'), (phi_slider, 'value'))\n",
    "\n",
    "# Create the individual controls for system properties\n",
    "Solar_radius = widgets.HTML('R<sub>&#x2609;</sub>', layout=widgets.Layout(overflow='visible'))\n",
    "AU_label = widgets.HTML('AU', layout=widgets.Layout(overflow='visible'))\n",
    "deg_label = widgets.HTML('&deg;', layout=widgets.Layout(overflow='visible'))\n",
    "semimajor_cntl = widgets.HBox([semimajor_slider, semimajor_readout, AU_label], \n",
    "                              layout=widgets.Layout(width=ControlColWidth, \n",
    "                                                    overflow='visible'))\n",
    "ecc_cntl = widgets.HBox([ecc_slider, ecc_readout], \n",
    "                        layout=widgets.Layout(width=ControlColWidth, \n",
    "                                              overflow='visible'))\n",
    "incl_cntl = widgets.HBox([incl_slider, incl_readout, deg_label], \n",
    "                         layout=widgets.Layout(width=ControlColWidth, \n",
    "                                               overflow='visible'))\n",
    "phi_cntl = widgets.HBox([phi_slider, phi_readout, deg_label], \n",
    "                        layout=widgets.Layout(width=ControlColWidth, \n",
    "                                              overflow='visible'))\n",
    "\n",
    "##\n",
    "## Orbital playback controls\n",
    "##\n",
    "phase_slider = widgets.IntSlider(value=0,\n",
    "                                   min=0,\n",
    "                                   max=N,\n",
    "                                   step=1,\n",
    "                                   description=\"Phase\",\n",
    "                                   style = {'description_width': 'initial'},\n",
    "                                   disabled=False,\n",
    "                                   continuous_update=True,\n",
    "                                   orientation='horizontal',\n",
    "                                   readout=False,\n",
    "                                   readout_format='.0f',\n",
    "                                   layout=widgets.Layout(width=slider_width, min_width=slider_minwidth,\n",
    "                                                         overflow='visible') )\n",
    "phase_play = widgets.Play(interval = 1, \n",
    "                          value = phase_slider.min, \n",
    "                          min=phase_slider.min, \n",
    "                          max=phase_slider.max, \n",
    "                          step=2, \n",
    "                          description=\"Press play\", \n",
    "                          disabled=False, \n",
    "                          show_repeat=True,\n",
    "                          layout=widgets.Layout(overflow='visible'))\n",
    "widgets.jslink((phase_play, 'value'), (phase_slider, 'value'))\n",
    "phase_cntl = widgets.HBox([phase_slider, phase_play], \n",
    "                          layout=widgets.Layout(width=SimWidth, \n",
    "                                                overflow='visible'))\n",
    "##\n",
    "## Create a viewer controls for scale, zoom, and grid drawing\n",
    "##\n",
    "\n",
    "zoom_slider = widgets.FloatSlider(value=1,\n",
    "                                  min=1,\n",
    "                                  max=300,\n",
    "                                  step=0.1,\n",
    "                                  description=\"View Magnification Factor\",\n",
    "                                  style = {'description_width': 'initial'},\n",
    "                                  disabled=False,\n",
    "                                  continuous_update=True,\n",
    "                                  orientation='horizontal',\n",
    "                                  readout=True,\n",
    "                                  readout_format='.1f',\n",
    "                                  layout=widgets.Layout(width=ControlColWidth, min_width=ControlColWidth, \n",
    "                                                        overflow='visible') )\n",
    "\n",
    "scale_force = widgets.Checkbox(value=False,\n",
    "                               description='Force objects\\' sizes to correct scale with orbit',\n",
    "                               disabled=False,\n",
    "                               layout=widgets.Layout(width=ControlColWidth,\n",
    "                                                     overflow='visible') )\n",
    "draw_grid = widgets.Checkbox(value=True,\n",
    "                               description='Draw grid and orbit',\n",
    "                               disabled=False,\n",
    "                               layout=widgets.Layout(width=ControlColWidth,\n",
    "                                                     overflow='visible') )\n",
    "\n",
    "##\n",
    "## Create text boxes for reporting certain system parameters\n",
    "##\n",
    "P_output = widgets.Text(value = str(0),\n",
    "                        description = 'Orbital Period (Days)',\n",
    "                        style = {'description_width': 'initial'},\n",
    "                        disabled = True, \n",
    "                        layout=widgets.Layout(width=inform_width, height='40px',\n",
    "                                              overflow='visible'))\n",
    "\n",
    "ap_output = widgets.Text(value = str(0),\n",
    "                         description = 'Periastron (R<sub>&#x2609;</sub>)',\n",
    "                         style = {'description_width': 'initial'},\n",
    "                         disabled = True, \n",
    "                         layout=widgets.Layout(width=inform_width, \n",
    "                                               overflow='visible'))\n",
    "\n",
    "aa_output = widgets.Text(value = str(0),\n",
    "                         description = 'Apastron ($R_\\odot$)',\n",
    "                         style = {'description_width': 'initial'},\n",
    "                         disabled = True, \n",
    "                         layout=widgets.Layout(width=inform_width, \n",
    "                                               overflow='visible'))\n",
    "\n",
    "gridsep_output = widgets.Text(value = str(0),\n",
    "                              description = 'Grid Spacing ($R_\\odot$)',\n",
    "                              style = {'description_width': 'initial'},\n",
    "                              disabled = True, \n",
    "                              layout=widgets.Layout(width=inform_width, \n",
    "                                                    overflow='visible'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initialize orbit control sliders\n",
    "semimajor_slider.value = init_a\n",
    "ecc_slider.value = init_ecc\n",
    "incl_slider.value = init_incl\n",
    "phi_slider.value = init_phi\n",
    "\n",
    "# Initialize stellar mass sliders\n",
    "starmass_slider.value = init_mass\n",
    "exomass_slider.value = init_mass\n",
    "\n",
    "# Set initial parameters based on stellar mass assuming main sequence stars\n",
    "(starradius.value, startemp.value, hexcolor1) = star.ConfigStar(starmass_slider.value)\n",
    "(exoradius.value, exotemp.value, hexcolor2) = ConfigPlanet(exomass_slider.value,\n",
    "                                                           semimajor_slider.value*AU2RSun,\n",
    "                                                           startemp.value,\n",
    "                                                           starradius.value)\n",
    "    \n",
    "# Set luminosities\n",
    "L1_output.value = str(nf.SigFig((startemp.value/star.Te_Sun)**4 * starradius.value**2, 3))\n",
    "L2_output.value = str(nf.SigFig((exotemp.value/star.Te_Sun)**4 * exoradius.value**2, 3))\n",
    "\n",
    "# Initialize time index\n",
    "t_idx = 0\n",
    "\n",
    "# Build the initial Binary Star Model with radial velocity and light curves turned off\n",
    "# and semimajor axis in AU\n",
    "bsm = star.BinaryStarModel(mass1=starmass_slider.value, \n",
    "                           mass2=exomass_in_MSun.value,\n",
    "                           rad1=starradius.value,\n",
    "                           rad2=exoradius_in_RSun.value,\n",
    "                           temp1=startemp.value,\n",
    "                           temp2=exotemp.value,\n",
    "                           a=semimajor_slider.value,\n",
    "                           e=ecc_slider.value,\n",
    "                           phi=phi_slider.value,\n",
    "                           N=N, \n",
    "                           Na=Na, \n",
    "                           Ntheta=Ntheta, \n",
    "                           rv_init=False, \n",
    "                           lc_init=False,\n",
    "                           a_in_AU=True)\n",
    "\n",
    "# Convert units\n",
    "P_output.value = \"{0:.1f}\".format(bsm.P)\n",
    "ap_output.value = \"{0:.1f}\".format(bsm.ap*AU2RSun)\n",
    "aa_output.value = \"{0:.1f}\".format(bsm.aa*AU2RSun)\n",
    "\n",
    "##\n",
    "## Set Up 3D Simulation and controls for left side\n",
    "##\n",
    "\n",
    "# creates the object that gets displayed to the screen (defaulting to showing orbital paths)\n",
    "exosystem_view = star.BinaryStarViewer(bsm=bsm, \n",
    "                                       t_idx0=t_idx,\n",
    "                                       lock_scale=scale_force.value,\n",
    "                                       enable_zoom=False,\n",
    "                                       draw_grid=draw_grid.value,\n",
    "                                       draw_orbits=draw_grid.value,\n",
    "                                       view_width=view_width, \n",
    "                                       view_height=view_height)\n",
    "exosystem_renderer = exosystem_view.renderer\n",
    "\n",
    "#\n",
    "# Construct left half controls\n",
    "#\n",
    "\n",
    "# Spacer widget\n",
    "spacer = widgets.HTML('<p>', layout=widgets.Layout(width='10px', overflow='visible'))\n",
    "\n",
    "# Create play button to control phase value automatically\n",
    "phase_title_default = '<b>Controls for Orbital Motion</b>:'\n",
    "phase_title_collision = '<B style=\"color:red\">COLLISION DETECTED! CONTROLS DISABLED!</B>'\n",
    "phase_title = widgets.HTML(phase_title_default)\n",
    "phase_controls = widgets.VBox([phase_title, phase_cntl], \n",
    "                              layout=widgets.Layout(width=SimWidth, \n",
    "                                                   overflow='visible'))\n",
    "\n",
    "# Creates System Parameter Controls (e.g. orbital property controls)\n",
    "sys_title_default = '<b>System Parameters</b>:'\n",
    "sys_title_scaling = '<b>System Parameters</b> (<B style=\"color:red\">Objects\\' sizes are NOT TO SCALE!</B>):'\n",
    "if (exosystem_view.multiplier > 1):\n",
    "    sys_title = widgets.HTML(value=sys_title_scaling, layout=widgets.Layout(overflow='visible'))\n",
    "else:\n",
    "    sys_title = widgets.HTML(value=sys_title_default, layout=widgets.Layout(overflow='visible'))\n",
    "    \n",
    "starorbit_controls = widgets.VBox([sys_title, \n",
    "                                   semimajor_cntl,\n",
    "                                   widgets.HBox([spacer, widgets.VBox([P_output])]),\n",
    "                                   ecc_cntl,  \n",
    "                                   incl_cntl, \n",
    "                                   phi_cntl, \n",
    "                                   spacer],\n",
    "                                  layout=widgets.Layout(overflow='visible'))\n",
    "\n",
    "# Assemble items in left column\n",
    "sim_view = widgets.HBox([widgets.HTML('<h3>Model View</h3>'), exosystem_renderer],  \n",
    "                                         layout=widgets.Layout(width=ControlColWidth,      \n",
    "                                                               flex_flow='column',\n",
    "                                                               align_items='center',\n",
    "                                                               align_contents='center',\n",
    "                                                               overflow='visible'))\n",
    "left_column = widgets.VBox([sim_view, phase_controls, starorbit_controls], \n",
    "                           layout=widgets.Layout(width=ControlColWidth,\n",
    "                                                 overflow='visible'))\n",
    "\n",
    "\n",
    "##\n",
    "## Build Controls on right side\n",
    "##\n",
    "\n",
    "# Compute light curve or radial velocity curve and select which to display initially\n",
    "graph_update()\n",
    "\n",
    "# Creates Stellar Mass Controls\n",
    "star_title = widgets.HTML('<b>Viewer Properties</b>:', layout=widgets.Layout(overflow='visible'))\n",
    "starmass_title = widgets.HTML(value=\"<b>Star Parameters</b>\", layout=widgets.Layout(overflow='visible'))\n",
    "exomass_title = widgets.HTML(value=\"<b>Planet Parameters</b>\", layout=widgets.Layout(overflow='visible'))\n",
    "star_controls = widgets.VBox([star_title, draw_grid, scale_force, zoom_slider, starmass_title, starmass_cntl, star_info, \n",
    "                               exomass_title, exomass_cntl, exo_info],\n",
    "                            layout=widgets.Layout(overflow='visible'))\n",
    "    \n",
    "# Creates a vertical box for the right control panel\n",
    "graph_box = widgets.VBox([fig_selector, graph_fig], \n",
    "                            layout=widgets.Layout(width=ControlColWidth,\n",
    "                                                  overflow='visible') )\n",
    "right_column = widgets.VBox([graph_box, star_controls], \n",
    "                            layout=widgets.Layout(width=ControlColWidth,\n",
    "                                                  overflow='visible') )\n",
    "\n",
    "# Places the figure, sliders, and output into a Vbox. The figure is \n",
    "# alone in the top, while the sliders and output are in a Hbox\n",
    "# inside the bottom of the Vbox.\n",
    "MainDisplay = widgets.HBox([left_column, spacer, right_column])\n",
    "\n",
    "# Sets the dimensions of the box. Sets the entire width and the height of \n",
    "# just the top.\n",
    "MainDisplay.layout.width = EntireWidth\n",
    "MainDisplay.layout.overflow = 'visible'\n",
    "display(MainDisplay)\n",
    "\n",
    "##\n",
    "## Turn on interactivity by linking sliders to functions changing the simulation.\n",
    "##\n",
    "\n",
    "# This control determines which graph to display\n",
    "fig_selector.observe(graph_update, names=['value'])\n",
    "\n",
    "# If mass changes, we also need to check for radius and temperature changes\n",
    "starmass_link = traitlets.directional_link((starmass_slider, 'value'), (bsm, 'mass1'))\n",
    "exomass_link = traitlets.directional_link((exomass_in_MSun, 'value'), (bsm, 'mass2'))\n",
    "\n",
    "# Handle changes in orbit directly in binary star model object\n",
    "a_link = traitlets.directional_link((semimajor_slider, 'value'), (bsm, 'a'))\n",
    "ecc_link = traitlets.directional_link((ecc_slider, 'value'), (bsm, 'e'))\n",
    "phi_link = traitlets.directional_link((phi_slider, 'value'), (bsm, 'phi'))\n",
    "\n",
    "# Call property_update if anything that potentially changes the orbit or luminosity is changed\n",
    "starmass_slider.observe(property_update, names=['value'])\n",
    "exomass_slider.observe(property_update, names=['value'])\n",
    "semimajor_slider.observe(property_update, names=['value'])\n",
    "ecc_slider.observe(property_update, names=['value'])\n",
    "phi_slider.observe(property_update, names=['value'])\n",
    "\n",
    "# Handle changes in inclination directly in binary star model and binary star viewer\n",
    "incl_link = traitlets.directional_link((incl_slider, 'value'), (bsm, 'incl'))\n",
    "incl_viewlink = traitlets.directional_link((incl_slider, 'value'), (exosystem_view, 'incl'))\n",
    "incl_slider.observe(inclination_update, names=['value'])\n",
    "\n",
    "# This slider just changes what phase of the orbit to display\n",
    "phase_link = traitlets.directional_link((phase_slider, 'value'), (exosystem_view, 't_idx'))\n",
    "phase_slider.observe(position_update, names=['value'])\n",
    "\n",
    "# If user wants to switch scaling, link that here\n",
    "orbit_change_link = traitlets.directional_link((scale_force, 'value'), (exosystem_view, 'lock_scale'))\n",
    "orbit_change_link = traitlets.directional_link((draw_grid, 'value'), (exosystem_view, 'draw_grid'))\n",
    "orbit_change_link = traitlets.directional_link((draw_grid, 'value'), (exosystem_view, 'draw_orbits'))\n",
    "scale_force.observe(scene_update, names=['value'])\n",
    "draw_grid.observe(scene_update, names=['value'])\n",
    "zoom_slider.observe(zoom_update, names=['value'])\n",
    "\n",
    "# if Binary Star Model orbit model changes, trigger the Binary Star Viewer to change\n",
    "orbit_change_link = traitlets.directional_link((bsm, 'mdl_counter'), (exosystem_view, 'mdl_counter'))"
   ]
  }
 ],
 "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.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}