{ "cells": [ { "cell_type": "markdown", "metadata": { "collapsed": false }, "source": [ "###### The cell below loads the visual style of the notebook when run." ] }, { "cell_type": "code", "execution_count": 0, "metadata": { "collapsed": false }, "outputs": [ ], "source": [ "from IPython.core.display import HTML\n", "css_file = '../../styles/styles.css'\n", "HTML(open(css_file, \"r\").read())" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false }, "source": [ "# Fitting an isochrone to your data" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false }, "source": [ "
\n", "
\n", "

Learning Objectives

\n", "
\n", "
\n", "\n", "> * Understand the location of co-eval stars in the HR diagram, and how it changes with age.\n", "> * Learn how to use the ```isochrones``` library to create a model of co-eval stars.\n", "> * Compare stellar models to your own open cluster data to find the best-fitting age and errors.\n" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false }, "source": [ "## Theory\n", "\n", "We work on the assumption that all the stars in a cluster are born at the same time (are *co-eval*), and have the same metallicity. Therefore, the main factor that accounts for the different properties of the individual stars is that they have different masses. \n", "\n", "Having produced a colour-magnitude diagram of our open cluster in previous lab sessions, in this session we compare to stellar models in order to calculate the age, distance and [interstellar extinction](http://slittlefair.staff.shef.ac.uk/teaching/phy241/lectures/L03/index.html) of the open cluster." ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "
\n", "
\n", "

Read in your data

\n", "
\n", "
\n", "\n", "> At the end of the previous session, you should have created a CSV file containing apparent $V$ magnitudes, and $B-V$ colours for all of the stars in your open cluster. \n", "\n", "> Upload that file to somewhere on CoCalc (probably easiest to load into the same folder as this notebook). If you haven't finished the previous session yet, use the example CSV file\n", "> included in this assignment ```example_data.csv```, which has only two columns ($V$ and $B-V$)\n", "\n", "> Write some code that reads in the CSV file into two arrays - one for the $V$ magnitudes and one for the $B-V$ colours." ] }, { "cell_type": "code", "execution_count": 0, "metadata": { "collapsed": false }, "outputs": [ ], "source": [ "# WRITE YOUR CODE HERE" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "
\n", "
\n", "

Plot your colour-magnitude diagram

\n", "
\n", "
\n", "\n", "> Write a function that takes these two arrays as the arguments (inputs) and plots a colour-magnitude diagram ($V$ magnitude on the y-axis, $B-V$ on the x-axis). You will want to invert the y-axis so that brighter stars are at the top. You can use the matplotlib function ```axis.invert_yaxis()``` to do this.\n", "\n", "> For use later on, your function will want to **return the matplotlib axis** that is created inside the function.\n", "\n", "> Use your function to plot the colour-magnitude diagram (CMD)." ] }, { "cell_type": "code", "execution_count": 0, "metadata": { "collapsed": false }, "outputs": [ ], "source": [ "# YOUR CODE HERE" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false }, "source": [ "## The `isochrones` library\n", "\n", "The `isochrones` library is not built into Python, but it is installed on the Python3 (System-Wide) kernel on CoCalc. If you are running on your own laptop, you will have to install it yourself (`pip install isochrones`). The `isochrones` library allows you to compute the location in a colour-magnitude diagram of a co-eval population of stars, based on pre-computed grids of evolutionary stellar models.\n", "\n", "You can use several sets of stellar models with the library, but we will use the so-called \"MIST\" stellar models, detailed in [Choi et al 2016](https://ui.adsabs.harvard.edu/abs/2016ApJ...823..102C/abstract).\n", "\n", "You can import the MIST Isochrones using the following code:" ] }, { "cell_type": "code", "execution_count": 0, "metadata": { "collapsed": false, "scrolled": true }, "outputs": [ ], "source": [ "from isochrones import get_ichrone" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false }, "source": [ "Once we've imported the library, you can create an isochrone object by specifying the filters you are interested in, like so:" ] }, { "cell_type": "code", "execution_count": 0, "metadata": { "collapsed": false }, "outputs": [ ], "source": [ "iso = get_ichrone('mist', bands=['Bessell_B', 'Bessell_V'])" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false }, "source": [ "Here we've stored an *isochrone object* in the variable named `iso`. This object contains a single function `isochrone`, that allows us to calculate the CMD location of stars of a given age. The single argument of this function is $\\log_{10}$ of the age in years, so, for a 1.7 Gyr cluster we would use:" ] }, { "cell_type": "code", "execution_count": 0, "metadata": { "collapsed": false }, "outputs": [ ], "source": [ "model = iso.isochrone(9.235) # note the argument is log10(1.7e9)\n", "print(type(model))" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false }, "source": [ "The `isochrone` function returns a `DataFrame` object from the `pandas` library. We haven't time or space to go into the `pandas` library in this course, except to say that many people love it for dealing with data tables. For now, all we need to know is we can access the theoretical $B$ and $V$ magnitudes like so:" ] }, { "cell_type": "code", "execution_count": 0, "metadata": { "collapsed": false }, "outputs": [ ], "source": [ "model_b = model.Bessell_B_mag\n", "model_v = model.Bessell_V_mag\n", "# calculate B-V for this model\n", "model_bv = model_b - model_v" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false }, "source": [ "
\n", "
\n", "

Plot the isochrone

\n", "
\n", "
\n", "\n", "> Using the isochrones library, and your function from earlier, plot an isochrone on your HR diagram. Use an isochrone of approximately the correct age for your cluster.\n" ] }, { "cell_type": "code", "execution_count": 0, "metadata": { "collapsed": false }, "outputs": [ ], "source": [ "iso.isochrone?" ] }, { "cell_type": "code", "execution_count": 0, "metadata": { "collapsed": false }, "outputs": [ ], "source": [ "# WRITE YOUR CODE HERE" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false }, "source": [ "# Distance, extinction and reddening\n", "\n", "If you've done the last step correctly, you should find the isochrone lies nowhere near your data. What's going on?\n", "\n", "The answer is that your data are in **apparent** magnitudes, **and** the light from the stars has been extincted and reddened by dust lying between us and the cluster. The `isochrones` library returns **absolute**, unreddened magnitudes.\n", "\n", "The conversion between absolute and apparent magnitudes depends upon the distance to the cluster:\n", "\n", "$$m - M = 5 \\log_{10} (d/10),$$\n", "\n", "where $d$ is the distance in parsecs. Notice that this will affect the magnitudes, but not the colour! This means the effect of distance is to move the isochrone vertically on the CMD.\n", "\n", "What about reddening and extinction? Just like extinction by dust in our own atmosphere, interstellar extinction will make the stars fainter, and redder. So, whilst a finite distance moves an star vertically in the CMD, dust along the line of sight will move a star down and to the right.\n", "\n", "We can write that the total extinction in the $V$-band, in magnitudes, is $A_V$. This means that if the true magnitude of a star is $V_0$, the observed magnitude is $V = V_0 + A_V$. However, the dust also makes the star redder, so we can write that the true colour of the star is related to the observed colour is $(B-V) = (B-V)_0 + E(B-V)$. The quantity $E(B-V)$ is known as the colour excess, or reddening.\n", "\n", "Since both $A_V$ and $E(B-V)$ are related to the amount of dust between us and the cluster, it should not surprise you that they are related. It turns out that $A_V = 3.1 E(B-V)$; see [Shultz & Wiemer (1975)](https://ui.adsabs.harvard.edu/abs/1975A%26A....43..133S/abstract), for example." ] }, { "cell_type": "markdown", "metadata": { "collapsed": false }, "source": [ "
\n", "
\n", "

Apply distance and reddening

\n", "
\n", "
\n", "\n", "> Write a new function that accepts four arguments: a matplotlib axis, the logarithm of the age, the distance and the extinction, and plots an isochrone on the provided axis. Correct the V magnitudes and the B-V colours of the isochrone for distance and reddening. I've provided a template to get you started.\n", "\n", "> Use your function in combination with others you've plotted above to plot your CMD with three isochrones on top of it. Plot one isochrone of approximately the correct age, distance and reddening. Then plot two more; one a Gyr older, the other a Gyr younger." ] }, { "cell_type": "code", "execution_count": 0, "metadata": { "collapsed": false }, "outputs": [ ], "source": [ "def plot_isochrone(axis, log_age, distance, a_v):\n", " \"\"\"\n", " Plots a MIST group isochrone on the axis supplied by the user\n", "\n", " Parameters\n", " ----------\n", " axis: matplotlib.axis\n", " the matplotlib axis object on which we want to plot\n", " log_age: float\n", " logarithm of the age in years (base 10)\n", " distance: float\n", " distance to cluster in parsecs\n", " a_v: float\n", " V-band extinction to cluster, in magnitudes\n", " \"\"\"\n", " # YOUR CODE HERE" ] }, { "cell_type": "code", "execution_count": 0, "metadata": { "collapsed": false }, "outputs": [ ], "source": [ "# USE YOUR FUNCTION HERE" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false }, "source": [ "---------\n", "
\n", "
\n", "

Homework #7

\n", "

Find the best fitting isochrone

\n", "
\n", "
\n", "\n", "Now it's time to use the code you wrote above to fit the isochrones to your data and find the best fitting model. Follow the instructions below to find the age, distance and reddening to your cluster, with errors.\n", "\n", "By trial and error, using the functions written above, you can alter the distance, reddening and age until the isochrone sequence lines up with the data.\n", "\n", "However, if you **only fit the position of the main-sequence stars**, you could do this for an isochrone of **any** age! Even if the isochrone is the wrong age there will be a value of the reddening and distance which will make the main sequence of the isochrone line up with the main sequence in the data. How then do we find the age of the cluster?\n", "\n", "The answer is that only an isochrone of the correct age will simultaneously fit the main sequence, the location of the M-S turn-off and the location of giant stars in your cluster. If you've chosen the right age, the location of stars in various evolutionary stages will line up with the isochrone - if you have the right distance and reddening.\n", "\n", "You might want to read the companion notebook `understanding_CMDs.ipynb`, which helps explain where stars at various stages of evolution appear in the CMD.\n", "\n", "To put all this together, you can imagine finding the age, distance and reddening all from one colour magnitude diagram. You would follow a process something like this:\n", "\n", "1. Pick an age. \n", "2. Adjust the distance and reddening until the isochrone lies near the data. \n", "3. Try a slightly younger, or older age and repeat step 2.\n", "4. If the fit has got worse, try changing the age in the other direction.\n", "\n", "When you find the right age, distance and redenning, the isochrone should match up with the data!\n", "\n", "
\n", "
\n", "

Q1: Find the best fitting isochrone (4 points)

\n", "
\n", "
\n", "\n", "> Using the method above, find the distance, reddening and age that best fits your cluster. Make a plot of the best fitting age." ] }, { "cell_type": "code", "execution_count": 0, "metadata": { "collapsed": false }, "outputs": [ ], "source": [ "# YOUR CODE HERE" ] }, { "cell_type": "markdown", "metadata": { "collapsed": false }, "source": [ "## Uncertainties\n", "\n", "How to judge our uncertainties in the three quantities? One way is to try out isochrones of different ages. Suppose we try an isochrone which is 200 Myr older than our best fit above. We'd need to change the distance and extinction again to get it as close as possible to the data, but once we'd done that, it might still be a \"reasonable\" fit; i.e it looks OK by-eye. \n", "\n", "What if we tried one 400 Myr older? Perhaps, even after tweaking the distance and reddening, we still had a poor fit to the data, again judged by-eye. In this case, we'd conclude that a cluster 200 Myr older was consistent with our data, but one 400 Myr older is not. Therefore, an uncertainty of 200 Myr on the age would seem reasonable.\n", "\n", "We could use the spreads in distance and extinction we found to estimate uncertainties in those properties, too.\n", "\n", "
\n", "
\n", "

Q2: Finding our errors (4 points)

\n", "
\n", "
\n", "\n", "> Use the approach above to find error bars for your age, distance and extinction to the cluster. Make a plot that has your best fitting isochrone on it, along with the isochrones that represent the worst fit you think is still OK. Clearly label which is which." ] }, { "cell_type": "code", "execution_count": 0, "metadata": { "collapsed": false }, "outputs": [ ], "source": [ "# YOUR CODE HERE." ] }, { "cell_type": "markdown", "metadata": { "collapsed": false }, "source": [ "
\n", "
\n", "

Q3: Systematic Errors (2 points)

\n", "
\n", "
\n", "\n", "> The technique above will give you a good idea of your statistical errors; i.e how the scatter in the CMD affects your ability to accurately judge the age. However, there may be sources of *systematic* error. A systematic error would pull **all the stars** around in the CMD **in the same direction** and therefore systematically affect the age, distance or extinction you derive.\n", "\n", "> In the markdown cell below, briefly outline what these may be, and what effect they would have. \n", "\n", "> Hint: you will want to think carefully about the steps you carried out to calculate the absolute photometry of the cluster." ] }, { "cell_type": "markdown", "metadata": { "collapsed": false }, "source": [ "**Write your answer here**" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (system-wide)", "language": "python", "metadata": { "cocalc": { "description": "Python 3 programming language", "priority": 100, "url": "https://www.python.org/" } }, "name": "python3", "resource_dir": "/ext/jupyter/kernels/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.6.9" } }, "nbformat": 4, "nbformat_minor": 4 }