Repackaging isochrones from MIST

The purpose of this notebook is to keep track of how we generate the file containing the isochrone data.

The source of the data is: https://waps.cfa.harvard.edu/MIST/

Download the synthetic CMDs for SDSS (both v/vcrit = 0 and 0.4) from the “Packaged Model Grids” tab, and decompress using the following lines:

wget https://waps.cfa.harvard.edu/MIST/data/tarballs_v1.2/MIST_v1.2_vvcrit0.0_SDSSugriz.txz

wget https://waps.cfa.harvard.edu/MIST/data/tarballs_v1.2/MIST_v1.2_vvcrit0.4_SDSSugriz.txz

tar xvf MIST_v1.2_vvcrit0.0_SDSSugriz.txz

tar xvf MIST_v1.2_vvcrit0.4_SDSSugriz.txz

To read these files, a module is provided, see https://waps.cfa.harvard.edu/MIST/resources.html .

Alternate link to this module: jieunchoi/MIST_codes

Note one difficulty: the number of rows (i.e., initial masses) differs for each age, making it not that trivial to store the whole data into a multidimensional array. We therfore opt for the list structure detailed below.

import read_mist_models
import pathlib
import numpy as np
import astropy.table
import pickle

# metalicities fe/H (as they appear in the filenames)
fehs = np.array([-4, -3.5, -3, -2.5, -2, -1.75, -1.5, -1.25, -1, -0.75, -0.5, -0.25, 0.0, 0.25, 0.5])
# rotation velocities (as they appear in the filenames)
vrots = np.array([0, 0.4])

# How the metalicity is written in the filename:
def feh_string(feh):
    if feh < -0.01:
        return f"m{-feh:.2f}"
    else:
        return f"p{feh:.2f}"
    
# How rotation is written in the filename:
def vrot_string(vrot):
    return f"vvcrit{vrot:.1f}"

# Filename:
def isocmd_filename(vrot, feh):
    return str(pathlib.Path(f"./MIST_v1.2_{vrot_string(vrot)}_SDSSugriz") / f"MIST_v1.2_feh_{feh_string(feh)}_afe_p0.0_{vrot_string(vrot)}_SDSSugriz.iso.cmd")

# Testing:
#print(isocmd_filename(rots[1], fehs[0]))

# Obtainign some metadata
isocmd = read_mist_models.ISOCMD(isocmd_filename(0.0, 0.0))

print( 'version: ', isocmd.version)
print( 'photometric system: ', isocmd.photo_sys)
print( 'abundances: ', isocmd.abun)
print( 'rotation: ', isocmd.rot)
print( 'ages: ', [round(x,2) for x in isocmd.ages])
print( 'number of ages: ', isocmd.num_ages)
print( 'available columns: ', isocmd.hdr_list)
print( 'Av extinction: ', isocmd.Av_extinction)

print(f"Log10 of ages in Gyr: {[float(round(x,2)) for x in isocmd.ages]}")
print(f"Available columns: {isocmd.hdr_list}")

#for feh in fehs:
#    isocmd = read_mist_models.ISOCMD(isocmd_filename(0.0))
#    print( 'number of ages: ', isocmd.num_ages)
#    print(f"Log10 of ages in Gyr: {[float(round(x,2)) for x in isocmd.ages]}")

# Ages are the same for all metallicities:
logages = np.array(isocmd.ages)

In principle, the data to draw an isochrone is a 2D table: columns are absolute magnitudes in g, r, i, and rows correspond to different stellar masses. We have such a table for each age, and for each metallicity. These table don’t have the same size, as the number of masses changes with age! We store them in a list of lists, with indices (metallicity, age).

# Assembling all the data:

columns_to_extract = ['SDSS_g', 'SDSS_r', 'SDSS_i', 'initial_mass']
corresponding_new_names = ['g', 'r', 'i', 'initial_mass']
assert len(columns_to_extract) == len(corresponding_new_names)

vrot_list = []
for vrot in vrots:
    feh_list = []
    for feh in fehs:#[:2]:
        age_list = []
        isocmd = read_mist_models.ISOCMD(isocmd_filename(vrot, feh))

        for age_index in range(isocmd.num_ages): # [0, 50, 106]
            column_list = []
            for column_to_extract in columns_to_extract:
                #print(column_to_extract)
                column_data = np.array(isocmd.isocmds[age_index][column_to_extract], dtype=np.float32)
                #print(data)
                column_list.append(column_data)
            isochrone = astropy.table.Table(column_list, names=corresponding_new_names)
            
            age_list.append(isochrone)
        feh_list.append(age_list)
    vrot_list.append(feh_list)

to_write = (vrot_list, vrots, fehs, logages)

with open('isochrones_MIST.pickle', 'wb') as f:
    pickle.dump(to_write, f, pickle.HIGHEST_PROTOCOL)

# how to read:

with open('isochrones_MIST.pickle', 'rb') as f:
    (isochrones, vrots, metallicities, logages) = pickle.load(f)

print(len(isochrones))
print(len(isochrones[0]))
print(len(isochrones[0][0]))
print(len(isochrones[0][0][0]))