Isochrone Color-Magnitude Diagram

This Cookbook recipe demonstrates how to create a color-magnitude diagram (CMD) for an isochrone, using the MSG Python interface.

Preparation

Download isochrone-cmd.tar.bz2 and unpack it into your working directory using the tar utility. This archive contains the following items:

MIST.iso — isochrone data file[1]
read_mist_models.py — Python module for reading MIST.iso

Initialization

To start, import modules and set other parameters:

[1]:

# Import standard modules

import os
import numpy as np
import matplotlib.pyplot as plt
import astropy.constants as con

# Import pymsg and read_mist_models

import pymsg
import read_mist_models

# Set constants

R_SUN = con.R_sun.cgs.value
PC = con.pc.cgs.value

# Set plot parameters

plt.rcParams.update({'font.size': 12})

Create PhotGrid Objects

Next, create a pair of pymsg.PhotGrid objects for interpolating fluxes in the Johnson B and V bands:

[2]:

# Create PhotGrid objects using the demo specgrid

MSG_DIR = os.environ['MSG_DIR']
GRID_DIR = os.path.join(MSG_DIR, 'data', 'grids')
PASS_DIR = os.path.join(MSG_DIR, 'data', 'passbands')

photgrid_B = pymsg.PhotGrid(os.path.join(GRID_DIR, 'sg-demo.h5'),
                            os.path.join(PASS_DIR, 'pb-Generic-Johnson.B-Vega.h5'))
photgrid_V = pymsg.PhotGrid(os.path.join(GRID_DIR, 'sg-demo.h5'),
                            os.path.join(PASS_DIR, 'pb-Generic-Johnson.V-Vega.h5'))

Read Isochrone Data

Now read in the isochrone data file and extract the stellar parameters:

[3]:

# Read isochrone data file

iso = read_mist_models.ISO('MIST.iso')

# Extract stellar parameters

Teff = 10**iso.isos[0]['log_Teff']
logg = iso.isos[0]['log_g']

R = 10**iso.isos[0]['log_R']*R_SUN

Reading in: MIST.iso

Evaluate Irradiances

Using these parameters, evaluate the photometric irradiance in each band:

[4]:

# Set the distance to the standard 10 parsecs used to define
# absolute magnitudes

d = 10*PC

# Evaluate irradiances

n = len(Teff)

F_obs_B = np.empty(n)
F_obs_V = np.empty(n)

for i in range(n):

    # Set up photospheric parameters dict

    x = {'Teff': Teff[i],
         'log(g)': logg[i]}

    # Evaluate irradiances. Use try/execpt clause to deal with
    # points that fall outside the grid

    try:
        F_obs_B[i] = (R[i]/d)**2 * photgrid_B.flux(x)
        F_obs_V[i] = (R[i]/d)**2 * photgrid_V.flux(x)
    except (ValueError, LookupError):
        F_obs_B[i] = np.NAN
        F_obs_V[i] = np.NAN

Plot the CMD

As a final step, convert fluxes to absolute magnitudes and plot the CMD:

[5]:

# Evaluate absolute magnitudes

M_B = -2.5*np.log10(F_obs_B)
M_V = -2.5*np.log10(F_obs_V)

# Plot the CMD

plt.figure()

plt.plot(M_B-M_V, M_V)

plt.xlim(0, 2)
plt.ylim(12.5,-2.5)

plt.xlabel('$B-V$')
plt.ylabel('$M_V$');

../../_images/cookbook_isochrone-cmd_isochrone-cmd_9_0.png

Footnotes