Interactives/HR_Diagram.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# HR Diagram Interactives"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This Jupyter Notebook was initially developed by Juan Cabanela over the Summer of 2018 for use in \n",
    "# Minnesota State University Moorhead's introductory astronomy lab activities.\n",
    "\n",
    "from IPython.display import display\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import bqplot as bq\n",
    "import ipywidgets as widgets\n",
    "import tempNcolor as t2c"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "##\n",
    "## Define some variables affecting all the plots here\n",
    "##\n",
    "\n",
    "OverlayColor = ['darkturquoise']\n",
    "PlotBkgStyle={'fill': 'black'}\n",
    "PlotLayout={'width': '400px', 'min_height': '400px'}\n",
    "\n",
    "# Import some data from MSUM Stellar Catalog (which is really a modified version of Hipparcos data catalog)\n",
    "MSUMdata = pd.read_csv('data/Hipparcos.csv')\n",
    "\n",
    "# Compute RGB colors for points based on B-V color\n",
    "hexcolors = t2c.rgb2hex(t2c.bv2rgb(MSUMdata['B-V']))\n",
    "hexlist = hexcolors.tolist()\n",
    "\n",
    "# Import data for four star clusters and subset it by cluster\n",
    "rawdata = pd.read_csv('data/clusterdata.csv')\n",
    "\n",
    "brightness = rawdata['Brightness']\n",
    "ColorIndex = rawdata['B-V']\n",
    "tuc = (rawdata['Cluster'] == '47tuc')\n",
    "pleiades = (rawdata['Cluster'] == 'pleiades')\n",
    "hyades = (rawdata['Cluster'] == 'hyades')\n",
    "m53 = (rawdata['Cluster'] == 'm53')\n",
    "\n",
    "# Compute RGB colors for stars in star clusters based on B-V color\n",
    "hexcolors_tuc = t2c.rgb2hex(t2c.bv2rgb(ColorIndex[tuc]))\n",
    "hexlist_tuc = hexcolors_tuc.tolist()\n",
    "hexcolors_m45 = t2c.rgb2hex(t2c.bv2rgb(ColorIndex[pleiades]))\n",
    "hexlist_m45 = hexcolors_m45.tolist()\n",
    "hexcolors_m53 = t2c.rgb2hex(t2c.bv2rgb(ColorIndex[m53]))\n",
    "hexlist_m53 = hexcolors_m53.tolist()\n",
    "hexcolors_hyades = t2c.rgb2hex(t2c.bv2rgb(ColorIndex[hyades]))\n",
    "hexlist_hyades = hexcolors_hyades.tolist()\n",
    "\n",
    "# Define Cluster information Pandas DataFrame\n",
    "\n",
    "cluster_info = pd.DataFrame(columns=['Name', 'LongName', 'Dist', 'FeH', 'minBright', 'maxBright'])\n",
    "cluster_info['index'] = ['m45', 'hyades', 'm53', '47tuc']\n",
    "cluster_info['Name'] = ['Pleiades', 'Hyades', 'M53', '47 Tuc']\n",
    "cluster_info['LongName'] = ['Pleiades (Solar Composition)', 'Hyades (Solar Composition)', 'M53 (Low in Metals)', '47 Tuc (Low in Metals)']\n",
    "cluster_info['Dist'] = [ 136, 47, 18000, 4000 ]\n",
    "cluster_info['FeH'] = [ 0, 0, -1.86 , -0.7 ]\n",
    "cluster_info['maxBright'] = [ 1e-2, 1e-2, 1e-3 , 1e-3 ]\n",
    "cluster_info['minBright'] = 1e-8*cluster_info['maxBright']\n",
    "cluster_info = cluster_info.set_index('index')\n",
    "\n",
    "# Set the initial cluster in each plot\n",
    "init_cluster = 'm45'\n",
    "\n",
    "# load the data for the interpolated stellar evolution models\n",
    "evol_data = pd.read_csv('data/interpolated_evolution.csv')\n",
    "\n",
    "# Generate a list of colors and add to that Pandas DataFrame\n",
    "evol_data['hexcolor'] = t2c.rgb2hex(t2c.bv2rgb(evol_data['B-V']))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Interactive Figure 1: A Comparison of an HR diagram to a Brightness diagram\n",
    "\n",
    "This interactive figure was built using the MSUM Stellar catalog of about 1000\n",
    "stars with measured B-V colors and known distances.  Recall that to measure the\n",
    "B-V color of a star, we have to have measured the brightness (also called\n",
    "\"flux\") of the star in the B and V filters.  To turn a brightness measurement\n",
    "into a luminosity estimate, we also need an estimate of the distance to the\n",
    "star.\n",
    "\n",
    "- The plot on the left is the Luminosity versus B-V color, that is, the **HR Diagram** for the MSUM Stellar Catalog.  \n",
    "- The plot on the right is those same stars, but plotted up based on their observed brightnesses and B-V color, a **Brightness Diagram**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you select a region on the HR diagram, you will see the same stars will be\n",
    "highlighted on the Brightness diagram.   This lets you see that the same stars\n",
    "can appear on very different regions in these two kinds of plots.  *Why is\n",
    "that?*\n",
    "\n",
    "The checkbox and slider on the far right allow you to highlight stars in a\n",
    "given distance range on the brightness diagram.  Examine the appearance of a\n",
    "group of stars in a narrow *range* of distances (like 100 pc) on the brightness\n",
    "diagram.  Select another group of stars at another distance but similar range.\n",
    "\n",
    "Can you explain why a group of stars close to each other on the HR diagram can\n",
    "appear at such varied positions on the Brightness Diagram?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "##\n",
    "## Define functions for first interactive to highlight a given distance range\n",
    "##\n",
    "def plot_distance_subset(mindist, maxdist):\n",
    "    global Dist_bright\n",
    "    \n",
    "    dist_subset = MSUMdata[(MSUMdata['distance']>=mindist) & (MSUMdata['distance']<=maxdist)]\n",
    "    Dist_bright.x = dist_subset['B-V']\n",
    "    Dist_bright.y = dist_subset['Brightness']\n",
    "    return \n",
    "\n",
    "def plot_distance_visibility(change):\n",
    "    global Dist_bright\n",
    "    \n",
    "    # Turn on distance group indicator\n",
    "    Dist_bright.visible = dist_toggle.value\n",
    "    if (Dist_bright.visible):\n",
    "        scat_bright.default_opacities = [0.9]\n",
    "        # Deselect any highlighted region when asking for distance range\n",
    "        Highlight_lum.selected = []\n",
    "        react_region_selector(change)\n",
    "        region_selector.brushing = True   # resets selector to be brushing (not sure if this is working)\n",
    "    else:\n",
    "        scat_bright.default_opacities = [1]\n",
    "    return \n",
    "\n",
    "# This function will respond to changes in slider by selecting subset of data.\n",
    "def dist_range_selection(change):\n",
    "    global min_dist0, max_dist0, dist_range, force_dist_range\n",
    "    max_range = dist_range\n",
    "    \n",
    "    if (change.new[1] != change.old[1]):\n",
    "        # Top slider changing\n",
    "        maxdist = dist_selector.value[1]\n",
    "        \n",
    "        # Adjust range is that is desired behavior\n",
    "        if force_dist_range:\n",
    "            mindist = maxdist - max_range\n",
    "            if (mindist < max_range):\n",
    "                mindist = 0\n",
    "                maxdist = max_range\n",
    "        else:\n",
    "            mindist = dist_selector.value[0]\n",
    "        \n",
    "    else:\n",
    "        # Bottom slider changing\n",
    "        mindist = dist_selector.value[0]\n",
    "        # Adjust range is that is desired behavior\n",
    "        if force_dist_range:\n",
    "            maxdist = mindist + max_range\n",
    "            if (maxdist > max_dist0):\n",
    "                mindist = max_dist0-max_range\n",
    "                maxdist = max_dist0\n",
    "        else:\n",
    "            maxdist = dist_selector.value[1]\n",
    "                \n",
    "    # move sliders if necessary\n",
    "    if force_dist_range:\n",
    "        dist_selector.value= (mindist, maxdist)\n",
    "    \n",
    "    plot_distance_subset(mindist, maxdist)\n",
    "    return\n",
    "\n",
    "# This function takes the information from the BrushSelector and acts on it\n",
    "def react_region_selector(change):\n",
    "    global transparent\n",
    "    \n",
    "    if Highlight_lum.selected is not None:\n",
    "        # Identify the selected stars\n",
    "        idx_selected = Highlight_lum.selected\n",
    "\n",
    "        # Assume zero opacity to start\n",
    "        opacities_hl = list(transparent)\n",
    "        for idx in idx_selected:\n",
    "            opacities_hl[idx]=1\n",
    "\n",
    "        # Turn on the highlighted points (make them opaque instead of transparent)\n",
    "        Highlight_bright.default_opacities = list(opacities_hl)\n",
    "        Highlight_lum.default_opacities = list(opacities_hl)\n",
    "        \n",
    "        # Disable distance subset if necessary\n",
    "        if (len(Highlight_lum.selected) > 0):\n",
    "            dist_toggle.value = False\n",
    "    else:\n",
    "        # Make all highlights transparent again\n",
    "        opacities_hl = list(transparent)\n",
    "        Highlight_bright.default_opacities = list(opacities_hl)\n",
    "        Highlight_lum.default_opacities = list(opacities_hl)        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "##\n",
    "## Define scales, axes, and tooltips, which could be similar between \n",
    "## various implementations of the displayed widgets\n",
    "##\n",
    "\n",
    "# Scales to transform the data\n",
    "min_loglum = np.floor(np.min(np.log10(MSUMdata['Lum'])))\n",
    "max_loglum = np.ceil(np.max(np.log10(MSUMdata['Lum'])))\n",
    "min_lum = 10**min_loglum\n",
    "max_lum = 10**max_loglum\n",
    "min_logbright = np.floor(np.min(np.log10(MSUMdata['Brightness'])))\n",
    "max_logbright = min_logbright + (max_loglum-min_loglum)  # Set to have same scale\n",
    "min_bright = 10**min_logbright\n",
    "max_bright = 10**max_logbright\n",
    "\n",
    "x_sc_color = bq.LinearScale(min =-1, max=2)\n",
    "y_sc_bright = bq.LogScale(min = min_bright, max=max_bright)\n",
    "y_sc_lum = bq.LogScale(min = min_lum, max=max_lum)\n",
    "\n",
    "# Labels and scales for Axes\n",
    "ax_x_color = bq.Axis(label='B-V', scale=x_sc_color, tick_format='.1f', num_ticks=7)\n",
    "ax_y_bright = bq.Axis(label='Brightness (L_sun/pc^2)', scale=y_sc_bright, orientation='vertical')\n",
    "ax_y_lum = bq.Axis(label='Luminosity (L_sun)', scale=y_sc_lum, orientation='vertical')\n",
    "\n",
    "# moving the label perpendicular to the axis\n",
    "ax_y_bright.label_offset = '3.5em'\n",
    "ax_y_lum.label_offset = '3em'\n",
    "\n",
    "# Define the distance limits for brightness plot\n",
    "min_dist0 = 0\n",
    "max_dist0 = 5000\n",
    "init_dist = 1000\n",
    "dist_range = 100\n",
    "force_dist_range = False   # If true, distance range is forced to dist_range by dynamically changing limits"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "##\n",
    "## Construct the actual plots and widgets for the first interactive\n",
    "##\n",
    "\n",
    "# Define a tooltip\n",
    "def_tt_bright = bq.Tooltip(fields=['x', 'y'], formats=['.2f', '.2e'], labels=['B-V', 'Brightness'])\n",
    "def_tt_lum = bq.Tooltip(fields=['x', 'y'], formats=['.2f', '.3f'], labels=['B-V', 'Luminosity'])\n",
    "\n",
    "# Designing the Brightness diagram\n",
    "opaque = list(1*MSUMdata['B-V']/MSUMdata['B-V'])\n",
    "scat_bright = bq.Scatter(x=MSUMdata['B-V'], y=MSUMdata['Brightness'], \n",
    "                         names=MSUMdata['HIP'], display_names=False, \n",
    "                         scales={'x': x_sc_color, 'y': y_sc_bright},\n",
    "                         marker='circle', colors=hexlist, \n",
    "                         default_size=2, default_opacities=opaque,\n",
    "                         tooltip=def_tt_bright, stroke=None, fill=True)\n",
    "scat_bright.tooltip = None\n",
    "\n",
    "# Designing the HR diagram\n",
    "scat_lum = bq.Scatter(x=MSUMdata['B-V'], y=MSUMdata['Lum'], \n",
    "                      names=MSUMdata['HIP'], display_names=False, \n",
    "                      scales={'x': x_sc_color, 'y': y_sc_lum},\n",
    "                      marker='circle', colors=hexlist, \n",
    "                      default_size=2, default_opacities=opaque,\n",
    "                      tooltip=def_tt_lum, stroke=None, fill=True)\n",
    "scat_lum.tooltip = None\n",
    "\n",
    "# This is a second copy of luminosity plot but with zero size points, we can then change the size of a point to\n",
    "# highlight it.\n",
    "transparent = list(0*MSUMdata['B-V'])\n",
    "Highlight_lum = bq.Scatter(x=MSUMdata['B-V'], y=MSUMdata['Lum'], names=MSUMdata['HIP'], display_names=False,\n",
    "                           default_size=20, default_opacities=transparent,\n",
    "                           scales={'x': x_sc_color, 'y': y_sc_lum}, \n",
    "                           marker='circle', colors=OverlayColor, visible = True\n",
    "                          )\n",
    "Highlight_bright = bq.Scatter(x=MSUMdata['B-V'], y=MSUMdata['Brightness'], names=MSUMdata['HIP'], display_names=False,\n",
    "                              default_size=20, default_opacities=transparent,\n",
    "                              scales={'x': x_sc_color, 'y': y_sc_bright}, \n",
    "                              marker='circle', colors=OverlayColor, visible = True\n",
    "                             )\n",
    "# Not sure why default_opacities value isn't sticking in above Scatter calls, but this fixes it - Juan 03.01.2021\n",
    "Highlight_lum.default_opacities = transparent\n",
    "Highlight_bright.default_opacities = transparent\n",
    "\n",
    "# Set up a Brush Selector to select regions of the HR Diagram and highlight the \n",
    "# same regions of the Brightness Diagram\n",
    "region_selector = bq.interacts.BrushSelector(x_scale=x_sc_color, y_scale=y_sc_bright, \n",
    "                                             marks=[Highlight_lum], color=OverlayColor[0])\n",
    "\n",
    "\n",
    "# Initial distance ranges to overlay on Brightness diagram\n",
    "dist_subset = MSUMdata[(MSUMdata['distance']>=init_dist) & (MSUMdata['distance']<=init_dist+dist_range)]\n",
    "Dist_bright = bq.Scatter(x=dist_subset['B-V'], y=dist_subset['Brightness'], scales={'x': x_sc_color, 'y': y_sc_bright},\n",
    "                         marker='circle', colors=OverlayColor, default_size=40, stroke=None, fill=False)\n",
    "Dist_bright.visible = False\n",
    "\n",
    "# Displaying the data\n",
    "fig_lum = bq.Figure(axes=[ax_x_color, ax_y_lum], marks=[scat_lum, Highlight_lum], \n",
    "                    title='HR Diagram of Sample Stars', \n",
    "                    background_style=PlotBkgStyle, layout=PlotLayout, interaction=region_selector)\n",
    "fig_bright = bq.Figure(axes=[ax_x_color, ax_y_bright], marks=[scat_bright, Dist_bright, Highlight_bright], \n",
    "                       title='Brightness of Sample Stars', \n",
    "                       background_style=PlotBkgStyle, layout=PlotLayout)\n",
    "\n",
    "# Add a distance selection slider and toggle switch\n",
    "dist_toggle = widgets.Checkbox(value=False, description='Highlight Stars in', disabled=False )\n",
    "dist_selector = widgets.IntRangeSlider(value=[init_dist, init_dist+dist_range], min=min_dist0, max=max_dist0, step=dist_range, \n",
    "                                       description='Distance Range', disabled=False,\n",
    "                                       continuous_update=True, orientation='vertical', \n",
    "                                       readout=True, readout_format='04d')\n",
    "dist_selector.layout.height = '300px'\n",
    "dist_selector.layout.width = '125px'\n",
    "dist_toggle.layout.width = '300px'\n",
    "                                  \n",
    "dist_control = widgets.VBox([dist_toggle, dist_selector], \n",
    "                            layout=widgets.Layout(align_content='center', align_items='center', \n",
    "                                          display='flex', \n",
    "                                          flex_flow='column', height='400px', max_height='500px', \n",
    "                                          max_width='250px', min_height='50px', min_width='150px', \n",
    "                                          overflow='hidden', width='150px'))\n",
    "\n",
    "# Display it all\n",
    "main_disp = widgets.HBox([fig_lum, fig_bright, dist_control], \n",
    "                     layout=widgets.Layout(align_content='center', align_items='center', \n",
    "                                           display='flex', \n",
    "                                           flex_flow='row', height='500px', max_height='500px', \n",
    "                                           max_width='1000px', min_height='400px', min_width='900px', \n",
    "                                           overflow='hidden', width='1000px'))\n",
    "\n",
    "display(main_disp)\n",
    "\n",
    "# Respond to changes in the distance range selected\n",
    "region_selector.observe(react_region_selector, names=['brushing'])\n",
    "dist_selector.observe(dist_range_selection, 'value')\n",
    "dist_toggle.observe(plot_distance_visibility, 'value')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Interactive Figure 2:  Distance Estimates via Main Sequence Fitting\n",
    "\n",
    "The interactive figure below shows the previous HR diagram on the left and the\n",
    "observed brightness diagrams of several star clusters on the right.\n",
    "\n",
    "Most stars we observe lie on the Main Sequence of the HR diagram, which means\n",
    "if we construct a brightness diagram for stars in a star cluster, we can\n",
    "compare the brightness of the star cluster's main sequece to the KNOWN\n",
    "luminosity of the main sequence stars to estimate the distance.\n",
    "\n",
    "This interactive allows *main sequence fitting*\" of the distance to a star\n",
    "cluster in two steps:\n",
    "\n",
    "1. The HR diagram has blue-green line representing a main sequence on the HR diagram.\n",
    "Initially this line has not been 'fit' to the main sequence.  Click and\n",
    "drag the blue-green dots until the blue-green line representing the main\n",
    "sequence matches the observed main sequence. You are essentially defining the\n",
    "luminosity of the main sequence for a series of B-V colors.\n",
    "\n",
    "2. Did you notice the corresponding blue-green line on the brightness diagram\n",
    "representing the brightness of your main sequence at 10 parcsecs changed as\n",
    "your made your previous adjustments? Adjusting the distance slider will change\n",
    "its brightness to correspond to the new distance.  Adjust the distance until\n",
    "you match the brightness of the main sequence to the observed main sequence for\n",
    "your selected star cluster.  When you do this, you will have estimate the\n",
    "distance to the star cluster via *main sequence fitting*.\n",
    "\n",
    "**Note**: The Pleiades and 47 Tuc data is fairly 'clean' whereas the Hyades and\n",
    "M53 data contains foreground stars contaminating the data.  This is meant to\n",
    "allow you to see a little of what real data looks like that astronomers have to\n",
    "deal with when attempting this method of distance estimation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "##\n",
    "## Define functions for second interactive to allow editable MS on HR and brightness diagrams\n",
    "##\n",
    "\n",
    "# Define function to generate brightness curve for a given distance\n",
    "def generateMScurve(dist):\n",
    "    global MSscat\n",
    "    return MSscat.y/(4*np.pi*dist*dist)\n",
    "\n",
    "# Define function to process distance change\n",
    "def MSdist_changed(change):\n",
    "    global MScurve\n",
    "    MScurve.y = generateMScurve(MSdist_control.value)\n",
    "\n",
    "def update_MS(change):\n",
    "    global MScurve, MSline, MSscat\n",
    "    # update line on screen\n",
    "    MSline.y = MSscat.y\n",
    "    \n",
    "    # Update displayed MS curve on brightness diagram\n",
    "    MScurve.y = generateMScurve(MSdist_control.value)\n",
    "\n",
    "# Define function to select the star cluster we are looking at\n",
    "def SC_changed(change):\n",
    "    global y_sc_cluster\n",
    "    \n",
    "    # Hide all the star clusters\n",
    "    TucData.visible = False\n",
    "    M45Data.visible = False\n",
    "    HyadesData.visible = False\n",
    "    M53Data.visible = False\n",
    "    \n",
    "    # Define new brightness scale for star cluster plots\n",
    "    idx = (cluster_info.Name == change.new)\n",
    "    y_sc_cluster.min = float(cluster_info[idx].minBright.tolist()[0])\n",
    "    y_sc_cluster.max = float(cluster_info[idx].maxBright.tolist()[0])\n",
    "\n",
    "    # Define what data to show\n",
    "    if (change.new == cluster_info.loc['m45'].Name):\n",
    "        M45Data.visible = True\n",
    "    elif (change.new == cluster_info.loc['hyades'].Name):\n",
    "        HyadesData.visible = True\n",
    "    elif (change.new == cluster_info.loc['47tuc'].Name):\n",
    "        TucData.visible = True\n",
    "    elif (change.new == cluster_info.loc['m53'].Name):\n",
    "        M53Data.visible = True\n",
    "    else:  # Shouldn't happen, but default to the Pleiades\n",
    "        M45Data.visible = True\n",
    "\n",
    "        \n",
    "# Create an editable main sequence line\n",
    "MScolor = np.linspace(0, 1.75, 8)\n",
    "MSlum   = np.ones_like(MScolor)    # Initial flat MS curve\n",
    "MSscat = bq.Scatter(x=MScolor, y=MSlum, marker='circle', default_size=50, \n",
    "                  scales={'x': x_sc_color, 'y': y_sc_lum}, colors=OverlayColor,  \n",
    "                  # Make it possible to move points\n",
    "                  enable_move=True, restrict_y = True, continuous_update=True)\n",
    "MSline = bq.Lines(x=MSscat.x, y=MSscat.y, scales={'x': x_sc_color, 'y': y_sc_lum}, colors=OverlayColor)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "##\n",
    "## Adjust scales/axes/tooltips as necessary from the previous plot\n",
    "##\n",
    "\n",
    "# Adjust brightness scale to the star cluster\n",
    "min_bright_cluster = cluster_info.loc[init_cluster].minBright\n",
    "max_bright_cluster = cluster_info.loc[init_cluster].maxBright\n",
    "\n",
    "# Define new brightness scale for star cluster plots\n",
    "y_sc_cluster = bq.LogScale(min = min_bright_cluster, max=max_bright_cluster)\n",
    "ax_y_cluster  = bq.Axis(label='Brightness (L_sun/pc^2)', scale=y_sc_cluster, orientation='vertical')\n",
    "\n",
    "# moving the label perpendicular to the axis\n",
    "ax_y_cluster.label_offset = '3.5em'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "##\n",
    "## This second pair of plots will preserve as much of the interface as the first pair for consistency\n",
    "##\n",
    "\n",
    "# Design the new HR diagram using previous Scatter object\n",
    "fig_lum_cluster = bq.Figure(axes=[ax_x_color, ax_y_lum], marks=[scat_lum, MSline, MSscat],\n",
    "                            animation_duration = 250, title='HR Diagram of Sample Stars',\n",
    "                            background_style=PlotBkgStyle, layout=PlotLayout)\n",
    "\n",
    "# Design brightness diagram for star clusters\n",
    "TucData = bq.Scatter(x=ColorIndex[tuc], y=brightness[tuc], scales={'x': x_sc_color, 'y': y_sc_cluster}, \n",
    "                     colors=hexlist_tuc, default_size=5, stroke=None, fill=True, labels=[cluster_info.loc['47tuc'].Name])\n",
    "TucData.visible = False\n",
    "M45Data = bq.Scatter(x=ColorIndex[pleiades], y=brightness[pleiades], scales={'x': x_sc_color, 'y': y_sc_cluster}, \n",
    "                     colors=hexlist_m45, default_size=5, stroke=None, fill=True,  labels=[cluster_info.loc['m45'].Name])\n",
    "M45Data.visible = True\n",
    "M53Data = bq.Scatter(x=ColorIndex[m53], y=brightness[m53], scales={'x': x_sc_color, 'y': y_sc_cluster}, \n",
    "                     colors=hexlist_m53, default_size=5, stroke=None, fill=True, labels=[cluster_info.loc['m53'].Name])\n",
    "M53Data.visible = False\n",
    "HyadesData = bq.Scatter(x=ColorIndex[hyades], y=brightness[hyades], scales={'x': x_sc_color, 'y': y_sc_cluster}, \n",
    "                     colors=hexlist_hyades, default_size=5, stroke=None, fill=True, labels=[cluster_info.loc['hyades'].Name])\n",
    "HyadesData.visible = False\n",
    "\n",
    "# Generate a brightness curve for MS\n",
    "initial_dist = 10\n",
    "brightnessMS = generateMScurve(initial_dist)\n",
    "MScurve = bq.Lines(x=MScolor, y=brightnessMS, scales={'x': x_sc_color, 'y': y_sc_cluster}, colors=OverlayColor)\n",
    "\n",
    "# Displaying the data\n",
    "fig_bright_cluster = bq.Figure(axes=[ax_x_color, ax_y_cluster], marks=[TucData, M45Data, M53Data, HyadesData, MScurve],\n",
    "                               legend_location='bottom-left', legend_style={'fill': 'white'},\n",
    "                               title='MS Brightness versus Star Cluster Brightness',\n",
    "                               background_style=PlotBkgStyle, layout=PlotLayout)\n",
    "\n",
    "# Select the star cluster\n",
    "SC_select = widgets.RadioButtons(options=list(cluster_info.Name), \n",
    "                                 value=cluster_info.loc[init_cluster].Name, description='Cluster:', disabled=False,\n",
    "                                 layout=widgets.Layout(align_content='center', align_items='center', \n",
    "                                          display='flex', \n",
    "                                          flex_flow='column', height='150px', max_height='200px', \n",
    "                                          max_width='300px', min_height='100px', min_width='125px', \n",
    "                                          overflow='hidden', width='175px'))\n",
    "\n",
    "# Add a slider to allow adjusting distance to main sequence on brightness diagram\n",
    "MSdist_control = widgets.FloatLogSlider(base=10, min=1, max=4.6, step=0.01, description='Distance (pc)', \n",
    "                                     value=initial_dist, orientation='vertical',\n",
    "                                     readout=True, readout_format='04d',\n",
    "                                     layout=widgets.Layout(align_content='center', align_items='center', \n",
    "                                          display='flex', \n",
    "                                          flex_flow='column', height='250px', max_height='400px', \n",
    "                                          max_width='200px', min_height='150px', min_width='50px', \n",
    "                                          overflow='hidden', width='125px'))\n",
    "\n",
    "# Observe the main sequence to see if something is going on\n",
    "MSscat.observe(update_MS, names=['y'])\n",
    "\n",
    "# Next plot the MS line scaled by distance\n",
    "MSdist_control.observe(MSdist_changed, 'value')\n",
    "\n",
    "# Select Cluster to display from widget\n",
    "SC_select.observe(SC_changed, 'value')\n",
    "\n",
    "# Display it all\n",
    "sec_disp = widgets.HBox([fig_lum_cluster, fig_bright_cluster, widgets.VBox([SC_select,MSdist_control])] , \n",
    "                     layout=widgets.Layout(align_content='center', align_items='center', \n",
    "                                           display='flex', \n",
    "                                           flex_flow='row', height='500px', max_height='500px', \n",
    "                                           max_width='1000px', min_height='400px', min_width='900px', \n",
    "                                           overflow='hidden', width='1000px'))\n",
    "\n",
    "display(sec_disp)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Interactive Figure 3: Estimating Age using Stellar Evolution Models\n",
    " \n",
    "The figure below is allows estimation of the age of star cluster by using stellar evolution models.  \n",
    "\n",
    "- The plot on the left shows you where the stars in a star cluster (which have widely varying mass) would appear on the HR diagram as the star cluster ages.\n",
    "- The plot on the right shows a brightness diagram of a real star cluster with the location of the stars in the HR diagram shown as a blue-green line.\n",
    "\n",
    "If you first set the distance to the cluster to the distance you previously estimated in Interactive Figure 2, then you can allow the cluster to age and use this interactive to determine both the distance and age of the star cluster.\n",
    "\n",
    "**Note**: This interactive starts very early on, before the star cluster has really left the gas cloud it formed in.  We never see star clusters only a few hundred thousand years old."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "##\n",
    "## Define functions to control the stellar evolution model\n",
    "##\n",
    "\n",
    "# Define a function to select data by age\n",
    "def generate_age_curve(age, feh):\n",
    "    global init_Z, data2plot, evol_data\n",
    "    \n",
    "    init_selector = (evol_data['Age'] == age) & (evol_data['FeH'] == feh)\n",
    "    return evol_data[init_selector]\n",
    "    \n",
    "# Define function to process distance change\n",
    "def age_changed(change):\n",
    "    global evol_data, age_lum, Age_label, age_bright, ModelDist_control, feh_value\n",
    "\n",
    "    # Update the label for the age\n",
    "    Age_label.value = age_label_str()\n",
    "        \n",
    "    # Update the curve plotted\n",
    "    data2plot = generate_age_curve(Ages[Age_control.value], feh_value)\n",
    "    \n",
    "    age_lum.x = data2plot['B-V']\n",
    "    age_lum.y = data2plot['Luminosity']\n",
    "    age_lum.colors = data2plot['hexcolor'].tolist()\n",
    "    \n",
    "    # Update the brightness curve\n",
    "    ModelDist_changed(ModelDist_control)\n",
    "    \n",
    "    return\n",
    "\n",
    "# Set age labe variably\n",
    "def age_label_str():\n",
    "    global Age_control\n",
    "    # Update the label for the age\n",
    "    if (Ages[Age_control.value] < 1e6):\n",
    "        return \"{0:.0f} years\".format(Ages[Age_control.value])\n",
    "    elif  (Ages[Age_control.value] < 1e9):\n",
    "        return\"{0:.2f} Million years\".format(Ages[Age_control.value]/1e6)\n",
    "    else:\n",
    "        return \"{0:.2f} Billion years\".format(Ages[Age_control.value]/1e9)\n",
    "    \n",
    "    return \"ERROR: You shouldn't see this.\"\n",
    "        \n",
    "# Define function to select the star cluster we are looking at in third figure\n",
    "def SC2_changed(change):\n",
    "    global y_sc2_cluster, TucData2, M45Data2, HyadesData2, M53Data2\n",
    "    # Hide all the star clusters\n",
    "    TucData2.visible = False\n",
    "    M45Data2.visible = False\n",
    "    HyadesData2.visible = False\n",
    "    M53Data2.visible = False\n",
    "    \n",
    "    # Adjust rightness scale for star cluster plots\n",
    "    idx = (cluster_info.LongName == change.new)\n",
    "    y_sc2_cluster.min = float(cluster_info[idx].minBright.tolist()[0])\n",
    "    y_sc2_cluster.max = float(cluster_info[idx].maxBright.tolist()[0])\n",
    "\n",
    "    if (change.new == cluster_info.loc['m45'].LongName):\n",
    "        M45Data2.visible = True\n",
    "    elif (change.new == cluster_info.loc['hyades'].LongName):\n",
    "        HyadesData2.visible = True\n",
    "    elif (change.new == cluster_info.loc['47tuc'].LongName):\n",
    "        TucData2.visible = True\n",
    "    elif (change.new == cluster_info.loc['m53'].LongName):\n",
    "        M53Data2.visible = True\n",
    "    else:  # Shouldn't happen, but default to the Pleiades\n",
    "        M45Data2.visible = True\n",
    "    \n",
    "    # Reset the metallicity\n",
    "    feh_setting(change.new)\n",
    "    \n",
    "\n",
    "# Define the metallicity dataset for the cluster\n",
    "def feh_setting(cluster):\n",
    "    global feh_value, SC2_select\n",
    "    \n",
    "    # Set the metallicity\n",
    "    idx = (cluster_info.LongName == SC2_select.value)\n",
    "    feh_value = float(cluster_info[idx].FeH.tolist()[0])\n",
    "    \n",
    "    # Change label\n",
    "    FeH_label.text = [\"[Fe/H]: {0:.2f}\".format(feh_value)]\n",
    "    \n",
    "    # Redraw the model data if necessary\n",
    "    age_changed(cluster)\n",
    "        \n",
    "# Define function to generate brightness curve for a given distance\n",
    "def generateBrightnessCurve(dist):\n",
    "    global age_lum\n",
    "    return age_lum.y/(4*np.pi*dist*dist)\n",
    "\n",
    "# Define function to deal with distance changes\n",
    "def ModelDist_changed(change):\n",
    "    global ModelDist_control, age_bright, ModelDist_control\n",
    "    \n",
    "    age_bright.x = age_lum.x\n",
    "    age_bright.y = generateBrightnessCurve(ModelDist_control.value)\n",
    "    return\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "##\n",
    "## This this pair of plots will preserve as much of the interface from the previous plots as possible\n",
    "##\n",
    "\n",
    "# Get the metallicities value to assume based on initial cluster\n",
    "feh_value = cluster_info.loc[init_cluster].FeH\n",
    "\n",
    "# Set initial distance\n",
    "initial_dist = 100\n",
    "\n",
    "# Labels and scales for Axes\n",
    "y_sc_lumVage = bq.LogScale(min = 10**(min_loglum+3.5), max=10**(max_loglum+3.5))\n",
    "ax_y_lumVage = bq.Axis(label='Luminosity (L_sun)', scale=y_sc_lumVage, orientation='vertical')\n",
    "ax_y_lumVage.label_offset = '3em'\n",
    "\n",
    "# Get the ages stored for these interpolated models\n",
    "Ages = np.sort(evol_data['Age'].unique())\n",
    "n_ages = len(Ages)\n",
    "\n",
    "# Select model data settings\n",
    "init_age = Ages[0]\n",
    "data2plot = generate_age_curve(init_age, feh_value)\n",
    "\n",
    "# Adjust brightness scale to the star cluster\n",
    "min_bright_cluster = cluster_info.loc[init_cluster].minBright\n",
    "max_bright_cluster = cluster_info.loc[init_cluster].maxBright\n",
    "y_sc2_cluster = bq.LogScale(min = min_bright_cluster, max=max_bright_cluster)\n",
    "\n",
    "# Design brightness diagram for star clusters\n",
    "TucData2 = bq.Scatter(x=ColorIndex[tuc], y=brightness[tuc], scales={'x': x_sc_color, 'y': y_sc2_cluster}, \n",
    "                     colors=hexlist_tuc, default_size=5, stroke=None, fill=True, labels=[cluster_info.loc['47tuc'].Name])\n",
    "TucData2.visible = False\n",
    "M45Data2 = bq.Scatter(x=ColorIndex[pleiades], y=brightness[pleiades], scales={'x': x_sc_color, 'y': y_sc2_cluster}, \n",
    "                     colors=hexlist_m45, default_size=5, stroke=None, fill=True,  labels=[cluster_info.loc['m45'].Name])\n",
    "M45Data2.visible = True\n",
    "M53Data2 = bq.Scatter(x=ColorIndex[m53], y=brightness[m53], scales={'x': x_sc_color, 'y': y_sc2_cluster}, \n",
    "                     colors=hexlist_m53, default_size=5, stroke=None, fill=True, labels=[cluster_info.loc['m53'].Name])\n",
    "M53Data2.visible = False\n",
    "HyadesData2 = bq.Scatter(x=ColorIndex[hyades], y=brightness[hyades], scales={'x': x_sc_color, 'y': y_sc_cluster}, \n",
    "                     colors=hexlist_hyades, default_size=5, stroke=None, fill=True, labels=[cluster_info.loc['hyades'].Name])\n",
    "HyadesData2.visible = False\n",
    "\n",
    "# Select the star cluster\n",
    "SC2_select = widgets.RadioButtons(options=cluster_info.LongName, \n",
    "                                 value=cluster_info.loc[init_cluster].LongName, description='Cluster:', disabled=False,\n",
    "                                 layout=widgets.Layout(align_content='center', align_items='center',\n",
    "                                                       display='flex', flex_flow='column',\n",
    "                                                       overflow='hidden',\n",
    "                                                       height='125px', max_height='200px', min_height='125px',\n",
    "                                                       width='250px', max_width='300px',min_width='100px'))\n",
    "\n",
    "# Define a tooltip\n",
    "tt_age_lum = bq.Tooltip(fields=['name'], formats=['5.2f'], labels=['Mass'])\n",
    "\n",
    "# Designing the HR diagram for a given age of data\n",
    "age_lum = bq.Scatter(x=data2plot['B-V'], y=data2plot['Luminosity'], names=data2plot['Mass'], display_names=False,\n",
    "                     scales={'x': x_sc_color, 'y': y_sc_lumVage}, tooltip=tt_age_lum,\n",
    "                     marker='circle', colors=data2plot['hexcolor'].tolist(), default_size=10, stroke=None, fill=True)\n",
    "\n",
    "# Metallicity label\n",
    "FeH_label = bq.Label(x=[0.75], y=[1e6], scales={'x': x_sc_color, 'y': y_sc_lumVage},\n",
    "                    text=[\"[Fe/H]: {0:.2f}\".format(cluster_info.loc[init_cluster].FeH)],\n",
    "                    default_size=15, font_weight='bolder',\n",
    "                    colors=['white'], update_on_move=False)\n",
    "\n",
    "# Design the new HR diagram using previous Scatter object\n",
    "aged_cluster = bq.Figure(axes=[ax_x_color, ax_y_lumVage], marks=[age_lum, FeH_label], \n",
    "                         title='MIST Model HR Diagram',\n",
    "                         background_style=PlotBkgStyle, layout=PlotLayout)\n",
    "\n",
    "# Compute brightness curve of the model data\n",
    "brightness2plot = generateBrightnessCurve(initial_dist)\n",
    "age_bright = bq.Lines(x=age_lum.x, y=brightness2plot, scales={'x': x_sc_color, 'y': y_sc2_cluster}, colors=OverlayColor)\n",
    "\n",
    "# Displaying the cluster data and the model data on the brightness curve\n",
    "fig2_bright_cluster = bq.Figure(axes=[ax_x_color, ax_y_cluster], \n",
    "                                marks=[TucData2, M45Data2, M53Data2, HyadesData2, age_bright],\n",
    "                               legend_location='bottom-left', legend_style={'fill': 'white'},\n",
    "                               title='Star Cluster Brightness',\n",
    "                               background_style=PlotBkgStyle, layout=PlotLayout)\n",
    "\n",
    "# Add a slider to allow adjusting distance to model data on brightness diagram\n",
    "ModelDist_control = widgets.FloatLogSlider(base=10, min=1, max=4.6, step=0.01, description='Distance (pc)', \n",
    "                                     value=initial_dist, orientation='vertical',\n",
    "                                     readout=True, readout_format='04d',\n",
    "                                     layout=widgets.Layout(align_content='center', align_items='center',\n",
    "                                                           display='flex', flex_flow='column',\n",
    "                                                           overflow='hidden',\n",
    "                                                           height='250px', max_height='400px', min_height='150px',\n",
    "                                                           width='125px', max_width='150px',min_width='50px'))\n",
    "\n",
    "# Add a slider to allow adjusting age of main sequence on HR diagram\n",
    "Age_control = widgets.IntSlider(value = 0, min=0, max=n_ages-1, description='Age of Model Stars',\n",
    "                                     orientation='vertical', readout=False, readout_format='.2e',\n",
    "                                     continuous_update=True,\n",
    "                                     layout=widgets.Layout(align_content='center', align_items='center',\n",
    "                                                           display='flex', flex_flow='column',\n",
    "                                                           overflow='hidden',\n",
    "                                                           height='250px', max_height='400px', min_height='150px',\n",
    "                                                           width='125px', max_width='150px',min_width='50px'))\n",
    "\n",
    "Age_label = widgets.Label(value=age_label_str(),\n",
    "                         layout=widgets.Layout(align_content='center', align_items='center',\n",
    "                                               display='flex', flex_flow='column',\n",
    "                                               height='25px', min_height='25px', max_height='50px',\n",
    "                                               width='125px', min_width='100px', max_width='200px',\n",
    "                                               overflow='hidden'))\n",
    "\n",
    "# Play animation widget\n",
    "play = widgets.Play(interval = 2, value = 0, min=0, max=n_ages-1, step=2, description=\"Press play\", \n",
    "                    disabled=False, show_repeat=False,\n",
    "                   layout=widgets.Layout(height='25px', min_height='25px', max_height='50px',\n",
    "                                         width='125px', min_width='100px', max_width='200px',   \n",
    "                                         overflow='hidden' ))\n",
    "widgets.jslink((play, 'value'), (Age_control, 'value'))\n",
    "\n",
    "# Direct functions to be called if various widgets change\n",
    "SC2_select.observe(SC2_changed, 'value')\n",
    "Age_control.observe(age_changed, 'value')\n",
    "ModelDist_control.observe(ModelDist_changed, 'value')\n",
    "\n",
    "# Build control panel\n",
    "Fig3Panel = widgets.VBox([SC2_select, widgets.HBox([Age_control, ModelDist_control]), Age_label, play],\n",
    "                        layout=widgets.Layout(overflow='hidden'))\n",
    "\n",
    "# Display it all\n",
    "third_disp = widgets.HBox([aged_cluster, fig2_bright_cluster, Fig3Panel] , \n",
    "                     layout=widgets.Layout(align_content='center', align_items='center', \n",
    "                                           display='flex', \n",
    "                                           flex_flow='row', height='500px', max_height='500px', \n",
    "                                           max_width='1000px', min_height='400px', min_width='900px', \n",
    "                                           overflow='hidden', width='1000px'))\n",
    "\n",
    "display(third_disp)\n",
    "\n"
   ]
  }
 ],
 "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": 4
}