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 2ae85cde4..eedd64a72 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()