The cell below loads the visual style of the notebook when run.
In [1]:
from IPython.core.display import HTML
css_file = '../styles/styles.css'
HTML(open(css_file, "r").read())
Out[1]:

Using the isochrones library


The isochrones library is not a standard Python library. If you are running on the managed desktop computers or your own laptop you will have to install it. Start a terminal (linux/OS X) or start an Anaconda command prompt on the managed desktop and type pip install git+https://github.com/timothydmorton/isochrones. This will install the library and its (many) dependencies.

The guide below shows how to use this library to create isochrones to plot on top of your colour-magnitude diagram for the PHY242 observing project. First we'll import the libraries we will need:

In [2]:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

I have downloaded some $B$ and $V$-band data for NGC 7789 from the WEBDA database. I'll read this in:

In [3]:
bv, v = np.loadtxt('ngc7789.txt',unpack=True)

And we can plot the data. There's a lot of data so I've used transparent plot markers to help us see the structure in the colour-magnitude diagram more clearly.

In [4]:
fig, axis = plt.subplots(figsize=(8,6))
axis.plot(bv,v,'r.',alpha=0.1)
axis.invert_yaxis()
axis.set_xlabel('B-V')
axis.set_ylabel('V')
plt.show()

Now let's try and plot an isochrone on top of it. The isochrones library has a nice interface to the Dartmouth group model isochrones (Dotter et al 2008). We import them using the following Python code:

In [5]:
from isochrones.dartmouth import Dartmouth_Isochrone
/usr/local/lib/python2.7/site-packages/IPython/kernel/__init__.py:13: ShimWarning: The `IPython.kernel` package has been deprecated. You should import from ipykernel or jupyter_client instead.
  "You should import from ipykernel or jupyter_client instead.", ShimWarning)
WARNING:root:PyMultiNest not available; only emcee fits will be possible.

Note that the first time you run this command it will download the Dartmouth isochrones over the internet - this can take a long time! Once we have imported the Dartmouth isochrones, we can create one in the filters we want like so:

In [6]:
iso = Dartmouth_Isochrone(bands=['B','V'])

This gives us an isochrone object. The function we are interested is the isochrone member. This has one mandatory argument, which is $log_{10}$ of the age in Myr. For NGC 7789 (1.7 Gyr) this is roughly 9.235.

In [7]:
model = iso.isochrone(9.235)
print(type(model))
<class 'pandas.core.frame.DataFrame'>

The isochrone function returns a DataFrame object from the pandas library. We haven't time or space to go into pandas except to say that many people love it for dealing with complex data tables. All we need to know is we can access the $B$ and $V$ magnitudes like so:

In [8]:
model_b = model.B_mag
model_v = model.V_mag
# calculate B-V for this model
model_bv = model_b - model_v

Alright! Let's plot this isochrone on top of our data and see how it looks!

In [9]:
# same plotting code as earlier
fig, axis = plt.subplots(figsize=(8,6))
axis.plot(bv,v,'r.',alpha=0.1)
axis.invert_yaxis()
axis.set_xlabel('B-V')
axis.set_ylabel('V')

# now I plot the isochrone
axis.plot(model_bv,model_v)
plt.show()

Distance, extinction and reddening

Wow - that looks like a lousy fit. That is because our data are in apparent, reddened magnitudes and the isochrones provide absolute, unreddened magnitudes.

As a first step in making our isochrone units the same as our data, let's apply the distance modulus. NGC 7789 is at roughly 2330 pc...

In [10]:
distance = 2300
distance_modulus = 5*np.log10(distance/10)
print ('dm = {}'.format(distance_modulus))

model_v = model_v + distance_modulus
# note: distance doesn't affect colour!
dm = 11.8086391801

Let's plot again:

In [11]:
# same plotting code as earlier
fig, axis = plt.subplots(figsize=(8,6))
axis.plot(bv,v,'r.',alpha=0.1)
axis.invert_yaxis()
axis.set_xlabel('B-V')
axis.set_ylabel('V')

# now I plot the isochrone
axis.plot(model_bv,model_v)
plt.show()

Still not right! The isochrone is to the left of the data! This is because of interstellar extinction that both absorbs light and reddens colours. The relationship between reddening $E(B-V)$ and extinction $A_V$ is $A_V = 3.2 E(B-V)$. For NGC 7789, $E(B-V)$ is roughly 0.2:

In [12]:
ebv = 0.2
av = 3.2*ebv

model_v = model_v + av
model_bv = model_bv + ebv

# same plotting code as earlier
fig, axis = plt.subplots(figsize=(8,6))
axis.plot(bv,v,'r.',alpha=0.1)
axis.invert_yaxis()
axis.set_xlabel('B-V')
axis.set_ylabel('V')

# now I plot the isochrone
axis.plot(model_bv,model_v)

axis.set_ylim([22,10])
axis.set_xlim([0.0,2.0])
plt.show()

Finding the age

The plot above doesn't look too bad. We could certainly do better by altering the distance and reddening until the main-sequence lines up with the data.

We could do this for an isochrone of any age of course: even if the isochrone is the wrong age there'll 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?

The answer is that only an isochrone of the correct age will simultaneously fit the main sequence and the location of the turn-off from the main sequence. Thus you can find the age, distance and reddening all from one colour-magnitude diagram by trying isochrones of different ages until you find an acceptable fit!

Naturally you'll have to fiddle a lot to get this right, so it makes sense to wrap the code above into a function that you can call which gives a model isochrone in apparent, reddened magnitudes, given an age, distance and reddening. Something like the code below:

In [13]:
def get_isochrone(log_age,distance,ebv):
    # your code here!
    raise NotImplementedException

BTW, as a hint - there is a very quick and easy way to do the above. Rather than tell you what it is, I'll leave it as an exercise in learning to read the documentation for libraries you use....