Interactives/Radioactivity.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Radioactive Decay Interactives\n",
    "\n",
    "## Interactive Figure 1: Model of Radioactive Decay\n",
    "\n",
    "This first figure takes a population of 900 atoms and models their radioactive decay."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This first interactive is designed to allow you to explore what happens during radioactive decay of some isotope. An *isotope* refers to a particular version of an element.  For example, the non-radioactive Carbon-12 isotope has 6 protons and 6 neutrons in its nucleus but the radioactive Carbon-14 isotope has 6 protons and 8 neutrons in its nucleus.  Different isotopes of the same element behave identically in terms of chemistry and bonding to other atoms, but their nuclear properties can differ.\n",
    "\n",
    "During radioactive decay, a *parent isotope* is said to decay into a *daughter isotope*.  So, for example, the parent isotope of Carbon-14 decays into Nitrogen-14.  There is no way to predict when any one nucleus of a *parent isotope* will decay into a *daughter isotope*.   That said, if you look at a large number of nuclei of a parent isotope, they exhibit a very simple property:\n",
    "\n",
    "> **The same *fraction* of a radioactive parent isotope will decay over the same amount of time.**\n",
    "    \n",
    "The time it takes for one-half of a population of parent isotopes to decay (on average) into daughter isotopes is called the *half-life* of that isotope.  Different parent isotopes can have very different half-lifes.   **NOTE:** This interactive shows a simulation of the decay of only 900 atoms.  The radioactive decay is still modelled as occurring randomly for any one atom, so the simulation will show slightly different results on different runs!!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.display import display\n",
    "import numpy as np\n",
    "import bqplot as bq\n",
    "import ipywidgets as widgets\n",
    "import random as random\n",
    "import pandas as pd\n",
    "import number_formatting as nf\n",
    "from math import ceil, floor, log10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Originally developed June 2018 by Samuel Holen\n",
    "##\n",
    "## Edits by Juan Cabanela October 2018 to allow changes in the GUI.\n",
    "## - Made the display of the precise number/faction of parent/daught atoms instructor\n",
    "##   configurable.\n",
    "## - Fixed a problem with the display of data, forced data to only be displayed for first\n",
    "##   10 half-lifes regardless of actual generated decay times (this allows hard-coding of)\n",
    "##   num_ticks and tick values.\n",
    "\n",
    "## Pre-construct model of radioactive decay of a population\n",
    "## of parent and daughter atoms.\n",
    "## \n",
    "\n",
    "# GUI Configuration Parameters\n",
    "show_counts = False\n",
    "max_half_lifes = 10   # maximum half-lives to graph\n",
    "time_ticks = 6  # number of ticks for the time axis\n",
    "# half-lives to place ticks on horizontal axis\n",
    "half_life_ticks = np.linspace(0, max_half_lifes, time_ticks)\n",
    "\n",
    "# Constants Related to decay of the parent species to the daughter species\n",
    "N_parent = 900          # initial number of parent atoms (should be a perfect square)\n",
    "N_daughter = 0          # initial number of daughter atoms\n",
    "tau = 1                 # placeholder for the half-life of the parent species \n",
    "h = 0.025               # time step (in half lives)\n",
    "mu = np.log(2.) / tau   # constant for decay time distribution \n",
    "Plot_all_times = True   # Plot all times (otherwise, selects only before time slider value)\n",
    "\n",
    "# Initialize tracking of number of atoms\n",
    "Parent_counts = []          # list of number of parent atoms \n",
    "Dauther_counts = []          # list of number of daughter atoms\n",
    "\n",
    "# Generate a uniform random distribution of N_parent numbers from 0 to 1\n",
    "z = np.random.rand(N_parent)\n",
    "\n",
    "# Function to convert uniform distribution of random numbers to\n",
    "# a distribution weighted to model radiactive decay. The times are in number of \n",
    "# half-lives. \n",
    "#\n",
    "# The unsorted data representing the number of half-lifes until the individual\n",
    "# decay of each atom.\n",
    "decay_times = -np.log(1 - z) / mu\n",
    "decay_times_sorted = np.sort( decay_times )\n",
    "\n",
    "# Genereate array of numbers of atoms left\n",
    "# Adjusted so that each count contains 0 and N_parent\n",
    "Parent_counts = np.arange(N_parent,-1, -1, dtype='int') # Number of parent atoms\n",
    "Daughter_counts = np.ones_like(Parent_counts)\n",
    "Daughter_counts = N_parent - Parent_counts   # Number of daughter atoms\n",
    "\n",
    "#\n",
    "# Construct Pandas data frames\n",
    "#\n",
    "\n",
    "# Time column adjusted to include t=0\n",
    "decay_data = pd.DataFrame()\n",
    "decay_data['time'] = np.concatenate((np.zeros(1), decay_times_sorted))\n",
    "decay_data['Parent'] = Parent_counts\n",
    "decay_data['Daughter'] = Daughter_counts\n",
    "\n",
    "# Data array for species\n",
    "species = pd.DataFrame()\n",
    "species['parent_long'] = ['Generic', 'Carbon', 'Thallium','Uranium','Rubidium']\n",
    "species['daughter_long'] = ['Generic', 'Nitrogen','Lead','Thorium','Strontium']\n",
    "species['parent_short'] = ['Parent', 'C-14','Tl-208','U-235','Rb-87']\n",
    "species['daughter_short'] = ['Daughter', 'N-14','Pb-208','Th-231','Sr-87']\n",
    "species['half-lives'] = [tau, 5730, 3.053 * 60, 703.8, 48.8]\n",
    "species['step-size'] = [h, 125, 0.05 * 60, 15, 1]\n",
    "species['timeunits'] = ['half-lives', 'years', 'seconds', 'million years', 'billion years']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "##\n",
    "## Define functions to respond to controls on population plots.\n",
    "##\n",
    "\n",
    "def UpdateSpecies(change=None):\n",
    "    ##\n",
    "    ## Deal with possible changes of species\n",
    "    ##\n",
    "    \n",
    "    # Generate a new uniform random distribution of N_parent numbers from 0 to 1 and\n",
    "    # then convert uniform distribution to one weighted to model radioactive decay.\n",
    "    z = np.random.rand(N_parent)\n",
    "    decay_times = -np.log(1 - z) / mu\n",
    "    decay_times_sorted = np.sort( decay_times )\n",
    "    decay_data['time'] = np.concatenate((np.zeros(1), decay_times_sorted))\n",
    "    \n",
    "    # Reset the time to zero\n",
    "    Time_slide.value = 0\n",
    "    \n",
    "    # Adjust half-life data and limits on plots appropriately \n",
    "    new_index = species.loc[species.parent_long == pick_Species.value].index[0]\n",
    "    \n",
    "    # Set the half-life based the selected species\n",
    "    hf = species['half-lives'][new_index]\n",
    "    \n",
    "    # Set the limit on the time sider to be 10 half-lifes\n",
    "    Time_slide.max = max_half_lifes*hf\n",
    "    \n",
    "    # Set the time step on slider\n",
    "    Time_slide.step = species['step-size'][new_index]\n",
    "\n",
    "    # Updates the time slider label/units\n",
    "    unit_label.value = species['timeunits'][new_index]\n",
    "    \n",
    "    # Updates the time axes on the plot\n",
    "    x_time.max = Time_slide.max\n",
    "    ax_x_time.scale = x_time\n",
    "    ax_x_time.label = unit_label.value\n",
    "    \n",
    "    # Update tick values\n",
    "    ax_x_time.tick_values = list(half_life_ticks*hf)\n",
    "\n",
    "    # Update the legend\n",
    "    parent_label_new = species['parent_short'][new_index]\n",
    "    daughter_label_new = species['daughter_short'][new_index]\n",
    "    line_parent.labels = [parent_label_new]\n",
    "    line_daughter.labels = [daughter_label_new]\n",
    "    \n",
    "    # Update the species in the box that shows how many are present    \n",
    "    parent_label.value = parent_label_new + ' produced'\n",
    "    daughter_label.value = daughter_label_new + ' remaining'\n",
    "\n",
    "    # Now call the function for updating the plot in the event of a time change\n",
    "    UpdateTimes()\n",
    "\n",
    "    \n",
    "def UpdateFraction(change=None):\n",
    "    ##\n",
    "    ## Deal with whether fraction view or number view selected to set scales for plot,\n",
    "    ## y values for plot, and label values\n",
    "    ##\n",
    "    if frac_or_num.value == False:\n",
    "        # Number count mode enabled\n",
    "\n",
    "        # Update axes and scales\n",
    "        fig_counts.axes = [ax_x_time, ax_y_number]\n",
    "        line_parent.scales={'x': x_time, 'y': y_number}\n",
    "        line_daughter.scales={'x': x_time, 'y': y_number}\n",
    "        pts_parent.scales={'x': x_time, 'y': y_number}\n",
    "        pts_daughter.scales={'x': x_time, 'y': y_number}\n",
    "        line_time.scales={'x': x_time, 'y': y_number}\n",
    "\n",
    "        # Update 'time line' limits\n",
    "        line_time.y = [0, N_parent]    \n",
    "    else:\n",
    "        # Fraction mode enabled\n",
    "\n",
    "        # Updated axes and scales\n",
    "        fig_counts.axes = [ax_x_time, ax_y_fraction]\n",
    "        line_parent.scales={'x': x_time, 'y': y_fraction}\n",
    "        line_daughter.scales={'x': x_time, 'y': y_fraction}\n",
    "        pts_parent.scales={'x': x_time, 'y': y_fraction}\n",
    "        pts_daughter.scales={'x': x_time, 'y': y_fraction}\n",
    "        line_time.scales={'x': x_time, 'y': y_fraction}\n",
    "\n",
    "        # Update 'time line' limits\n",
    "        line_time.y = [0, 1]\n",
    "\n",
    "    UpdateTimes()\n",
    "\n",
    "    \n",
    "def UpdateTimes(change=None):\n",
    "    ##\n",
    "    ## Deal with changes in the time slider\n",
    "    ##\n",
    "\n",
    "    # Recall half-life of this species\n",
    "    species_idx = species.loc[species.parent_long == pick_Species.value].index[0]\n",
    "    hf = species['half-lives'][species_idx]\n",
    "\n",
    "    # Update time label with 3 significant figures\n",
    "    Time_label.value = str(nf.SigFig(Time_slide.value, 3))\n",
    "        \n",
    "    # Get the array of times\n",
    "    time_arr = hf * decay_data['time']\n",
    "    \n",
    "    # Set the times to be x values\n",
    "    line_daughter.x = time_arr\n",
    "    line_parent.x = line_daughter.x\n",
    "    pts_daughter.x = line_daughter.x\n",
    "    pts_parent.x = line_daughter.x\n",
    "    line_time.x = [Time_slide.value, Time_slide.value]\n",
    "    \n",
    "    # Changes the color of population to reflect the decays\n",
    "    for i in range(N_parent):\n",
    "        if Time_slide.value >= decay_times[i]*hf:\n",
    "            Colors[i] = 'blue'\n",
    "        else:\n",
    "            Colors[i] = 'red'\n",
    "    population_scat.colors = Colors\n",
    "    \n",
    "    # Identify where we are in the data set\n",
    "    i = 0\n",
    "    while i < N_parent + 1 and hf*decay_data['time'][i] < Time_slide.value:        \n",
    "        i += 1\n",
    "    if i > 0:\n",
    "        i -= 1\n",
    "        \n",
    "    # Update selection of the parent and daughter data to be plotted\n",
    "    if Plot_all_times:\n",
    "        daughter_decay = decay_data['Daughter']\n",
    "        parent_decay = decay_data['Parent']\n",
    "    else:\n",
    "        daughter_decay = decay_data['Daughter'][0:i+1]\n",
    "        parent_decay = decay_data['Parent'][0:i+1] \n",
    "    \n",
    "    line_parent.y = parent_decay\n",
    "    line_daughter.y = daughter_decay\n",
    "        \n",
    "    ##\n",
    "    ## Deal with whether fraction view or number view selected to set scales for plot,\n",
    "    ## y values for plot, and label values\n",
    "    ##\n",
    "    if frac_or_num.value == False:\n",
    "        # Number count mode enabled\n",
    "\n",
    "        # Update number of parent and daughter in labels\n",
    "        parent_present.value = str(decay_data['Parent'][i])\n",
    "        daughter_present.value = str(decay_data['Daughter'][i])\n",
    "    else:\n",
    "        # Fraction mode enabled\n",
    "        \n",
    "        # Update the x and y arrays for the parent and daughter lines\n",
    "        line_parent.y = (1/N_parent)*parent_decay\n",
    "        line_daughter.y = (1/N_parent)*daughter_decay\n",
    "    \n",
    "        # Update number of parent and daughter in labels\n",
    "        parent_present.value = '{:.3f}'.format((1/N_parent)*decay_data['Parent'][i])\n",
    "        daughter_present.value = '{:.3f}'.format((1/N_parent)*decay_data['Daughter'][i])\n",
    "    \n",
    "    # Update parent and daughter lines\n",
    "    pts_parent.y = line_parent.y\n",
    "    pts_daughter.y = line_daughter.y\n",
    "    \n",
    "    # Update the tooltip labels and formats (Bug? Doesn't appear to be working)\n",
    "    #pts_parent.tooltip.labels[1] = parent_label_new\n",
    "    #pts_daughter.tooltip.labels[1] = daughter_label_new\n",
    "    #if frac_or_num.value == False:\n",
    "    #    pts_parent.tooltip.formats[1] = '3.0f'\n",
    "    #    pts_daughter.tooltip.formats[1] = '3.0f'\n",
    "    #else:\n",
    "    #    pts_parent.tooltip.formats[1] = '0.3f'\n",
    "    #    pts_daughter.tooltip.formats[1] = '0.3f'\n",
    "        \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "##\n",
    "## Set up counts versus time plot\n",
    "##\n",
    "\n",
    "# Set up Species to be Generic\n",
    "init_species_ind = 0  \n",
    "\n",
    "# Set up initial time\n",
    "init_time = 0\n",
    "init_time_idx = 0\n",
    "\n",
    "# Set up axes\n",
    "x_time = bq.LinearScale(min = 0, max=max_half_lifes)\n",
    "y_number = bq.LinearScale(min = 0, max=N_parent)\n",
    "y_fraction = bq.LinearScale(min = 0, max=1)\n",
    "\n",
    "# Labels and scales for Axes\n",
    "ax_x_time = bq.Axis(label=species['timeunits'][init_species_ind], scale=x_time, num_ticks = time_ticks,\n",
    "                    tick_values = half_life_ticks )\n",
    "ax_y_number = bq.Axis(label='Number of atoms', scale=y_number, orientation='vertical')\n",
    "ax_y_fraction = bq.Axis(label='Fraction of atoms', scale=y_fraction, orientation='vertical')\n",
    "\n",
    "# Define tooltip (Bug: doesn't allow relabeling Tooltips, also needs to apply to Scatter, not Lines)\n",
    "#def_tt_parent = bq.Tooltip(fields=['x', 'y'], formats=['.2f', '3.0f'], labels=['time', species['parent_short'][init_species_ind]])\n",
    "#def_tt_daughter = bq.Tooltip(fields=['x', 'y'], formats=['.2f', '3.0f'], labels=['time', species['daughter_short'][init_species_ind]])\n",
    "def_tt_parent = bq.Tooltip(fields=['x', 'y'], formats=['.2f', '.3f'], labels=['time', 'amount of parent isotope'])\n",
    "def_tt_daughter = bq.Tooltip(fields=['x', 'y'], formats=['.2f', '.3f'], labels=['time', 'amount of daughter isotope'])\n",
    "\n",
    "# Define the Lines and Scatter plots\n",
    "# NOTE: Scatter only necessary to allow tooltips to function.\n",
    "init_x = decay_data['time']*species['half-lives'][init_species_ind]\n",
    "if Plot_all_times:\n",
    "    init_parent = decay_data['Parent']\n",
    "    init_daughter = decay_data['Daughter']\n",
    "else:\n",
    "    init_parent = decay_data['Parent'][0:init_time_idx]\n",
    "    init_daughter = decay_data['Daughter'][0:init_time_idx]\n",
    "    \n",
    "pts_parent = bq.Scatter(x=init_x, y=init_parent,\n",
    "                      scales={'x': x_time, 'y': y_number}, marker='circle', default_size=2,\n",
    "                      display_legend=False, colors=['red'], labels=[species['parent_short'][init_species_ind]], \n",
    "                      tooltip=def_tt_parent)\n",
    "pts_daughter = bq.Scatter(x=init_x, y=init_daughter,\n",
    "                        scales={'x': x_time, 'y': y_number}, marker='circle', default_size=2,\n",
    "                        display_legend=False, colors=['blue'], labels=[species['daughter_short'][init_species_ind]],\n",
    "                        tooltip=def_tt_daughter)\n",
    "line_parent = bq.Lines(x=init_x, y=init_parent, \n",
    "                       scales={'x': x_time, 'y': y_number}, display_legend=True, colors=['red'], \n",
    "                       labels=[species['parent_short'][init_species_ind]], )\n",
    "line_daughter = bq.Lines(x=init_x, y=init_daughter, \n",
    "                         scales={'x': x_time, 'y': y_number}, display_legend=True, colors=['blue'], \n",
    "                         labels=[species['daughter_short'][init_species_ind]] )\n",
    "\n",
    "# Set up a vertical line on this plot to indicate the current time\n",
    "times_x = [init_time, init_time]\n",
    "times_y = [0, N_parent]\n",
    "line_time = bq.Lines(x=times_x, y=times_y,\n",
    "                     scales={'x': x_time, 'y': y_number}, \n",
    "                     colors=['greenyellow'])\n",
    "\n",
    "# Creates figure for plot\n",
    "fig_counts = bq.Figure(axes=[ax_x_time, ax_y_number], marks=[line_parent, line_daughter, line_time, pts_parent, pts_daughter], \n",
    "                       legend_location='right', legend_style={'fill': 'white'}, \n",
    "                       title='Counts versus Time', background_style={'fill': 'black'}, \n",
    "                       layout={'width': '500px', 'min_height': '400px'},\n",
    "                      animation=1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Slider widget to control the amount of time that has passed\n",
    "# Set for generic half-life situation initially (0 to 10 half-lifes)\n",
    "Time_slide = widgets.FloatSlider(\n",
    "    value=0.,\n",
    "    description='Time',\n",
    "    min=0.,\n",
    "    max=max_half_lifes,\n",
    "    step=h,\n",
    "    disabled=False,\n",
    "    continuous_update=False,\n",
    "    orientation='horizontal',\n",
    "    readout=False,\n",
    "    readout_format='.1f',\n",
    "    layout=widgets.Layout(overflow='visible',\n",
    "                          width='400px',\n",
    "                          max_width='500px',\n",
    "                          min_width='250px')\n",
    ")\n",
    "\n",
    "# Widget to display the number of parent atoms present\n",
    "parent_present = widgets.Text(\n",
    "    value = str(N_parent),\n",
    "    style = {'description_width': 'initial'},\n",
    "    #description = species['parent_short'][0]+' remaining',\n",
    "    disabled = True,\n",
    "    layout=widgets.Layout(overflow='visible',\n",
    "                          width='175px',\n",
    "                          max_width='300px',\n",
    "                          min_width='125px')\n",
    ")\n",
    "\n",
    "# Widget to display the number of daughter atoms present\n",
    "daughter_present = widgets.Text(\n",
    "    value = str(0),\n",
    "    style = {'description_width': 'initial'},\n",
    "    #description = species['daughter_short'][0]+' produced',\n",
    "    disabled = True,\n",
    "    layout=widgets.Layout(overflow='visible',\n",
    "                          width='175px',\n",
    "                          max_width='300px',\n",
    "                          min_width='125px')\n",
    ")\n",
    "\n",
    "# Widgets to label the time slider with units\n",
    "Time_label = widgets.Label(value=str(Time_slide.value))\n",
    "unit_label = widgets.Label(value=str(species['timeunits'][0]))\n",
    "Composite_Time_Label = widgets.HBox([Time_label, unit_label],\n",
    "                                   layout=widgets.Layout(width='200px'))\n",
    "\n",
    "# Labels for the parent/daughter present displays\n",
    "parent_label = widgets.Label(value=species['parent_short'][0]+' remaining',\n",
    "                             layout=widgets.Layout(width='200px', \n",
    "                                                   overflow='visible') )\n",
    "daughter_label = widgets.Label(value=species['daughter_short'][0]+' produced',\n",
    "                             layout=widgets.Layout(width='200px', \n",
    "                                                   overflow='visible') )\n",
    "\n",
    "# Checkbox to choose whether to display the number of each species\n",
    "# or the fraction of each\n",
    "frac_or_num = widgets.Checkbox(value=False, description='Display as fraction of atoms')\n",
    "\n",
    "# Widget to allow one to choose which species to work with\n",
    "pick_Species = widgets.RadioButtons(options=species['parent_long'][:],\n",
    "                                    value='Generic', description='Species:', disabled=False,\n",
    "                                    layout=widgets.Layout(width='250px'))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Scale for population figure\n",
    "x_sc = bq.LinearScale(min=1, max=np.sqrt(N_parent))\n",
    "y_sc = bq.LinearScale(min=1, max=np.sqrt(N_parent))\n",
    "\n",
    "# Axes for population figure\n",
    "ax_x = bq.Axis(scale=x_sc, num_ticks=0)\n",
    "ax_y = bq.Axis(scale=y_sc, orientation='vertical', num_ticks=0)\n",
    "\n",
    "# Creates an array of x values: [1,2,...,30,1,2,...30,.....,1,2,...,30]\n",
    "x_ls = []\n",
    "for i in range(1,int(np.sqrt(N_parent))+1):\n",
    "    x_ls.append(float(i))\n",
    "x_ls = x_ls * int(np.sqrt(N_parent))\n",
    "x_arr = np.array(x_ls)\n",
    "\n",
    "# Creates an array of y values: [1,1,...,1,2,2,...2,......,30,30,...,30]\n",
    "y_ls = []\n",
    "for i in range(1,int(np.sqrt(N_parent))+1):\n",
    "    y_ls += [float(i)] * int(np.sqrt(N_parent))\n",
    "y_arr = np.array(y_ls)    \n",
    "\n",
    "# Creates a color array with the same number of entries as the number of atoms in\n",
    "# the sample\n",
    "Colors = ['red'] * N_parent\n",
    "\n",
    "# Plot the population model\n",
    "population_scat = bq.Scatter(x=x_arr, y=y_arr, scales={'x': x_sc, 'y': y_sc}, colors =['red'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Picking a new species resets everything\n",
    "pick_Species.observe(UpdateSpecies, names=['value'])\n",
    "\n",
    "# Update view from fraction to/from number\n",
    "frac_or_num.observe(UpdateFraction, names=['value'])\n",
    "\n",
    "# Update times\n",
    "Time_slide.observe(UpdateTimes, names=['value'])\n",
    "Time_label.observe(UpdateTimes, names=['value'])\n",
    "\n",
    "# Figure for the population\n",
    "fig_population = bq.Figure(title='Population of Atoms', marks=[population_scat], axes=[ax_x, ax_y], \n",
    "                background_style={'fill' : 'black'},padding_x = 0.025,\n",
    "                min_aspect_ratio=1, max_aspect_ratio=1)\n",
    "\n",
    "# Boxes to organize display\n",
    "parent_box = widgets.HBox([parent_label, parent_present],\n",
    "                          layout=widgets.Layout(overflow='visible'))\n",
    "daughter_box = widgets.HBox([daughter_label, daughter_present],\n",
    "                            layout=widgets.Layout(overflow='visible'))\n",
    "# Set visibility of the exact counts/fractions\n",
    "if (show_counts == False):\n",
    "    parent_box.layout.visibility = 'hidden'\n",
    "    daughter_box.layout.visibility = 'hidden'\n",
    "\n",
    "value_box = widgets.VBox([parent_box, daughter_box])\n",
    "species_box = widgets.HBox([value_box, pick_Species])\n",
    "\n",
    "slide_box = widgets.HBox([Time_slide, Composite_Time_Label])\n",
    "slide_check_box = widgets.VBox([slide_box, frac_or_num])\n",
    "\n",
    "top_box = widgets.HBox([fig_counts, fig_population],\n",
    "                      layout=widgets.Layout(width='900px'))\n",
    "top_box.children[0].layout.width = '450px'\n",
    "top_box.children[1].layout.width = '450px'\n",
    "\n",
    "bottom_box = widgets.HBox([species_box, slide_check_box],\n",
    "                      layout=widgets.Layout(width='900px'))\n",
    "bottom_box.children[0].layout.width = '450px'\n",
    "bottom_box.children[1].layout.width = '450px'\n",
    "\n",
    "# Final display\n",
    "Final = widgets.VBox([top_box, bottom_box])\n",
    "Final.layout.overflow = 'hidden'\n",
    "display(Final)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Interactive Figure 2: Geochron Plot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Assuming a non-radiogenic isotope (that is, an isotope that is not the result of radioactive decay) that also will not decay, its amount should be constant.  This means that for different mineral samples we can measure the ratio of parent isotope versus the non-radiogenic isotope ($P/D_i$) and daughter isotope ($D$) versus the non-radiogenic isotope ($D/D_i$) to build an geochron plot.  For example, using the following isotopes\n",
    "\n",
    "- $D_i$ (non-radiogenic isotope of daughter element)\n",
    "- $D$ (Daughter Isoptope)\n",
    "- $P$ (Parent isotope)\n",
    "\n",
    "an geochron plot could plot $D/D_i$ versus $P/D_i$.  \n",
    "\n",
    "What sets the *geochron method* (also known as the *isochron method*) apart from the just measuring parent and daughter abundances is the use of the non-radiogenic isotope of the daughter element.  This avoids the assumption of no initial daughter isotope before the rock solidified (radioactive decay can occur while rock is molten).\n",
    "\n",
    "Some minerals in the rock incorporate the parent better than daughter which is why the initial amount of parent \n",
    "isotope versus daughter isotope can vary.  We expect daughter versus non-radiogenic isotope ratio to be constant\n",
    "if we pick the non-radiogenic isotope to be the same element as the daughter isotope.\n",
    "\n",
    "With all this said, it is actually often not this simple as many daughter isotopes are themselves radioactive and decay, leading to a chain of reactions, so comparing abundances of parent to daughter isotopes is not simple.\n",
    "\n",
    "*Note:* The idea for the geochron dating interactive came from a Isochron Diagram Java app at *ScienceCourseware.org*.  However that app had some issues in that it didn't divide by a non-radiogenic isotope (or at least didn't mention it).  In fact, they used $D_i$ for the initial amount of daughter isotope instead of the non-radiogenic isotope of the same element as the daughter isotope."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "##\n",
    "## Define the various isotopes we want to consider\n",
    "##\n",
    "\n",
    "isotope_info = pd.DataFrame(columns=['Name', 'PName', 'PAbbrev', 'DName', 'DAbbrev', 'DiName', 'DiAbbrev', 'HalfLife', 'HLUnits'])\n",
    "isotope_info['index'] = ['generic', 'Rb87']\n",
    "isotope_info['Name'] = ['Generic', 'Rb-87->Sr-87']\n",
    "isotope_info['PName'] = ['Parent', 'Rubidium-87']\n",
    "isotope_info['PAbbrev'] = ['P', 'Rb-87']\n",
    "isotope_info['DName'] = ['Daughter', 'Strontium-87']\n",
    "isotope_info['DAbbrev'] = ['D', 'Sr-87']\n",
    "isotope_info['DiName'] = ['Non-Radiogenic Isotope of Daughter Element', 'Strontium-86']\n",
    "isotope_info['DiAbbrev'] = ['D_i', 'Sr-86']\n",
    "isotope_info['HalfLife'] = [ 1, 48.8 ]\n",
    "isotope_info['HLUnits'] = [ 'half-lives', 'Billion years']\n",
    "isotope_info = isotope_info.set_index('index')\n",
    "\n",
    "# Set initial isotope to plot\n",
    "init_isotope = 'generic'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "##\n",
    "## Define the initial amounts of parent and daughter in the sample.\n",
    "##\n",
    "## In principle, I would change this depending on the isotopes we plot.  But I am only plotting\n",
    "## Rb87 --> Sr-87, since that is the most classical use of this Geochron approach.\n",
    "##\n",
    "\n",
    "# Range of P to D_i fractions and initial amounts of D to D_i to consider\n",
    "P2Di_min = 0.05\n",
    "P2Di_max = 0.40\n",
    "D2Di0_min = 0.05\n",
    "D2Di0_max = 0.75\n",
    "\n",
    "# Generate three mineral samples in different thirds of the entire range\n",
    "range_P2Di  = (P2Di_max-P2Di_min)\n",
    "\n",
    "# Create sample amounts\n",
    "n_samples = 4\n",
    "nums = np.array(list(range(1, n_samples+1)))\n",
    "initial_samples = pd.DataFrame(index=nums)\n",
    "initial_D2Di0 = D2Di0_min + (D2Di0_max - D2Di0_min) * np.random.random()\n",
    "initial_samples['P2Di'] = P2Di_min + (range_P2Di/n_samples) * (nums - np.random.random(n_samples))\n",
    "initial_samples['D2Di'] = initial_D2Di0*np.ones_like(nums)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "##\n",
    "## Define functions to call when building interactive plot\n",
    "##\n",
    "\n",
    "def amt_left(sample_in, taus):\n",
    "    # Generate a sample DataFrame after tau half-lifes given an initial DataFrame\n",
    "    sample = sample_in.copy(deep = True)\n",
    "    sample['P2Di'] = sample_in['P2Di']*((1/2)**(taus))\n",
    "    sample['D2Di'] = sample_in['D2Di'] + sample_in['P2Di']*(1 - (1/2)**(taus))\n",
    "    return sample\n",
    "\n",
    "def line_points(sample):\n",
    "    global x_min, x_max, y_min, y_max, initial_D2Di0\n",
    "    \n",
    "    # Determine the end points of a line going through the sample points.\n",
    "    x_range = x_max - x_min\n",
    "    y_range = y_max - y_min\n",
    "    \n",
    "    # Slope (extrapolate from first two points - could be done by a fit to the points)\n",
    "    slope = (sample['D2Di'][2]-sample['D2Di'][1])/(sample['P2Di'][2]-sample['P2Di'][1])\n",
    "    y_final = initial_D2Di0 + slope*x_range\n",
    "    x_points = (x_min, x_max)\n",
    "    y_points = (initial_D2Di0, y_final)\n",
    "    return x_points, y_points, slope\n",
    "\n",
    "def init2current(samples0, samples):\n",
    "    # Compute the lines connecting initital and final points for plotting\n",
    "    n_pts = len(samples0)\n",
    "\n",
    "    xlist = []\n",
    "    ylist = []\n",
    "    for pt in range(1, n_pts+1):\n",
    "        x = np.array([ samples0['P2Di'][pt], samples['P2Di'][pt] ])\n",
    "        y = np.array([ samples0['D2Di'][pt], samples['D2Di'][pt] ])\n",
    "        xlist.append(x)\n",
    "        ylist.append(y)\n",
    "    \n",
    "    return(xlist, ylist)\n",
    "    \n",
    "def HL_changed(change):\n",
    "    global isotope, sample, initial_samples, dots_current, line_current, connectors, slope_label\n",
    "    \n",
    "    # Determine half-life of this isotope\n",
    "    idx = (isotope_info.Name == isotope.value)\n",
    "    HL = float(isotope_info[idx].HalfLife.tolist()[0])\n",
    "    \n",
    "    # How many half-lives have passed?  Use this to get new sample and line info\n",
    "    this_tau = HL_slider.value / HL\n",
    "    sample = amt_left(initial_samples, this_tau)\n",
    "    x_sample, y_sample, slope =  line_points(sample)\n",
    "    \n",
    "    # Update plot\n",
    "    dots_current.x = sample['P2Di']\n",
    "    dots_current.y = sample['D2Di']\n",
    "    line_current.x = x_sample\n",
    "    line_current.y = y_sample\n",
    "    slope_label.value = 'Slope: {0:0.4f}'.format(slope)\n",
    "    xlist, ylist = init2current(initial_samples, sample)\n",
    "    connectors.x = xlist\n",
    "    connectors.y = ylist\n",
    "    \n",
    "    \n",
    "def isotope_changed(change):\n",
    "    global ax_x_P2Di, ax_y_D2Di, HL_slider, HLlabel, UnitsText, Max_half_lives\n",
    "\n",
    "    # Extract the necessary isotope descriptors from the Pandas DataFrame\n",
    "    idx = (isotope_info.Name == change.new)\n",
    "    HL = float(isotope_info[idx].HalfLife.tolist()[0])\n",
    "    HLUnits = isotope_info[idx].HLUnits.tolist()[0]\n",
    "    PAbbrev = isotope_info[idx].PAbbrev.tolist()[0]\n",
    "    DAbbrev = isotope_info[idx].DAbbrev.tolist()[0]\n",
    "    DiAbbrev = isotope_info[idx].DiAbbrev.tolist()[0]\n",
    "\n",
    "    # Get old half-life\n",
    "    idx_old = (isotope_info.Name == change.old)\n",
    "    HL_old = float(isotope_info[idx_old].HalfLife.tolist()[0])\n",
    "\n",
    "    # Determine current age reading from slider and adjust to new units\n",
    "    init_age = HL_slider.value \n",
    "    \n",
    "    # Hard code generic versus others\n",
    "    if (change.new != isotope_info.loc['generic'].Name):\n",
    "        HL_slider.description = \"Time\"\n",
    "    else: \n",
    "        HL_slider.description = \"Half-lives\"    \n",
    "\n",
    "    # Adjust time scales\n",
    "    if (HL_old < HL):\n",
    "        # Adjust maximum limits first before adjusting values (since new HL > old HL)\n",
    "        HL_slider.max = Max_half_lives*HL\n",
    "        HLlabel.max = HL_slider.max \n",
    "        HL_slider.value = HL*(init_age/HL_old)\n",
    "        HLlabel.value = HL_slider.value       \n",
    "    else:\n",
    "        # Adjust maximum limits after adjusting values (since new HL < old HL)\n",
    "        HL_slider.value = HL*(init_age/HL_old)\n",
    "        HLlabel.value = HL_slider.value\n",
    "        HL_slider.max = Max_half_lives*HL\n",
    "        HLlabel.max = HL_slider.max \n",
    "                \n",
    "    # Set the axes and other labels to display\n",
    "    UnitsText.value = HLUnits\n",
    "    ax_x_P2Di.label = '{0} / {1}'.format(PAbbrev, DiAbbrev)\n",
    "    ax_y_D2Di.label = '{0} / {1}'.format(DAbbrev, DiAbbrev)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "##\n",
    "## Set up isochron plot\n",
    "##\n",
    "\n",
    "# Largest possible fraction of decay (only go out to 5 half-lives)\n",
    "Max_half_lives = 5\n",
    "Max_decay_fraction = 1 - (1/2)**(Max_half_lives)\n",
    "\n",
    "# detemine maximum and minimum values of X and Y axes\n",
    "x_step = 0.05\n",
    "x_min = 0\n",
    "x_max = x_step * ceil(initial_samples['P2Di'][n_samples] / x_step)\n",
    "y_step = 0.04\n",
    "y_min = y_step * floor(initial_D2Di0 / y_step)\n",
    "y_max = y_step * ceil((initial_D2Di0 + initial_samples['P2Di'][n_samples] * Max_decay_fraction) / y_step)\n",
    "\n",
    "# Labels and scales for Axes\n",
    "x_P2Di = bq.LinearScale(min = x_min, max = x_max)\n",
    "y_D2Di = bq.LinearScale(min = y_min, max = y_max)\n",
    "ax_x_P2Di = bq.Axis(label='P / D_i', scale=x_P2Di)\n",
    "ax_y_D2Di = bq.Axis(label='D / D_i', scale=y_D2Di, orientation='vertical')\n",
    "\n",
    "# Set up initial conditions\n",
    "taus = 0    # zero half lives past\n",
    "sample = amt_left(initial_samples, taus)\n",
    "\n",
    "##\n",
    "## Define the lines\n",
    "##\n",
    "\n",
    "# Initial amount of daughter line (with dots for initial amounts of parent)\n",
    "x_init, y_init, slope_init =  line_points(initial_samples)\n",
    "line_initial = bq.Lines(x=x_init, y=y_init, scales={'x': x_P2Di, 'y': y_D2Di}, \n",
    "                   line_style='dashed', colors=['red'], labels=['Initial Sample'])\n",
    "dots_initial = bq.Scatter(x=initial_samples['P2Di'], y=initial_samples['D2Di'], scales={'x': x_P2Di, 'y': y_D2Di}, \n",
    "                   colors=['white'], stroke='red', fill= True, labels=['Initial Isochron'])\n",
    "\n",
    "# Current quantities on isochron line\n",
    "x_sample, y_sample, slope =  line_points(sample)\n",
    "line_current = bq.Lines(x=x_sample, y=y_sample, scales={'x': x_P2Di, 'y': y_D2Di}, \n",
    "                   line_style='solid', colors=['red'], labels=['Current Isochron'])\n",
    "dots_current = bq.Scatter(x=sample['P2Di'], y=sample['D2Di'], scales={'x': x_P2Di, 'y': y_D2Di}, \n",
    "                   colors=['red'], stroke='red', fill= True, labels=['Current Isochron'])\n",
    "\n",
    "# Connect Initial and Current quantities on isochron line\n",
    "xlist, ylist = init2current(initial_samples, sample)\n",
    "connectors = bq.Lines(x=xlist, y=ylist, scales={'x': x_P2Di, 'y': y_D2Di}, \n",
    "                   line_style='dotted', colors=['black'])\n",
    "\n",
    "##\n",
    "## Construct plot\n",
    "##\n",
    "isochron = bq.Figure(axes=[ax_x_P2Di, ax_y_D2Di], \n",
    "                     marks=[connectors, line_initial, dots_initial, line_current, dots_current],\n",
    "                     title='Geochron Diagram', \n",
    "                     layout={'width': '700px', 'height': '500px', \n",
    "                             'max_width': '700px', 'max_height': '500px',\n",
    "                             'min_width': '600px', 'min_height': '400px'})\n",
    "\n",
    "##\n",
    "## Construct controls\n",
    "##\n",
    "\n",
    "# Select Generic or Specific Isotopes\n",
    "isotope = widgets.RadioButtons(options=list(isotope_info.Name), \n",
    "                               value=isotope_info.loc[init_isotope].Name, description='Isotope:', \n",
    "                               disabled=False, \n",
    "                               layout=widgets.Layout(height='75px', max_height='100px', min_height='50px', \n",
    "                                                    width='200px', max_width='300px',  min_width='100px'))\n",
    "isotope.observe(isotope_changed, 'value')\n",
    "\n",
    "# Slider and text field controling age\n",
    "idx = (isotope_info.Name == isotope.value)\n",
    "HL = float(isotope_info[idx].HalfLife.tolist()[0])\n",
    "HLUnits = isotope_info[idx].HLUnits.tolist()[0]\n",
    "\n",
    "HL_slider = widgets.FloatSlider(value=0, min=0, max=Max_half_lives*HL, step=0.02,\n",
    "                                description='Half-lives', disabled=False,\n",
    "                                continuous_update=False, orientation='horizontal',\n",
    "                                readout=False, readout_format='.2f',\n",
    "                                layout=widgets.Layout(height='75px', max_height='100px', min_height='50px', \n",
    "                                                    width='200px', max_width='300px',  min_width='100px'))\n",
    "HL_slider.observe(HL_changed, 'value')\n",
    "\n",
    "# Get units value and units for age label, then apply them\n",
    "HLlabel = widgets.BoundedFloatText(value = HL_slider.value, min = HL_slider.min, max = HL_slider.max, \n",
    "                                   step = HL_slider.step,\n",
    "                                       layout={'width': '75px', 'height': '50px', \n",
    "                                               'max_width': '75px', 'max_height': '75px',\n",
    "                                               'min_width': '50px', 'min_height': '50px'})\n",
    "UnitsText = widgets.Label(value=HLUnits)\n",
    "age_label = widgets.HBox([HLlabel, UnitsText])\n",
    "# Link HL slider with this text\n",
    "widgets.jslink((HL_slider, 'value'), (HLlabel, 'value'))\n",
    "\n",
    "# Describe slope\n",
    "slope_label = widgets.Label(value = 'Slope: {0:0.4f}'.format(slope),\n",
    "                                       layout={'align_items':'center','align_content':'center', \n",
    "                                               'justify_content':'center', \n",
    "                                               'width': '100px', 'height': '50px', \n",
    "                                               'max_width': '100px', 'max_height': '75px',\n",
    "                                               'min_width': '50px', 'min_height': '50px'})\n",
    "\n",
    "\n",
    "controls = widgets.VBox( [isotope, HL_slider, age_label, slope_label], \n",
    "                        layout=widgets.Layout(align_content='center', align_items='center', \n",
    "                                              justify_content='center', \n",
    "                                              width='300px', height='500px', \n",
    "                                              max_width='300px', max_height='500px',\n",
    "                                              min_width='100px', min_height='400px',\n",
    "                                              overflow='hidden') )\n",
    "\n",
    "display(widgets.HBox( [isochron, controls] ) )"
   ]
  }
 ],
 "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
}