WIP: Evolutionary tracks for physicsgrid

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("~")

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(

beast/physicsmodel/stars/ET_misc/conv_mist_ascii_to_fits.py

@@ -0,0 +1,82 @@
+"""
+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__":
+    # 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}_vvcrit0.4.fits"
+        convert_ascii_to_fits_one_met(cfiles, cFeH, outfile)

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)