From 9c52e9a596ddb67b709af84d089da0524009af0e Mon Sep 17 00:00:00 2001
From: Karl Gordon <kgordon@stsci.edu>
Date: Fri, 18 Oct 2019 13:54:58 -0400
Subject: [PATCH 1/6] Manually rebasing from old branch

---
 .../padova_evolutionary_track_import.py       |  72 ++++
 beast/physicsmodel/stars/evoltracks.py        | 319 ++++++++++++++++++
 beast/physicsmodel/stars/isochrone.py         |  96 ++++++
 beast/plotting/plot_evoltracks.py             |  75 ++++
 beast/plotting/plot_isochrones.py             |  46 +++
 5 files changed, 608 insertions(+)
 create mode 100644 beast/physicsmodel/stars/ET_misc/padova_evolutionary_track_import.py
 create mode 100644 beast/physicsmodel/stars/evoltracks.py
 create mode 100644 beast/plotting/plot_evoltracks.py
 create mode 100644 beast/plotting/plot_isochrones.py

diff --git a/beast/physicsmodel/stars/ET_misc/padova_evolutionary_track_import.py b/beast/physicsmodel/stars/ET_misc/padova_evolutionary_track_import.py
new file mode 100644
index 000000000..fb35a7c51
--- /dev/null
+++ b/beast/physicsmodel/stars/ET_misc/padova_evolutionary_track_import.py
@@ -0,0 +1,72 @@
+import glob
+import os
+import numpy as np
+from astropy.table import Table, vstack, Column
+
+data = Table()
+phil_phase_data = Table()
+directories = []
+
+data.meta['comments'] = [
+    'Interpolated Padova evolutionary tracks',
+    'For more information on data see http://philrosenfield.github.io/padova_tracks/',
+    'Created for use in the BEAST code (https://github.com/BEAST-Fitting)',
+    'Bolometric magnitudes have been calculated using: 4.77 - 2.5 * log L',
+    'Surface gravities have been calculated using:  -10.616 + log mass + 4.0 * log Teff - log L',
+    "0 : Beginning of Pre Main Sequence", "1 : Beginning of Main Sequence (MS)*",
+    "2 : log Teff minimum on the Main Sequence*", "3 : Main Sequence turn off (MSTO)*",
+    "4 : RG tip: Tip of the RGB (TRGB)*",
+    "5 : Beginning Core He Fusion: Start central He-burning phase*", "6 : End of Core He Fusion*",
+    "7 : Horizonzal Branch", "8 : Beginning of TP-AGB (if the star has the a phase)", "",
+    "* EEP definitions followed Dotter et al. 2016,ApJ,222,8D as close as possible"]
+
+# equivalent evolutionary points
+eeps = np.array(([0] * 200) + [1] * 200 + [2] * 200 + [3] * 500 + [4] * 30
+                + [5] * 500 + [6] * 100 + [8] * 200)
+eeps_hb = np.array([7] * 600 + [8] * 200)
+
+# gets data from directories in current directory
+for dir in os.listdir('.'):
+    directories.append(dir)
+
+# combines data
+for dir in directories:
+    print(dir)
+    z = dir.split('Y')[0][1:]
+    files = glob.glob(dir + '/*')
+    for file in files:
+        new_line = Table.read(file, format='ascii')
+        z_col = Column(np.array([float(z)] * len(new_line)).tolist(), name="Z",
+                       description='Metallicity in solar metallicity').T
+        index = Column(np.arange(0, len(z_col)), name="index")
+        if 'HB' in file:
+            phase_col = eeps_hb[0:len(new_line)]
+        else:
+            phase_col = eeps[0:len(new_line)]
+        # print(phase_col)
+        phil_phase_array = Column(phase_col, name="evolutionary_point",
+                                  description='rosenfield eep').T
+        new_line.add_column(z_col)
+        new_line.add_column(phil_phase_array)
+        new_line.add_column(index)
+        data = vstack([data, new_line])
+
+# specifies header information
+data.rename_column('logAge', 'logA')
+data.rename_column('Mass', 'M_ini')
+data.rename_column('logTe', 'logT')
+data.rename_column('Mbol', 'mbolmag')
+
+data['logA'] = Column(data['logA'], unit='yr', description='Age')
+data['M_ini'] = Column(data['M_ini'], unit='Msun', description='Initial mass')
+data['logT'] = Column(data['logT'], unit='K',
+                      description='Effective temperature')
+data['mbolmag'] = Column(data['mbolmag'],
+                         description='Bolomertric luminosities (magnitudes),\
+                                      calculated using: 4.77 - 2.5 * log L')
+data['logg'] = Column(data['logg'], unit='cm/s**2',
+                      description='surface gravity, calculated using:\
+                                  -10.616 + log mass + 4.0 * log Teff - log L')
+data['C/O'] = Column(data['C/O'], description='Carbon to oxygen ratio')
+
+data.write('Padova_ev_tracks.fits', overwrite=True)
diff --git a/beast/physicsmodel/stars/evoltracks.py b/beast/physicsmodel/stars/evoltracks.py
new file mode 100644
index 000000000..9fa86c71b
--- /dev/null
+++ b/beast/physicsmodel/stars/evoltracks.py
@@ -0,0 +1,319 @@
+# use evolutionary tracks instead of isochrones as the basis of
+# the stellar physicsgrid
+
+import numpy as np
+from scipy.interpolate import interp1d
+
+from astropy.table import Table
+
+
+__all__ = ['EvolTracks', 'ETParsec', 'ETMist']
+
+
+class EvolTracks(object):
+    """
+    Stores the evolutionary tracks interpolated to input mass and
+    metallicity grids with the age grid set by the evolutionary track
+    calculations (may wish to revisit and allow for condensation).
+
+    Attributes
+    ----------
+    logmass : log10 of the masses (solar masses)
+
+    Z : numpy array of floats
+        metallicites of tracks
+
+    data: table
+        columns include:
+        M_ini - initial mass [Msun]
+        M_act - actual mass [Msun]
+        Z - metallicity [??]
+        logL - log luminosity [Lsun]
+        logg - log surface gravity [cm/s^2]
+        logA - log age [years]
+        logT - log surface effective temperature [K]
+        stage - evolutionary stage [index]
+    """
+    def __init__(self):
+        self.name = '<auto>'
+
+        # define axis labels for plotting
+        self.alabels = {'logT': 'log(Teff)',
+                        'logg': 'log(g)',
+                        'logL': 'log(L)',
+                        'logA': 'log(age)',
+                        'phase': 'evol phase',
+                        'M_act': 'log(current mass)',
+                        'M_ini': 'log(initial mass)'}
+
+    def load_orig_tables(self, source):
+        """
+        Read the tracks from  the original files (stub)
+        """
+        print('not implemented')
+
+    def plot_tracks(self, ax, xval='logT', yval='logL',
+                    linestyle='-'):
+        """
+        Plot the tracks with the input x, y choices
+
+        Parameters
+        ----------
+        ax : matplotlib axis
+            where to plot
+
+        xval : str, optional
+            what data for x
+
+        xval : str, optional
+            what data for y
+
+        linestyle : string
+            matplotlib linestyle
+        """
+        if xval not in self.data.keys():
+            raise ValueError("xval choice not in data table")
+        if yval not in self.data.keys():
+            raise ValueError("yval choice not in data table")
+
+        # get uniq M_ini values
+        uvals, indices = np.unique(self.data['M_ini'], return_inverse=True)
+        for k, cval in enumerate(uvals):
+            cindxs = np.where(k == indices)
+            ax.plot(self.data[xval][cindxs], self.data[yval][cindxs],
+                    linestyle=linestyle)
+
+        ax.set_xlabel(self.alabels[xval])
+        ax.set_ylabel(self.alabels[yval])
+
+        if xval == 'logT':
+            xmin, xmax = ax.get_xlim()
+            if xmin < xmax:
+                ax.set_xlim(xmax, xmin)
+
+    def grid_metrics(self,
+                     check_keys=['logL', 'logT', 'logg']):
+        """
+        Compute metrics of the grid
+        Primarily to determine how well parameter space is covered
+
+        Parameters
+        ----------
+        check_keys : string array
+            keys in grid to generage metrics for
+
+        Returns
+        -------
+        metrics : dictonary
+            each entry has an array with [min, max, median, mean] deltas
+            for that grid paramters
+        """
+        # loop over eep values accumulating deltas
+        dvals = {}
+        for cname in check_keys:
+            dvals[cname] = []
+
+        for gparam in ['M_ini']:
+            uvals, indices = np.unique(self.data[gparam], return_inverse=True)
+            for k, cval in enumerate(uvals):
+                cindxs, = np.where(k == indices)
+                for cname in check_keys:
+                    dvals[cname] = np.concatenate(
+                        (dvals[cname],
+                         np.absolute(np.diff(self.data[cname][cindxs]))))
+
+        # compute the metrics
+        metrics = {}
+        for cname in check_keys:
+            metrics[cname] = [np.min(dvals[cname]),
+                              np.max(dvals[cname]),
+                              np.median(dvals[cname]),
+                              np.mean(dvals[cname])]
+
+        return metrics
+
+    def regrid(self,
+               logmass_range=[-1., 2.],
+               logmass_delta=0.05,
+               logL_delta=0.05,
+               logT_delta=0.05):
+        """
+        Interpolate a set of evolutionary tracks to a uniform grid
+        in log(initial mass) and variable grid in stellar age.
+        Use Equivalent Evolutionary Points (EEPs) values to do the
+        mass interpolation.  EEPs are provide as part of the evolutionary
+        tracks.
+
+        Parameters
+        ----------
+        logmass_range : (float, float)
+            range of new mass grid
+            default is -1 to 2 (0.1 to 100 M_sun)
+
+        logmass_delta : float
+            log(mass) delta for new mass grid
+            default is 0.05
+        """
+        # get the unique mass values
+        uvals, indices = np.unique(self.data['M_ini'], return_inverse=True)
+        # print(10**uvals)
+
+        # ensure there are evolutionary tracks spanning
+        #   the min/max of the new grid --> no extrapolation
+        new_min_mass = max([min(uvals), logmass_range[0]])
+        new_max_mass = min([max(uvals), logmass_range[1]])
+
+        n_new_masses = int((new_max_mass - new_min_mass)/logmass_delta) + 1
+        new_mass_vals = np.linspace(new_min_mass, new_max_mass,
+                                    num=n_new_masses)
+        # print(n_new_masses, len(uvals))
+        # print(10**new_mass_vals)
+
+        # setup the new grid
+        new_grid = {}
+        for cname in self.data.keys():
+            new_grid[cname] = np.array([])
+
+        # loop over eep values and interopolate to new mass grid
+        # along constant eep tracks
+        uvals, indices = np.unique(self.data['eep'], return_inverse=True)
+        for k, cval in enumerate(uvals):
+            cindxs, = np.where(k == indices)
+            cur_masses = self.data['M_ini'][cindxs]
+
+            # only interpolate for masses defined for the current eep
+            new_gindxs, = np.where(
+                np.logical_and(min(cur_masses) <= new_mass_vals,
+                               new_mass_vals <= max(cur_masses)))
+
+            for cname in self.data.keys():
+                if cname == 'eep':
+                    vals = np.full((len(new_gindxs)), cval)
+                elif cname == 'M_ini':
+                    vals = new_mass_vals[new_gindxs]
+                else:
+                    f = interp1d(cur_masses, self.data[cname][cindxs])
+                    vals = f(new_mass_vals[new_gindxs])
+
+                new_grid[cname] = np.concatenate((new_grid[cname],
+                                                  vals))
+
+        # update the grid
+        self.data = new_grid
+
+        # setup the new grid
+        new_grid = {}
+        for cname in self.data.keys():
+            new_grid[cname] = np.array([])
+
+        # loop over each mass track and condense
+        uvals, indices = np.unique(self.data['M_ini'], return_inverse=True)
+        for k, cval in enumerate(uvals):
+            cindxs, = np.where(k == indices)
+            delta_logL = np.absolute(np.diff(self.data['logL'][cindxs]))
+            delta_logT = np.absolute(np.diff(self.data['logT'][cindxs]))
+            nindxs = [0]
+            cdelt_logL = 0.0
+            cdelt_logT = 0.0
+            for i in range(len(delta_logL)):
+                cdelt_logL += delta_logL[i]
+                cdelt_logT += delta_logT[i]
+                if ((cdelt_logL > logL_delta) or (cdelt_logT > logT_delta)):
+                    nindxs.append(i)
+                    cdelt_logL = delta_logL[i]
+                    cdelt_logT = delta_logT[i]
+
+            if not max(nindxs) == len(delta_logL) - 1:
+                nindxs.append(len(delta_logL) - 1)
+
+            for cname in self.data.keys():
+                new_grid[cname] = np.concatenate(
+                    (new_grid[cname],
+                     self.data[cname][cindxs][nindxs]))
+
+        # update the grid
+        self.data = new_grid
+
+
+class ETParsec(EvolTracks):
+    """
+    EvolTracks specific to PARSEC calculations
+    """
+
+    source = 'PARSEC'
+
+    def load_orig_tables(self, files):
+        """
+        Read the tracks from  the original files
+
+        files : str
+            file or files with the evolutionary track calculations
+            often each file is for a single initial mass
+        """
+        if not type(files) is list:
+            files = [files]
+
+        mass_act = np.array([])
+        logL = np.array([])
+        logT = np.array([])
+        for cfile in files:
+            a = Table.read(cfile, format='ascii')
+            mass_act = np.concatenate((mass_act, a['MASS'].data))
+            logL = np.concatenate((logL, a['LOG_L'].data))
+            logT = np.concatenate((logT, a['LOG_TE'].data))
+
+        self.data = {}
+        self.data['M_act'] = np.array(mass_act)
+        self.data['logL'] = np.array(logL)
+        self.data['logT'] = np.array(logT)
+
+
+class ETMist(EvolTracks):
+    """
+    MIST evolutionary tracks
+    """
+    source = 'MIST'
+
+    def load_orig_tables(self, files):
+        """
+        Read the tracks from  the original files
+
+        files : str
+            file or files with the evolutionary track calculations
+            often each file is for a single initial mass
+        """
+        if type(files) is not list:
+            files = [files]
+
+        mass_act = np.array([])
+        mass_ini = np.array([])
+        logA = np.array([])
+        logL = np.array([])
+        logT = np.array([])
+        logg = np.array([])
+        phase = np.array([])
+        eep = np.array([])
+        for cfile in files:
+            a = Table.read(cfile, format='ascii', header_start=11)
+            tmass = a['star_mass'].data
+            mass_act = np.concatenate((mass_act, tmass))
+            mass_ini = np.concatenate((mass_ini,
+                                       np.full((len(tmass)), max(tmass))))
+            logA = np.concatenate((logA, np.log10(a['star_age'].data)))
+            logL = np.concatenate((logL, a['log_L'].data))
+            logT = np.concatenate((logT, a['log_Teff'].data))
+            logg = np.concatenate((logg, a['log_g'].data))
+            phase = np.concatenate((phase, a['phase'].data))
+            eep = np.concatenate((eep, range(len(a))))
+
+        self.data = {}
+        self.data['M_act'] = np.log10(mass_act)
+        self.data['M_ini'] = np.log10(mass_ini)
+        self.data['logA'] = logA
+        self.data['logL'] = logL
+        self.data['logT'] = logT
+        self.data['logg'] = logg
+        self.data['phase'] = phase
+        self.data['eep'] = eep
+
+        self.alabels['eep'] = 'EEP'
diff --git a/beast/physicsmodel/stars/isochrone.py b/beast/physicsmodel/stars/isochrone.py
index 05caca5de..6b6b74e8c 100644
--- a/beast/physicsmodel/stars/isochrone.py
+++ b/beast/physicsmodel/stars/isochrone.py
@@ -26,6 +26,17 @@ class Isochrone(object):
     def __init__(self, name="", *args, **kwargs):
         self.name = name
 
+        # define axis labels for plotting
+        self.alabels = {
+            "logT": "log(Teff)",
+            "logg": "log(g)",
+            "logL": "log(L)",
+            "logA": "log(age)",
+            "stage": "evol phase",
+            "M_act": "log(current mass)",
+            "M_ini": "log(initial mass)",
+        }
+
     def metalToFeH(self, metal):
         """
         Convert Z to [Fe/H] values
@@ -126,6 +137,91 @@ def _get_continuous_isochrone(self, *args, **kwargs):
 
         return table
 
+    def plot(self, ax, xval="logT", yval="logL", linestyle="-"):
+        """
+        Plot the isochrones with the input x, y choices
+
+        Parameters
+        ----------
+        ax : matplotlib axis
+            where to plot
+
+        xval : str, optional
+            what data for x
+
+        xval : str, optional
+            what data for y
+
+        linestyle : string
+            matplotlib linestyle
+        """
+        if xval not in self.data.keys():
+            raise ValueError("xval choice not in data table")
+        if yval not in self.data.keys():
+            raise ValueError("yval choice not in data table")
+
+        # get uniq logA values
+        uvals, indices = np.unique(self.data["logA"], return_inverse=True)
+        for k, cval in enumerate(uvals):
+            cindxs = np.where(k == indices)
+            ax.plot(
+                self.data[xval][cindxs], self.data[yval][cindxs], linestyle=linestyle
+            )
+
+        if xval in ["M_ini", "M_act"]:
+            ax.set_xscale("log")
+        if yval in ["M_ini", "M_act"]:
+            ax.set_yscale("log")
+
+        ax.set_xlabel(self.alabels[xval])
+        ax.set_ylabel(self.alabels[yval])
+
+        if xval == "logT":
+            xmin, xmax = ax.get_xlim()
+            ax.set_xlim(xmax, xmin)
+
+    def grid_metrics(self, check_keys=["logL", "logT", "logg"]):
+        """
+        Compute metrics of the grid
+        Primarily to determine how well parameter space is covered
+
+        Parameters
+        ----------
+        check_keys : string array
+            keys in grid to generage metrics for
+
+        Returns
+        -------
+        metrics : dictonary
+            each entry has an array with [min, max, median, mean] deltas
+            for that grid paramters
+        """
+        # loop over eep values accumulating deltas
+        dvals = {}
+        for cname in check_keys:
+            dvals[cname] = []
+
+        for gparam in ["logA"]:
+            uvals, indices = np.unique(self.data[gparam], return_inverse=True)
+            for k, cval in enumerate(uvals):
+                cindxs, = np.where(k == indices)
+                for cname in check_keys:
+                    dvals[cname] = np.concatenate(
+                        (dvals[cname], np.absolute(np.diff(self.data[cname][cindxs])))
+                    )
+
+        # compute the metrics
+        metrics = {}
+        for cname in check_keys:
+            metrics[cname] = [
+                np.min(dvals[cname]),
+                np.max(dvals[cname]),
+                np.median(dvals[cname]),
+                np.mean(dvals[cname]),
+            ]
+
+        return metrics
+
 
 class padova2010(Isochrone):
     def __init__(self):
diff --git a/beast/plotting/plot_evoltracks.py b/beast/plotting/plot_evoltracks.py
new file mode 100644
index 000000000..c167a3b10
--- /dev/null
+++ b/beast/plotting/plot_evoltracks.py
@@ -0,0 +1,75 @@
+"""
+Make a plot of the evolutionary tracks
+"""
+from __future__ import (absolute_import, division, print_function,
+                        unicode_literals)
+
+import matplotlib.pyplot as plt
+
+from beast.physicsmodel.stars.evoltracks import (ETParsec, ETMist)
+from beast.plotting.beastplotlib import initialize_parser
+
+
+if __name__ == '__main__':
+    parser = initialize_parser()
+    parser.add_argument('filename', type=str, nargs='*',
+                        help="name of file(s) with evolutionary track")
+    parser.add_argument('--type', default='mist', choices=['mist', 'parsec'],
+                        help="source of tracks")
+    args = parser.parse_args()
+
+    fig, ax = plt.subplots(2, 2, figsize=(10, 10))
+
+    # switch between the types of tracks
+    if args.type == 'parsec':
+        et = ETParsec()
+    else:
+        et = ETMist()
+
+    # read in the tracks
+    et.load_orig_tables(args.filename)
+
+    orig_metrics = et.grid_metrics()
+
+    # plot, plot, plot
+    nls = ':'
+    et.plot_tracks(ax[0, 0], xval='logT', yval='logL', linestyle=nls)
+    et.plot_tracks(ax[1, 1], xval='logA', yval='M_ini', linestyle=nls)
+    et.plot_tracks(ax[0, 1], xval='phase', yval='M_ini', linestyle=nls)
+    et.plot_tracks(ax[1, 0], xval='logT', yval='logg', linestyle=nls)
+
+    # regrid the evolutionary tracks to uniform log(mass) and variable age
+    print('size orig = ', len(et.data['M_ini']))
+
+    et.regrid(logmass_range=[-1., 3.],
+              logmass_delta=0.05,
+              logT_delta=0.05,
+              logL_delta=0.05)
+
+    print('size regrid = ', len(et.data['M_ini']))
+
+    regrid_metrics = et.grid_metrics()
+    for ckey in regrid_metrics.keys():
+        print(ckey)
+        print(orig_metrics[ckey])
+        print(regrid_metrics[ckey])
+
+    # get the new grid metrics
+    # et.grid_metrics()
+
+    # plot, plot, plot
+    nls = '-'
+    et.plot_tracks(ax[0, 0], xval='logT', yval='logL', linestyle=nls)
+    et.plot_tracks(ax[1, 1], xval='logA', yval='M_ini', linestyle=nls)
+    et.plot_tracks(ax[0, 1], xval='phase', yval='M_ini', linestyle=nls)
+    et.plot_tracks(ax[1, 0], xval='logT', yval='logg', linestyle=nls)
+
+    fig.tight_layout()
+
+    save_name = 'evoltracks'
+    if args.tex:
+        plt.rc({'usetex': True})
+    if args.savefig:
+        fig.savefig('{}.{}'.format(save_name, args.savefig))
+    else:
+        plt.show()
diff --git a/beast/plotting/plot_isochrones.py b/beast/plotting/plot_isochrones.py
new file mode 100644
index 000000000..512aab1c0
--- /dev/null
+++ b/beast/plotting/plot_isochrones.py
@@ -0,0 +1,46 @@
+"""
+Make a plot of the evolutionary tracks
+"""
+import matplotlib.pyplot as plt
+
+from beast.physicsmodel.model_grid import make_iso_table
+
+# from beast.physicsmodel.stars.evoltracks import (ETParsec, ETMist)
+from beast.plotting.beastplotlib import initialize_parser
+
+
+if __name__ == '__main__':
+    parser = initialize_parser()
+    args = parser.parse_args()
+
+    fig, ax = plt.subplots(2, 2, figsize=(10, 10))
+
+    # download the file live from the website
+    savename = '/tmp/padova_iso.csv'
+    iso_fname, oiso = make_iso_table('test', iso_fname=savename,
+                                     logtmin=6.0, logtmax=10.13, dlogt=0.15,
+                                     z=[0.019])
+
+    # plot, plot, plot
+    nls = ':'
+    oiso.plot(ax[0, 0], xval='logT', yval='logL', linestyle=nls)
+    oiso.plot(ax[1, 1], xval='logA', yval='M_ini', linestyle=nls)
+    oiso.plot(ax[0, 1], xval='stage', yval='M_ini', linestyle=nls)
+    oiso.plot(ax[1, 0], xval='logT', yval='logg', linestyle=nls)
+
+    print('size grid = ', len(oiso.data['M_ini']))
+
+    grid_metrics = oiso.grid_metrics()
+    for ckey in grid_metrics.keys():
+        print(ckey)
+        print(grid_metrics[ckey])
+
+    fig.tight_layout()
+
+    save_name = 'evoltracks'
+    if args.tex:
+        plt.rc({'usetex': True})
+    if args.savefig:
+        fig.savefig('{}.{}'.format(save_name, args.savefig))
+    else:
+        plt.show()

From de9415f2bfd9e22e21f844e1ad928ca7f23c6c2c Mon Sep 17 00:00:00 2001
From: Karl Gordon <kgordon@stsci.edu>
Date: Fri, 15 Nov 2019 14:31:10 -0500
Subject: [PATCH 2/6] minor updates

---
 beast/physicsmodel/stars/evoltracks.py | 231 +++++++++++++------------
 1 file changed, 125 insertions(+), 106 deletions(-)

diff --git a/beast/physicsmodel/stars/evoltracks.py b/beast/physicsmodel/stars/evoltracks.py
index 9fa86c71b..840587781 100644
--- a/beast/physicsmodel/stars/evoltracks.py
+++ b/beast/physicsmodel/stars/evoltracks.py
@@ -6,8 +6,10 @@
 
 from astropy.table import Table
 
+from beast.config import __ROOT__
 
-__all__ = ['EvolTracks', 'ETParsec', 'ETMist']
+
+__all__ = ["EvolTracks", "ETParsec", "ETMist"]
 
 
 class EvolTracks(object):
@@ -34,26 +36,28 @@ class EvolTracks(object):
         logT - log surface effective temperature [K]
         stage - evolutionary stage [index]
     """
+
     def __init__(self):
-        self.name = '<auto>'
+        self.name = "<auto>"
 
         # define axis labels for plotting
-        self.alabels = {'logT': 'log(Teff)',
-                        'logg': 'log(g)',
-                        'logL': 'log(L)',
-                        'logA': 'log(age)',
-                        'phase': 'evol phase',
-                        'M_act': 'log(current mass)',
-                        'M_ini': 'log(initial mass)'}
+        self.alabels = {
+            "logT": "log(Teff)",
+            "logg": "log(g)",
+            "logL": "log(L)",
+            "logA": "log(age)",
+            "phase": "evol phase",
+            "M_act": "log(current mass)",
+            "M_ini": "log(initial mass)",
+        }
 
     def load_orig_tables(self, source):
         """
         Read the tracks from  the original files (stub)
         """
-        print('not implemented')
+        print("not implemented")
 
-    def plot_tracks(self, ax, xval='logT', yval='logL',
-                    linestyle='-'):
+    def plot_tracks(self, ax, xval="logT", yval="logL", linestyle="-"):
         """
         Plot the tracks with the input x, y choices
 
@@ -77,22 +81,22 @@ def plot_tracks(self, ax, xval='logT', yval='logL',
             raise ValueError("yval choice not in data table")
 
         # get uniq M_ini values
-        uvals, indices = np.unique(self.data['M_ini'], return_inverse=True)
+        uvals, indices = np.unique(self.data["M_ini"], return_inverse=True)
         for k, cval in enumerate(uvals):
             cindxs = np.where(k == indices)
-            ax.plot(self.data[xval][cindxs], self.data[yval][cindxs],
-                    linestyle=linestyle)
+            ax.plot(
+                self.data[xval][cindxs], self.data[yval][cindxs], linestyle=linestyle
+            )
 
         ax.set_xlabel(self.alabels[xval])
         ax.set_ylabel(self.alabels[yval])
 
-        if xval == 'logT':
+        if xval == "logT":
             xmin, xmax = ax.get_xlim()
             if xmin < xmax:
                 ax.set_xlim(xmax, xmin)
 
-    def grid_metrics(self,
-                     check_keys=['logL', 'logT', 'logg']):
+    def grid_metrics(self, check_keys=["logL", "logT", "logg"]):
         """
         Compute metrics of the grid
         Primarily to determine how well parameter space is covered
@@ -113,30 +117,35 @@ def grid_metrics(self,
         for cname in check_keys:
             dvals[cname] = []
 
-        for gparam in ['M_ini']:
+        for gparam in ["M_ini"]:
             uvals, indices = np.unique(self.data[gparam], return_inverse=True)
             for k, cval in enumerate(uvals):
-                cindxs, = np.where(k == indices)
+                (cindxs,) = np.where(k == indices)
                 for cname in check_keys:
                     dvals[cname] = np.concatenate(
-                        (dvals[cname],
-                         np.absolute(np.diff(self.data[cname][cindxs]))))
+                        (dvals[cname], np.absolute(np.diff(self.data[cname][cindxs])))
+                    )
 
         # compute the metrics
         metrics = {}
         for cname in check_keys:
-            metrics[cname] = [np.min(dvals[cname]),
-                              np.max(dvals[cname]),
-                              np.median(dvals[cname]),
-                              np.mean(dvals[cname])]
+            metrics[cname] = [
+                np.min(dvals[cname]),
+                np.max(dvals[cname]),
+                np.median(dvals[cname]),
+                np.mean(dvals[cname]),
+            ]
 
         return metrics
 
-    def regrid(self,
-               logmass_range=[-1., 2.],
-               logmass_delta=0.05,
-               logL_delta=0.05,
-               logT_delta=0.05):
+    def regrid(
+        self,
+        logmass_range=[-1.0, 2.0],
+        logmass_delta=0.05,
+        condense=False,
+        logL_delta=0.05,
+        logT_delta=0.05,
+    ):
         """
         Interpolate a set of evolutionary tracks to a uniform grid
         in log(initial mass) and variable grid in stellar age.
@@ -150,12 +159,16 @@ def regrid(self,
             range of new mass grid
             default is -1 to 2 (0.1 to 100 M_sun)
 
-        logmass_delta : float
-            log(mass) delta for new mass grid
-            default is 0.05
+        condense : boolean
+            set to condense grid using delta values
+            default = False
+
+        logmass_delta, logL_delta, logT_delta : float
+            deltas for the condensed grid
+            default is 0.05 for all 3
         """
         # get the unique mass values
-        uvals, indices = np.unique(self.data['M_ini'], return_inverse=True)
+        uvals, indices = np.unique(self.data["M_ini"], return_inverse=True)
         # print(10**uvals)
 
         # ensure there are evolutionary tracks spanning
@@ -163,9 +176,8 @@ def regrid(self,
         new_min_mass = max([min(uvals), logmass_range[0]])
         new_max_mass = min([max(uvals), logmass_range[1]])
 
-        n_new_masses = int((new_max_mass - new_min_mass)/logmass_delta) + 1
-        new_mass_vals = np.linspace(new_min_mass, new_max_mass,
-                                    num=n_new_masses)
+        n_new_masses = int((new_max_mass - new_min_mass) / logmass_delta) + 1
+        new_mass_vals = np.linspace(new_min_mass, new_max_mass, num=n_new_masses)
         # print(n_new_masses, len(uvals))
         # print(10**new_mass_vals)
 
@@ -176,63 +188,65 @@ def regrid(self,
 
         # loop over eep values and interopolate to new mass grid
         # along constant eep tracks
-        uvals, indices = np.unique(self.data['eep'], return_inverse=True)
+        uvals, indices = np.unique(self.data["eep"], return_inverse=True)
         for k, cval in enumerate(uvals):
-            cindxs, = np.where(k == indices)
-            cur_masses = self.data['M_ini'][cindxs]
+            (cindxs,) = np.where(k == indices)
+            cur_masses = self.data["M_ini"][cindxs]
 
             # only interpolate for masses defined for the current eep
-            new_gindxs, = np.where(
-                np.logical_and(min(cur_masses) <= new_mass_vals,
-                               new_mass_vals <= max(cur_masses)))
+            (new_gindxs,) = np.where(
+                np.logical_and(
+                    min(cur_masses) <= new_mass_vals, new_mass_vals <= max(cur_masses)
+                )
+            )
 
             for cname in self.data.keys():
-                if cname == 'eep':
+                if cname == "eep":
                     vals = np.full((len(new_gindxs)), cval)
-                elif cname == 'M_ini':
+                elif cname == "M_ini":
                     vals = new_mass_vals[new_gindxs]
                 else:
                     f = interp1d(cur_masses, self.data[cname][cindxs])
                     vals = f(new_mass_vals[new_gindxs])
 
-                new_grid[cname] = np.concatenate((new_grid[cname],
-                                                  vals))
+                new_grid[cname] = np.concatenate((new_grid[cname], vals))
 
         # update the grid
         self.data = new_grid
 
-        # setup the new grid
-        new_grid = {}
-        for cname in self.data.keys():
-            new_grid[cname] = np.array([])
-
-        # loop over each mass track and condense
-        uvals, indices = np.unique(self.data['M_ini'], return_inverse=True)
-        for k, cval in enumerate(uvals):
-            cindxs, = np.where(k == indices)
-            delta_logL = np.absolute(np.diff(self.data['logL'][cindxs]))
-            delta_logT = np.absolute(np.diff(self.data['logT'][cindxs]))
-            nindxs = [0]
-            cdelt_logL = 0.0
-            cdelt_logT = 0.0
-            for i in range(len(delta_logL)):
-                cdelt_logL += delta_logL[i]
-                cdelt_logT += delta_logT[i]
-                if ((cdelt_logL > logL_delta) or (cdelt_logT > logT_delta)):
-                    nindxs.append(i)
-                    cdelt_logL = delta_logL[i]
-                    cdelt_logT = delta_logT[i]
-
-            if not max(nindxs) == len(delta_logL) - 1:
-                nindxs.append(len(delta_logL) - 1)
-
+        # setup the condensed grid
+        if condense:
+            new_grid = {}
             for cname in self.data.keys():
-                new_grid[cname] = np.concatenate(
-                    (new_grid[cname],
-                     self.data[cname][cindxs][nindxs]))
+                new_grid[cname] = np.array([])
 
-        # update the grid
-        self.data = new_grid
+            # loop over each mass track and condense
+            uvals, indices = np.unique(self.data["M_ini"], return_inverse=True)
+            for k, cval in enumerate(uvals):
+                (cindxs,) = np.where(k == indices)
+                delta_logL = np.absolute(np.diff(self.data["logL"][cindxs]))
+                delta_logT = np.absolute(np.diff(self.data["logT"][cindxs]))
+                nindxs = [0]
+                cdelt_logL = 0.0
+                cdelt_logT = 0.0
+                for i in range(len(delta_logL)):
+                    cdelt_logL += delta_logL[i]
+                    cdelt_logT += delta_logT[i]
+                    if (cdelt_logL > logL_delta) or (cdelt_logT > logT_delta):
+                        nindxs.append(i)
+                        cdelt_logL = delta_logL[i]
+                        cdelt_logT = delta_logT[i]
+
+                if not max(nindxs) == len(delta_logL) - 1:
+                    nindxs.append(len(delta_logL) - 1)
+
+                for cname in self.data.keys():
+                    new_grid[cname] = np.concatenate(
+                        (new_grid[cname], self.data[cname][cindxs][nindxs])
+                    )
+
+            # update the grid
+            self.data = new_grid
 
 
 class ETParsec(EvolTracks):
@@ -240,7 +254,7 @@ class ETParsec(EvolTracks):
     EvolTracks specific to PARSEC calculations
     """
 
-    source = 'PARSEC'
+    source = "PARSEC"
 
     def load_orig_tables(self, files):
         """
@@ -257,22 +271,28 @@ def load_orig_tables(self, files):
         logL = np.array([])
         logT = np.array([])
         for cfile in files:
-            a = Table.read(cfile, format='ascii')
-            mass_act = np.concatenate((mass_act, a['MASS'].data))
-            logL = np.concatenate((logL, a['LOG_L'].data))
-            logT = np.concatenate((logT, a['LOG_TE'].data))
+            a = Table.read(cfile, format="ascii")
+            mass_act = np.concatenate((mass_act, a["MASS"].data))
+            logL = np.concatenate((logL, a["LOG_L"].data))
+            logT = np.concatenate((logT, a["LOG_TE"].data))
 
         self.data = {}
-        self.data['M_act'] = np.array(mass_act)
-        self.data['logL'] = np.array(logL)
-        self.data['logT'] = np.array(logT)
+        self.data["M_act"] = np.array(mass_act)
+        self.data["logL"] = np.array(logL)
+        self.data["logT"] = np.array(logT)
 
 
 class ETMist(EvolTracks):
     """
     MIST evolutionary tracks
     """
-    source = 'MIST'
+
+    def __init__(self):
+        super().__init()
+        self.source = "MIST"
+
+        files = f"{__ROOT__}/MIST/MIST_v1.2_feh_p0.00_afe_p0.0_vvcrit0.4_EEPS/*.eep"
+        self.load_orig_tables(files)
 
     def load_orig_tables(self, files):
         """
@@ -294,26 +314,25 @@ def load_orig_tables(self, files):
         phase = np.array([])
         eep = np.array([])
         for cfile in files:
-            a = Table.read(cfile, format='ascii', header_start=11)
-            tmass = a['star_mass'].data
+            a = Table.read(cfile, format="ascii", header_start=11)
+            tmass = a["star_mass"].data
             mass_act = np.concatenate((mass_act, tmass))
-            mass_ini = np.concatenate((mass_ini,
-                                       np.full((len(tmass)), max(tmass))))
-            logA = np.concatenate((logA, np.log10(a['star_age'].data)))
-            logL = np.concatenate((logL, a['log_L'].data))
-            logT = np.concatenate((logT, a['log_Teff'].data))
-            logg = np.concatenate((logg, a['log_g'].data))
-            phase = np.concatenate((phase, a['phase'].data))
+            mass_ini = np.concatenate((mass_ini, np.full((len(tmass)), max(tmass))))
+            logA = np.concatenate((logA, np.log10(a["star_age"].data)))
+            logL = np.concatenate((logL, a["log_L"].data))
+            logT = np.concatenate((logT, a["log_Teff"].data))
+            logg = np.concatenate((logg, a["log_g"].data))
+            phase = np.concatenate((phase, a["phase"].data))
             eep = np.concatenate((eep, range(len(a))))
 
         self.data = {}
-        self.data['M_act'] = np.log10(mass_act)
-        self.data['M_ini'] = np.log10(mass_ini)
-        self.data['logA'] = logA
-        self.data['logL'] = logL
-        self.data['logT'] = logT
-        self.data['logg'] = logg
-        self.data['phase'] = phase
-        self.data['eep'] = eep
-
-        self.alabels['eep'] = 'EEP'
+        self.data["M_act"] = np.log10(mass_act)
+        self.data["M_ini"] = np.log10(mass_ini)
+        self.data["logA"] = logA
+        self.data["logL"] = logL
+        self.data["logT"] = logT
+        self.data["logg"] = logg
+        self.data["phase"] = phase
+        self.data["eep"] = eep
+
+        self.alabels["eep"] = "EEP"

From 6ede39db066b86bc0cbbb218a70cdc9b2cb24a87 Mon Sep 17 00:00:00 2001
From: Karl Gordon <kgordon@stsci.edu>
Date: Thu, 21 Nov 2019 14:22:10 -0500
Subject: [PATCH 3/6] code to convert MIST ascii to FITS

---
 .../stars/ET_misc/conv_mist_ascii_to_fits.py  |  74 ++++++++++++
 beast/physicsmodel/stars/evoltracks.py        | 110 +++++++++---------
 2 files changed, 131 insertions(+), 53 deletions(-)
 create mode 100644 beast/physicsmodel/stars/ET_misc/conv_mist_ascii_to_fits.py

diff --git a/beast/physicsmodel/stars/ET_misc/conv_mist_ascii_to_fits.py b/beast/physicsmodel/stars/ET_misc/conv_mist_ascii_to_fits.py
new file mode 100644
index 000000000..fceabe517
--- /dev/null
+++ b/beast/physicsmodel/stars/ET_misc/conv_mist_ascii_to_fits.py
@@ -0,0 +1,74 @@
+"""
+Convert the MIST ascii files to FITS tables
+"""
+from glob import glob
+
+import numpy as np
+from tqdm import tqdm
+from astropy.table import Table, Column
+
+from beast.config import __ROOT__
+
+
+def convert_ascii_to_fits_one_met(infiles, FeH, outfilename):
+    """
+    Convert the ascii files into FITS tables to speed up repeated reading
+
+    Parameters
+    ----------
+    infiles : list
+        ascii files names (one per intial mass)
+    FeH : float
+        Fe/H value for this set of tracks
+    outfilename : str
+        output filename
+    """
+    if type(infiles) is not list:
+        files = [infiles]
+    else:
+        files = infiles
+
+    mass_act = np.array([])
+    mass_ini = np.array([])
+    logA = np.array([])
+    logL = np.array([])
+    logT = np.array([])
+    logg = np.array([])
+    phase = np.array([])
+    eep = np.array([])
+    for cfile in tqdm(files, desc=outfilename):
+        a = Table.read(cfile, format="ascii", header_start=11)
+        tmass = a["star_mass"].data
+        mass_act = np.concatenate((mass_act, tmass))
+        mass_ini = np.concatenate((mass_ini, np.full((len(tmass)), max(tmass))))
+        logA = np.concatenate((logA, np.log10(a["star_age"].data)))
+        logL = np.concatenate((logL, a["log_L"].data))
+        logT = np.concatenate((logT, a["log_Teff"].data))
+        logg = np.concatenate((logg, a["log_g"].data))
+        phase = np.concatenate((phase, a["phase"].data))
+        eep = np.concatenate((eep, range(len(a))))
+
+    data = Table()
+    data["M_act"] = Column(np.log10(mass_act))
+    data["M_ini"] = Column(np.log10(mass_ini))
+    data["logA"] = Column(logA)
+    data["logL"] = Column(logL)
+    data["logT"] = Column(logT)
+    data["logg"] = Column(logg)
+    data["phase"] = Column(phase)
+    data["eep"] = Column(eep)
+
+    met = np.full((len(eep)), FeH)
+    data["met"] = Column(met, description="Fe/H values")
+
+    data.write(outfilename, overwrite=True)
+
+
+if __name__ == "__main__":
+    FeH = [0.0, 0.25, 0.5]
+    FeH_str = ["p0.00", "p0.25", "p0.50"]
+
+    for cFeH, cFeH_str in zip(FeH, FeH_str):
+        cfiles = glob(f"{__ROOT__}/MIST/MIST_v1.2_feh_{cFeH_str}_afe_p0.0_vvcrit0.4_EEPS/*.eep")
+        outfile = f"MIST_FeH{cFeH:.2f}.fits"
+        convert_ascii_to_fits_one_met(cfiles, cFeH, outfile)
diff --git a/beast/physicsmodel/stars/evoltracks.py b/beast/physicsmodel/stars/evoltracks.py
index 840587781..96a81e905 100644
--- a/beast/physicsmodel/stars/evoltracks.py
+++ b/beast/physicsmodel/stars/evoltracks.py
@@ -4,7 +4,7 @@
 import numpy as np
 from scipy.interpolate import interp1d
 
-from astropy.table import Table
+from astropy.table import Table, vstack
 
 from beast.config import __ROOT__
 
@@ -49,6 +49,7 @@ def __init__(self):
             "phase": "evol phase",
             "M_act": "log(current mass)",
             "M_ini": "log(initial mass)",
+            "eep": "EEP",
         }
 
     def load_orig_tables(self, source):
@@ -249,50 +250,71 @@ def regrid(
             self.data = new_grid
 
 
-class ETParsec(EvolTracks):
+class ETMist(EvolTracks):
     """
-    EvolTracks specific to PARSEC calculations
+    MIST evolutionary tracks
     """
 
-    source = "PARSEC"
+    def __init__(self):
+        super().__init__()
+        self.source = "MIST"
+
+        self.orig_FeH = [0.0, 0.25, 0.5]
+        self.orig_files = ["MIST_FeH0.00.fits", "MIST_FeH0.25.fits", "MIST_FeH0.50.fits"]
+        # self.load_orig_tables(self.orig_files)
 
     def load_orig_tables(self, files):
         """
-        Read the tracks from  the original files
+        Read the tracks from the original files
 
+        Parameters
+        ----------
         files : str
             file or files with the evolutionary track calculations
-            often each file is for a single initial mass
+            often each for a single metallicity
+
+        Returns
+        -------
+        self.orig_tracks : astropy Table
+            Table with evolutionary track info as columns for all metallicities
         """
-        if not type(files) is list:
+        if type(files) is not list:
             files = [files]
 
-        mass_act = np.array([])
-        logL = np.array([])
-        logT = np.array([])
-        for cfile in files:
-            a = Table.read(cfile, format="ascii")
-            mass_act = np.concatenate((mass_act, a["MASS"].data))
-            logL = np.concatenate((logL, a["LOG_L"].data))
-            logT = np.concatenate((logT, a["LOG_TE"].data))
+        itables = [Table.read(cfile) for cfile in files]
+        if len(itables) > 1:
+            self.orig_tracks = vstack(itables)
+        else:
+            self.orig_tracks = itables[0]
 
-        self.data = {}
-        self.data["M_act"] = np.array(mass_act)
-        self.data["logL"] = np.array(logL)
-        self.data["logT"] = np.array(logT)
+    def get_evoltracks(self, masses, ages, metals=None, FeHs=None):
+        """
+        Get the evolutionary tracks for the specified ages, initial masses,
+        and metallicities.
 
+        Parameters
+        ----------
+        masses : list
+            Initial masses for grid
+        age : list
+            Ages for grid
+        metal, FeH : list
+            At least one needs to be set for the grid metallicities
 
-class ETMist(EvolTracks):
+        Returns
+        -------
+        type
+            Description of returned object.
+        """
+        pass
+
+
+class ETParsec(EvolTracks):
     """
-    MIST evolutionary tracks
+    EvolTracks specific to PARSEC calculations
     """
 
-    def __init__(self):
-        super().__init()
-        self.source = "MIST"
-
-        files = f"{__ROOT__}/MIST/MIST_v1.2_feh_p0.00_afe_p0.0_vvcrit0.4_EEPS/*.eep"
-        self.load_orig_tables(files)
+    source = "PARSEC"
 
     def load_orig_tables(self, files):
         """
@@ -302,37 +324,19 @@ def load_orig_tables(self, files):
             file or files with the evolutionary track calculations
             often each file is for a single initial mass
         """
-        if type(files) is not list:
+        if not type(files) is list:
             files = [files]
 
         mass_act = np.array([])
-        mass_ini = np.array([])
-        logA = np.array([])
         logL = np.array([])
         logT = np.array([])
-        logg = np.array([])
-        phase = np.array([])
-        eep = np.array([])
         for cfile in files:
-            a = Table.read(cfile, format="ascii", header_start=11)
-            tmass = a["star_mass"].data
-            mass_act = np.concatenate((mass_act, tmass))
-            mass_ini = np.concatenate((mass_ini, np.full((len(tmass)), max(tmass))))
-            logA = np.concatenate((logA, np.log10(a["star_age"].data)))
-            logL = np.concatenate((logL, a["log_L"].data))
-            logT = np.concatenate((logT, a["log_Teff"].data))
-            logg = np.concatenate((logg, a["log_g"].data))
-            phase = np.concatenate((phase, a["phase"].data))
-            eep = np.concatenate((eep, range(len(a))))
+            a = Table.read(cfile, format="ascii")
+            mass_act = np.concatenate((mass_act, a["MASS"].data))
+            logL = np.concatenate((logL, a["LOG_L"].data))
+            logT = np.concatenate((logT, a["LOG_TE"].data))
 
         self.data = {}
-        self.data["M_act"] = np.log10(mass_act)
-        self.data["M_ini"] = np.log10(mass_ini)
-        self.data["logA"] = logA
-        self.data["logL"] = logL
-        self.data["logT"] = logT
-        self.data["logg"] = logg
-        self.data["phase"] = phase
-        self.data["eep"] = eep
-
-        self.alabels["eep"] = "EEP"
+        self.data["M_act"] = np.array(mass_act)
+        self.data["logL"] = np.array(logL)
+        self.data["logT"] = np.array(logT)

From 7b7f741d6d86e002481c2cfed615faaab7e7e83f Mon Sep 17 00:00:00 2001
From: Karl Gordon <kgordon@stsci.edu>
Date: Thu, 21 Nov 2019 19:19:16 -0500
Subject: [PATCH 4/6] Reading orig files

---
 beast/physicsmodel/stars/evoltracks.py | 48 +++++++++++++-------------
 1 file changed, 24 insertions(+), 24 deletions(-)

diff --git a/beast/physicsmodel/stars/evoltracks.py b/beast/physicsmodel/stars/evoltracks.py
index 96a81e905..071343c61 100644
--- a/beast/physicsmodel/stars/evoltracks.py
+++ b/beast/physicsmodel/stars/evoltracks.py
@@ -139,7 +139,7 @@ def grid_metrics(self, check_keys=["logL", "logT", "logg"]):
 
         return metrics
 
-    def regrid(
+    def regrid_one_met(
         self,
         logmass_range=[-1.0, 2.0],
         logmass_delta=0.05,
@@ -148,11 +148,10 @@ def regrid(
         logT_delta=0.05,
     ):
         """
-        Interpolate a set of evolutionary tracks to a uniform grid
-        in log(initial mass) and variable grid in stellar age.
-        Use Equivalent Evolutionary Points (EEPs) values to do the
-        mass interpolation.  EEPs are provide as part of the evolutionary
-        tracks.
+        Interpolate a set of evolutionary tracks for a single metallicity
+        to a uniform grid in log(initial mass) and variable grid in stellar age.
+        Use Equivalent Evolutionary Points (EEPs) values to do the mass
+        interpolation.  EEPs are provide as part of the evolutionary tracks.
 
         Parameters
         ----------
@@ -260,34 +259,33 @@ def __init__(self):
         self.source = "MIST"
 
         self.orig_FeH = [0.0, 0.25, 0.5]
-        self.orig_files = ["MIST_FeH0.00.fits", "MIST_FeH0.25.fits", "MIST_FeH0.50.fits"]
-        # self.load_orig_tables(self.orig_files)
+        self.orig_files = [
+            f"{__ROOT__}MIST/MIST_FeH{cstr:.2f}.fits" for cstr in self.orig_FeH
+        ]
 
-    def load_orig_tables(self, files):
+    def get_orig_tables(self):
         """
         Read the tracks from the original files
 
-        Parameters
-        ----------
-        files : str
-            file or files with the evolutionary track calculations
-            often each for a single metallicity
-
         Returns
         -------
-        self.orig_tracks : astropy Table
+        orig_tracks : astropy Table
             Table with evolutionary track info as columns for all metallicities
         """
-        if type(files) is not list:
-            files = [files]
+        if isinstance(self.orig_files, list):
+            files = self.orig_files
+        else:
+            files = [self.orig_files]
 
         itables = [Table.read(cfile) for cfile in files]
         if len(itables) > 1:
-            self.orig_tracks = vstack(itables)
+            orig_tracks = vstack(itables)
         else:
-            self.orig_tracks = itables[0]
+            orig_tracks = itables[0]
+
+        return orig_tracks
 
-    def get_evoltracks(self, masses, ages, metals=None, FeHs=None):
+    def get_evoltracks(self, masses, metals=None, FeHs=None):
         """
         Get the evolutionary tracks for the specified ages, initial masses,
         and metallicities.
@@ -296,8 +294,6 @@ def get_evoltracks(self, masses, ages, metals=None, FeHs=None):
         ----------
         masses : list
             Initial masses for grid
-        age : list
-            Ages for grid
         metal, FeH : list
             At least one needs to be set for the grid metallicities
 
@@ -306,7 +302,11 @@ def get_evoltracks(self, masses, ages, metals=None, FeHs=None):
         type
             Description of returned object.
         """
-        pass
+        orig_tracks = self.get_orig_tracks()
+
+        # first interpolate for mass spacing
+
+        # then interpolate for metallicity spacing
 
 
 class ETParsec(EvolTracks):

From 98328d9f9982f7219520063206bf1e37ed7c0f69 Mon Sep 17 00:00:00 2001
From: "Karl D. Gordon" <kgordon@stsci.edu>
Date: Fri, 26 May 2023 12:03:55 -0400
Subject: [PATCH 5/6] better evoltrack plot and details

---
 beast/physicsmodel/stars/evoltracks.py | 22 ++++--
 beast/plotting/plot_evoltracks.py      | 92 ++++++++++++++++----------
 2 files changed, 74 insertions(+), 40 deletions(-)

diff --git a/beast/physicsmodel/stars/evoltracks.py b/beast/physicsmodel/stars/evoltracks.py
index 071343c61..9425851d8 100644
--- a/beast/physicsmodel/stars/evoltracks.py
+++ b/beast/physicsmodel/stars/evoltracks.py
@@ -58,7 +58,7 @@ def load_orig_tables(self, source):
         """
         print("not implemented")
 
-    def plot_tracks(self, ax, xval="logT", yval="logL", linestyle="-"):
+    def plot_tracks(self, ax, xval="logT", yval="logL", linestyle="-", color="k"):
         """
         Plot the tracks with the input x, y choices
 
@@ -75,6 +75,9 @@ def plot_tracks(self, ax, xval="logT", yval="logL", linestyle="-"):
 
         linestyle : string
             matplotlib linestyle
+
+        color : string
+            matplotlib color
         """
         if xval not in self.data.keys():
             raise ValueError("xval choice not in data table")
@@ -85,8 +88,11 @@ def plot_tracks(self, ax, xval="logT", yval="logL", linestyle="-"):
         uvals, indices = np.unique(self.data["M_ini"], return_inverse=True)
         for k, cval in enumerate(uvals):
             cindxs = np.where(k == indices)
+            # ax.plot(
+            #     self.data[xval][cindxs], self.data[yval][cindxs], linestyle=linestyle, color=color,
+            # )
             ax.plot(
-                self.data[xval][cindxs], self.data[yval][cindxs], linestyle=linestyle
+                self.data[xval][cindxs], self.data[yval][cindxs], "o", color=color, markersize=2
             )
 
         ax.set_xlabel(self.alabels[xval])
@@ -263,7 +269,7 @@ def __init__(self):
             f"{__ROOT__}MIST/MIST_FeH{cstr:.2f}.fits" for cstr in self.orig_FeH
         ]
 
-    def get_orig_tables(self):
+    def load_orig_tables(self, filename=None):
         """
         Read the tracks from the original files
 
@@ -272,10 +278,13 @@ def get_orig_tables(self):
         orig_tracks : astropy Table
             Table with evolutionary track info as columns for all metallicities
         """
-        if isinstance(self.orig_files, list):
-            files = self.orig_files
+        if filename is None:
+            if isinstance(self.orig_files, list):
+                files = self.orig_files
+            else:
+                files = [self.orig_files]
         else:
-            files = [self.orig_files]
+            files = filename
 
         itables = [Table.read(cfile) for cfile in files]
         if len(itables) > 1:
@@ -303,6 +312,7 @@ def get_evoltracks(self, masses, metals=None, FeHs=None):
             Description of returned object.
         """
         orig_tracks = self.get_orig_tracks()
+        self.data = orig_tracks
 
         # first interpolate for mass spacing
 
diff --git a/beast/plotting/plot_evoltracks.py b/beast/plotting/plot_evoltracks.py
index c167a3b10..aafe65aca 100644
--- a/beast/plotting/plot_evoltracks.py
+++ b/beast/plotting/plot_evoltracks.py
@@ -1,75 +1,99 @@
 """
 Make a plot of the evolutionary tracks
 """
-from __future__ import (absolute_import, division, print_function,
-                        unicode_literals)
-
 import matplotlib.pyplot as plt
 
-from beast.physicsmodel.stars.evoltracks import (ETParsec, ETMist)
+from beast.physicsmodel.stars.evoltracks import ETParsec, ETMist
 from beast.plotting.beastplotlib import initialize_parser
 
 
-if __name__ == '__main__':
+if __name__ == "__main__":
     parser = initialize_parser()
-    parser.add_argument('filename', type=str, nargs='*',
-                        help="name of file(s) with evolutionary track")
-    parser.add_argument('--type', default='mist', choices=['mist', 'parsec'],
-                        help="source of tracks")
+    parser.add_argument(
+        "filename", type=str, nargs="*", help="name of file(s) with evolutionary track"
+    )
+    parser.add_argument(
+        "--type", default="mist", choices=["mist", "parsec"], help="source of tracks"
+    )
+    parser.add_argument(
+        "--condense",
+        help="condense grid based on logM, logT, logL deltas",
+        action="store_true",
+    )
     args = parser.parse_args()
 
     fig, ax = plt.subplots(2, 2, figsize=(10, 10))
 
     # switch between the types of tracks
-    if args.type == 'parsec':
+    if args.type == "parsec":
         et = ETParsec()
     else:
         et = ETMist()
 
     # read in the tracks
-    et.load_orig_tables(args.filename)
+    et.data = et.load_orig_tables(filename=args.filename)
 
     orig_metrics = et.grid_metrics()
 
     # plot, plot, plot
-    nls = ':'
-    et.plot_tracks(ax[0, 0], xval='logT', yval='logL', linestyle=nls)
-    et.plot_tracks(ax[1, 1], xval='logA', yval='M_ini', linestyle=nls)
-    et.plot_tracks(ax[0, 1], xval='phase', yval='M_ini', linestyle=nls)
-    et.plot_tracks(ax[1, 0], xval='logT', yval='logg', linestyle=nls)
+    nls = ":"
+    color = "k"
+    et.plot_tracks(ax[0, 0], xval="logT", yval="logL", linestyle=nls, color=color)
+    et.plot_tracks(ax[1, 1], xval="logA", yval="M_ini", linestyle=nls, color=color)
+    et.plot_tracks(ax[0, 1], xval="eep", yval="M_ini", linestyle=nls, color=color)
+    et.plot_tracks(ax[1, 0], xval="logT", yval="logg", linestyle=nls, color=color)
 
     # regrid the evolutionary tracks to uniform log(mass) and variable age
-    print('size orig = ', len(et.data['M_ini']))
-
-    et.regrid(logmass_range=[-1., 3.],
-              logmass_delta=0.05,
-              logT_delta=0.05,
-              logL_delta=0.05)
-
-    print('size regrid = ', len(et.data['M_ini']))
+    print("size orig = ", len(et.data["M_ini"]))
+
+    logmass_range = [-1.0, 3.0]
+    condense = args.condense
+    logmass_delta = 0.05
+    logT_delta = 0.05
+    logL_delta = 0.05
+    et.regrid_one_met(
+        logmass_range=logmass_range,
+        condense=condense,
+        logmass_delta=logmass_delta,
+        logT_delta=logT_delta,
+        logL_delta=logL_delta,
+    )
+
+    print("size regrid = ", len(et.data["M_ini"]))
+
+    print("logM range:", logmass_range)
+    print("logM:", logmass_delta)
+
+    print("condense:", condense)
+    if condense:
+        print("condense deltas")
+        print("logT:", logT_delta)
+        print("logL:", logL_delta)
 
     regrid_metrics = et.grid_metrics()
+    print("grid metrics: diffs [min,  max, median, mean]")
     for ckey in regrid_metrics.keys():
         print(ckey)
-        print(orig_metrics[ckey])
-        print(regrid_metrics[ckey])
+        print("  orig:", orig_metrics[ckey])
+        print("regrid:", regrid_metrics[ckey])
 
     # get the new grid metrics
     # et.grid_metrics()
 
     # plot, plot, plot
-    nls = '-'
-    et.plot_tracks(ax[0, 0], xval='logT', yval='logL', linestyle=nls)
-    et.plot_tracks(ax[1, 1], xval='logA', yval='M_ini', linestyle=nls)
-    et.plot_tracks(ax[0, 1], xval='phase', yval='M_ini', linestyle=nls)
-    et.plot_tracks(ax[1, 0], xval='logT', yval='logg', linestyle=nls)
+    nls = "-"
+    color = "b"
+    et.plot_tracks(ax[0, 0], xval="logT", yval="logL", linestyle=nls, color=color)
+    et.plot_tracks(ax[1, 1], xval="logA", yval="M_ini", linestyle=nls, color=color)
+    et.plot_tracks(ax[0, 1], xval="eep", yval="M_ini", linestyle=nls, color=color)
+    et.plot_tracks(ax[1, 0], xval="logT", yval="logg", linestyle=nls, color=color)
 
     fig.tight_layout()
 
-    save_name = 'evoltracks'
+    save_name = "evoltracks"
     if args.tex:
-        plt.rc({'usetex': True})
+        plt.rc({"usetex": True})
     if args.savefig:
-        fig.savefig('{}.{}'.format(save_name, args.savefig))
+        fig.savefig("{}.{}".format(save_name, args.savefig))
     else:
         plt.show()

From 1e19c38d49728fd9aa1a04fbce443999d8bb8b2f Mon Sep 17 00:00:00 2001
From: "Karl D. Gordon" <kgordon@stsci.edu>
Date: Fri, 26 May 2023 16:39:35 -0400
Subject: [PATCH 6/6] partial progress to get evolutionary tracks fully
 regridded for physics model

---
 beast/config.py                               |   3 +
 beast/physicsmodel/model_grid.py              |  87 +++++++-------
 .../stars/ET_misc/conv_mist_ascii_to_fits.py  |  16 ++-
 beast/physicsmodel/stars/evoltracks.py        | 107 ++++++++++++------
 beast/plotting/plot_evoltracks.py             |   2 +-
 beast/tools/run/create_physicsmodel.py        |  11 +-
 beast/tools/verify_beast_settings.py          |  17 +++
 7 files changed, 160 insertions(+), 83 deletions(-)

diff --git a/beast/config.py b/beast/config.py
index ef919c2ea..3a5351392 100644
--- a/beast/config.py
+++ b/beast/config.py
@@ -13,6 +13,9 @@
 # __NTHREADS__ = 6
 __NTHREADS__ = 1
 
+# useful constant(s)
+solar_metalicity = 0.0153
+
 # library directory
 beast_envvar = "BEAST_LIBS"
 userhome = expanduser("~")
diff --git a/beast/physicsmodel/model_grid.py b/beast/physicsmodel/model_grid.py
index 626c0f8ae..80f9f17df 100644
--- a/beast/physicsmodel/model_grid.py
+++ b/beast/physicsmodel/model_grid.py
@@ -6,6 +6,7 @@
 from beast.physicsmodel.grid import SpectralGrid, SEDGrid
 from beast.physicsmodel import creategrid
 from beast.physicsmodel.stars import isochrone, stellib
+from beast.physicsmodel.stars import evoltracks
 from beast.physicsmodel.stars.isochrone import ezIsoch
 from beast.physicsmodel.dust import extinction
 from beast.physicsmodel.grid_and_prior_weights import (
@@ -14,50 +15,48 @@
 from beast.tools.beast_info import add_to_beast_info_file
 
 __all__ = [
-    "make_iso_table",
+    "make_evoltrack_table",
     "make_spectral_grid",
     "add_stellar_priors",
     "make_extinguished_sed_grid",
 ]
 
 
-def make_iso_table(
+def make_evoltrack_table(
     project,
-    oiso=None,
-    logtmin=6.0,
-    logtmax=10.13,
-    dlogt=0.05,
+    oet=None,
+    age_info=[6.0, 10.13, 0.05],
+    mass_info=[-1.0, 3.0, 0.05],
     z=[0.0152],
-    iso_fname=None,
+    et_fname=None,
     info_fname=None,
 ):
     """
-    The isochrone tables are loaded (downloading if necessary)
+    The evolutionary track table are loaded.
+    Determined based on evolutionary tracks directly or isochrones derived
+    from evolutionary tracks.
 
     Parameters
     ----------
     project : str
         project name
 
-    oiso : isochrone.Isochrone object
-        contains the full isochrones information
+    oet : evoltracks.EvolTracks or isochrone.Isochrone object
+        contains the full evolutionary track information
 
-    logtmin : float
-        log-age min
+    age_info : list
+        age information needed for isochrones [logtmin, logtmax, dlogt]
 
-    logtmax : float
-        log-age max
-
-    dlogt : float
-        log-age step to request
+    mass_info : list
+        mass information needed for evolutionary tracks [logmmin, logmmax, dlogm]
 
     z : float or sequence
         list of metalicity values, where default (Z=0.152) is adopted Z_sun
         for PARSEC/COLIBRI models
 
-    iso_fname : str
-        Set to specify the filename to save the isochrones to, otherwise
-        saved to project/project_iso.csv
+    et_fname : str
+        Set to specify the filename to save the gridded evolutionary tracks
+        to, otherwise saved to project/project_et.csv
 
     info_fname : str
         Set to specify the filename to save beast info to, otherwise
@@ -68,32 +67,44 @@ def make_iso_table(
     fname: str
        name of saved file
 
-    oiso: isochrone.Isochrone object
-        contains the full isochrones information
+    oet: evoltracks.EvolTracks or isochrone.Isochrone object
+        contains the full evolutionary track information
     """
-    if iso_fname is None:
-        iso_fname = "%s/%s_iso.csv" % (project, project)
-    if not os.path.isfile(iso_fname):
-        if oiso is None:
-            oiso = isochrone.PadovaWeb()
-
-        t = oiso._get_t_isochrones(max(5.0, logtmin), min(10.13, logtmax), dlogt, z)
-        t.header["NAME"] = "{0} Isochrones".format("_".join(iso_fname.split("_")[:-1]))
-        print("{0} Isochrones".format("_".join(iso_fname.split("_")[:-1])))
+    if et_fname is None:
+        et_fname = "%s/%s_et.csv" % (project, project)
+    if not os.path.isfile(et_fname):
+        if oet is None:
+            oet = isochrone.PadovaWeb()
+
+        if isinstance(oet, isochrone.Isochrone):
+            logtmin, logtmax, dlogt = age_info
+            t = oet._get_t_isochrones(max(5.0, logtmin), min(10.13, logtmax), dlogt, z)
+            t.header["NAME"] = "{0} Isochrones".format("_".join(et_fname.split("_")[:-1]))
+            print("{0} Isochrones".format("_".join(et_fname.split("_")[:-1])))
+            info = {"project": project, "logt_input": age_info, "z_input": z}
+            t.write(et_fname)
+            # maybe needed as ezIsoch is a proxy for a Table
+            # maybe we can just use a table????
+            oet = ezIsoch(et_fname)
+        elif isinstance(oet, evoltracks.EvolTracks):
+            tab = oet.get_evoltracks(mass_info, z)
+            tab.header["NAME"] = "{0} EvolTracks".format("_".join(et_fname.split("_")[:-1]))
+            print("{0} EvolTracks".format("_".join(et_fname.split("_")[:-1])))
+            info = {"project": project, "logm_input": mass_info, "z_input": z}
+        else:
+            print(f"Type {type(oet)} of evolutionary track not supported")
 
-        t.write(iso_fname)
+    else:
+        # read in the isochrone data from the file
+        #   not sure why this is needed, but reproduces previous ezpipe method
+        oet = ezIsoch(et_fname)
 
     # save info to the beast info file
-    info = {"project": project, "logt_input": [logtmin, logtmax, dlogt], "z_input": z}
     if info_fname is None:
         info_fname = f"{project}/{project}_beast_info.asdf"
     add_to_beast_info_file(info_fname, info)
 
-    # read in the isochrone data from the file
-    #   not sure why this is needed, but reproduces previous ezpipe method
-    oiso = ezIsoch(iso_fname)
-
-    return (iso_fname, oiso)
+    return (et_fname, oet)
 
 
 def make_spectral_grid(
diff --git a/beast/physicsmodel/stars/ET_misc/conv_mist_ascii_to_fits.py b/beast/physicsmodel/stars/ET_misc/conv_mist_ascii_to_fits.py
index fceabe517..23e374398 100644
--- a/beast/physicsmodel/stars/ET_misc/conv_mist_ascii_to_fits.py
+++ b/beast/physicsmodel/stars/ET_misc/conv_mist_ascii_to_fits.py
@@ -65,10 +65,18 @@ def convert_ascii_to_fits_one_met(infiles, FeH, outfilename):
 
 
 if __name__ == "__main__":
-    FeH = [0.0, 0.25, 0.5]
-    FeH_str = ["p0.00", "p0.25", "p0.50"]
+    # fmt: off
+    FeH = [-4.00, -3.50, -3.00, -2.50, -2.00, -1.75,
+           -1.50, -1.25, -1.00, -0.75, -0.25, 0.0,
+           0.25, 0.5]
+    FeH_str = ["m4.00", "m3.50", "m3.00", "m2.50", "m2.00", "m1.75",
+               "m1.50", "m1.25", "m1.00", "m0.75", "m0.25", "p0.00",
+               "p0.25", "p0.50"]
+    # fmt: on
 
     for cFeH, cFeH_str in zip(FeH, FeH_str):
-        cfiles = glob(f"{__ROOT__}/MIST/MIST_v1.2_feh_{cFeH_str}_afe_p0.0_vvcrit0.4_EEPS/*.eep")
-        outfile = f"MIST_FeH{cFeH:.2f}.fits"
+        cfiles = glob(
+            f"{__ROOT__}/MIST/MIST_v1.2_feh_{cFeH_str}_afe_p0.0_vvcrit0.4_EEPS/*.eep"
+        )
+        outfile = f"MIST_FeH{cFeH:.2f}_vvcrit0.4.fits"
         convert_ascii_to_fits_one_met(cfiles, cFeH, outfile)
diff --git a/beast/physicsmodel/stars/evoltracks.py b/beast/physicsmodel/stars/evoltracks.py
index 9425851d8..13286c201 100644
--- a/beast/physicsmodel/stars/evoltracks.py
+++ b/beast/physicsmodel/stars/evoltracks.py
@@ -6,7 +6,7 @@
 
 from astropy.table import Table, vstack
 
-from beast.config import __ROOT__
+from beast.config import __ROOT__, solar_metalicity
 
 
 __all__ = ["EvolTracks", "ETParsec", "ETMist"]
@@ -92,7 +92,11 @@ def plot_tracks(self, ax, xval="logT", yval="logL", linestyle="-", color="k"):
             #     self.data[xval][cindxs], self.data[yval][cindxs], linestyle=linestyle, color=color,
             # )
             ax.plot(
-                self.data[xval][cindxs], self.data[yval][cindxs], "o", color=color, markersize=2
+                self.data[xval][cindxs],
+                self.data[yval][cindxs],
+                "o",
+                color=color,
+                markersize=2,
             )
 
         ax.set_xlabel(self.alabels[xval])
@@ -145,8 +149,9 @@ def grid_metrics(self, check_keys=["logL", "logT", "logg"]):
 
         return metrics
 
+    @staticmethod
     def regrid_one_met(
-        self,
+        one_track,
         logmass_range=[-1.0, 2.0],
         logmass_delta=0.05,
         condense=False,
@@ -174,7 +179,7 @@ def regrid_one_met(
             default is 0.05 for all 3
         """
         # get the unique mass values
-        uvals, indices = np.unique(self.data["M_ini"], return_inverse=True)
+        uvals, indices = np.unique(one_track["M_ini"], return_inverse=True)
         # print(10**uvals)
 
         # ensure there are evolutionary tracks spanning
@@ -189,15 +194,15 @@ def regrid_one_met(
 
         # setup the new grid
         new_grid = {}
-        for cname in self.data.keys():
+        for cname in one_track.keys():
             new_grid[cname] = np.array([])
 
         # loop over eep values and interopolate to new mass grid
         # along constant eep tracks
-        uvals, indices = np.unique(self.data["eep"], return_inverse=True)
+        uvals, indices = np.unique(one_track["eep"], return_inverse=True)
         for k, cval in enumerate(uvals):
             (cindxs,) = np.where(k == indices)
-            cur_masses = self.data["M_ini"][cindxs]
+            cur_masses = one_track["M_ini"][cindxs]
 
             # only interpolate for masses defined for the current eep
             (new_gindxs,) = np.where(
@@ -206,32 +211,32 @@ def regrid_one_met(
                 )
             )
 
-            for cname in self.data.keys():
+            for cname in one_track.keys():
                 if cname == "eep":
                     vals = np.full((len(new_gindxs)), cval)
                 elif cname == "M_ini":
                     vals = new_mass_vals[new_gindxs]
                 else:
-                    f = interp1d(cur_masses, self.data[cname][cindxs])
+                    f = interp1d(cur_masses, one_track[cname][cindxs])
                     vals = f(new_mass_vals[new_gindxs])
 
                 new_grid[cname] = np.concatenate((new_grid[cname], vals))
 
         # update the grid
-        self.data = new_grid
+        one_track = new_grid
 
         # setup the condensed grid
         if condense:
             new_grid = {}
-            for cname in self.data.keys():
+            for cname in one_track.keys():
                 new_grid[cname] = np.array([])
 
             # loop over each mass track and condense
-            uvals, indices = np.unique(self.data["M_ini"], return_inverse=True)
+            uvals, indices = np.unique(one_track["M_ini"], return_inverse=True)
             for k, cval in enumerate(uvals):
                 (cindxs,) = np.where(k == indices)
-                delta_logL = np.absolute(np.diff(self.data["logL"][cindxs]))
-                delta_logT = np.absolute(np.diff(self.data["logT"][cindxs]))
+                delta_logL = np.absolute(np.diff(one_track["logL"][cindxs]))
+                delta_logT = np.absolute(np.diff(one_track["logT"][cindxs]))
                 nindxs = [0]
                 cdelt_logL = 0.0
                 cdelt_logT = 0.0
@@ -246,13 +251,15 @@ def regrid_one_met(
                 if not max(nindxs) == len(delta_logL) - 1:
                     nindxs.append(len(delta_logL) - 1)
 
-                for cname in self.data.keys():
+                for cname in one_track.keys():
                     new_grid[cname] = np.concatenate(
-                        (new_grid[cname], self.data[cname][cindxs][nindxs])
+                        (new_grid[cname], one_track[cname][cindxs][nindxs])
                     )
 
             # update the grid
-            self.data = new_grid
+            one_track = new_grid
+
+        return one_track
 
 
 class ETMist(EvolTracks):
@@ -262,14 +269,26 @@ class ETMist(EvolTracks):
 
     def __init__(self):
         super().__init__()
+        self.name = "MIST EvolTracks"
         self.source = "MIST"
 
-        self.orig_FeH = [0.0, 0.25, 0.5]
+        # fmt: off
+        self.orig_FeH = np.array([-4.00, -3.50, -3.00, -2.50, -2.00, -1.75,
+                                  -1.50, -1.25, -1.00, -0.75, -0.25, 0.0,
+                                  0.25, 0.5])
+        # fmt: on
         self.orig_files = [
-            f"{__ROOT__}MIST/MIST_FeH{cstr:.2f}.fits" for cstr in self.orig_FeH
+            f"{__ROOT__}MIST/MIST_FeH{cstr:.2f}_vvcrit0.4.fits"
+            for cstr in self.orig_FeH
         ]
 
-    def load_orig_tables(self, filename=None):
+        self.logmass_range = np.log10(np.array([0.1, 300.0]))
+        self.z_range = (
+            10 ** (np.array([min(self.orig_FeH), max(self.orig_FeH)]))
+            * solar_metalicity
+        )
+
+    def load_orig_tables(self, FeH=None, filename=None):
         """
         Read the tracks from the original files
 
@@ -287,37 +306,57 @@ def load_orig_tables(self, filename=None):
             files = filename
 
         itables = [Table.read(cfile) for cfile in files]
-        if len(itables) > 1:
-            orig_tracks = vstack(itables)
-        else:
-            orig_tracks = itables[0]
 
-        return orig_tracks
+        return itables
 
-    def get_evoltracks(self, masses, metals=None, FeHs=None):
+    def get_evoltracks(
+        self, mass_info, metals, condense=False, logT_delta=0.05, logL_delta=0.05
+    ):
         """
         Get the evolutionary tracks for the specified ages, initial masses,
         and metallicities.
 
         Parameters
         ----------
-        masses : list
-            Initial masses for grid
-        metal, FeH : list
-            At least one needs to be set for the grid metallicities
+        mass_info : list
+            info for mass range [logminmass, logmaxmass, deltalogmass]
+        metals : list
+            metallicities
+        logL_delta, logT_delta : float
+            deltas for the condensed grid
+            default is 0.05 for both
 
         Returns
         -------
         type
             Description of returned object.
         """
-        orig_tracks = self.get_orig_tracks()
-        self.data = orig_tracks
-
-        # first interpolate for mass spacing
+        # determine the min/max metallilcities needed from original grid
+        # FeH = np.log10(np.array(metals) / solar_metalicity)
+        # 0.5 is make sure the orig FeH includes one metallicity beyond the values
+        # gvals = (self.orig_FeH >= min(FeH) - 0.5) & (self.orig_FeH <= max(FeH) + 0.5)
+
+        # get the as computed evolutionary tracks
+        orig_data = self.load_orig_tables()
+        print(orig_data[0])
+
+        # interpolate for requested mass spacing
+        for k, ctrack in enumerate(orig_data):
+            orig_data[k] = self.regrid_one_met(
+                ctrack,
+                logmass_range=mass_info[0:2],
+                logmass_delta=mass_info[2],
+                condense=condense,
+                logT_delta=logT_delta,
+                logL_delta=logL_delta,
+            )
+        print(len(orig_data[0]["M_act"]))
+        exit()
 
         # then interpolate for metallicity spacing
 
+        # condense grid if requested
+
 
 class ETParsec(EvolTracks):
     """
diff --git a/beast/plotting/plot_evoltracks.py b/beast/plotting/plot_evoltracks.py
index aafe65aca..0b38a17ca 100644
--- a/beast/plotting/plot_evoltracks.py
+++ b/beast/plotting/plot_evoltracks.py
@@ -46,7 +46,7 @@
     # regrid the evolutionary tracks to uniform log(mass) and variable age
     print("size orig = ", len(et.data["M_ini"]))
 
-    logmass_range = [-1.0, 3.0]
+    logmass_range = [-1.0, 2.47]
     condense = args.condense
     logmass_delta = 0.05
     logT_delta = 0.05
diff --git a/beast/tools/run/create_physicsmodel.py b/beast/tools/run/create_physicsmodel.py
index 125c45d00..ee54587c7 100644
--- a/beast/tools/run/create_physicsmodel.py
+++ b/beast/tools/run/create_physicsmodel.py
@@ -10,7 +10,7 @@
 # BEAST imports
 from beast.physicsmodel.create_project_dir import create_project_dir
 from beast.physicsmodel.model_grid import (
-    make_iso_table,
+    make_evoltrack_table,
     make_spectral_grid,
     add_stellar_priors,
     make_extinguished_sed_grid,
@@ -69,12 +69,11 @@ def create_physicsmodel(beast_settings_info, nsubs=1, nprocs=1, subset=[None, No
     create_project_dir(settings.project)
 
     # download and load the isochrones
-    (iso_fname, oiso) = make_iso_table(
+    (iso_fname, oiso) = make_evoltrack_table(
         settings.project,
-        oiso=settings.oiso,
-        logtmin=settings.logt[0],
-        logtmax=settings.logt[1],
-        dlogt=settings.logt[2],
+        oet=settings.oiso,
+        age_info=settings.logt,
+        mass_info=settings.logmass,
         z=settings.z,
     )
 
diff --git a/beast/tools/verify_beast_settings.py b/beast/tools/verify_beast_settings.py
index 3abd112c9..0f0504b01 100644
--- a/beast/tools/verify_beast_settings.py
+++ b/beast/tools/verify_beast_settings.py
@@ -164,6 +164,23 @@ def verify_input_format(settings):
             [1.0, 7.0],
             [0.0, 1.0],
         ]
+    if settings.oiso.name == "MIST EvolTracks":
+        print('Working on the MIST Evolutionary Tracks')
+        parameters[3] = settings.logmass
+        parameters_names[3] = "logmass"
+
+        parameters_limits = [
+            settings.oiso.z_range,
+            None,
+            None,
+            settings.oiso.logmass_range,
+            [0.0, inf],
+            [1.0, 7.0],
+            [0.0, 1.0],
+        ]
+    else:
+        print("setup needed for ", settings.oiso.name)
+        exit()
 
     for i, param_ in enumerate(parameters):
         verify_one_input_format(