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 new file mode 100644 index 000000000..23e374398 --- /dev/null +++ b/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) 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..13286c201 --- /dev/null +++ b/beast/physicsmodel/stars/evoltracks.py @@ -0,0 +1,391 @@ +# 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, vstack + +from beast.config import __ROOT__, solar_metalicity + + +__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 = "" + + # 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)", + "eep": "EEP", + } + + 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="-", color="k"): + """ + 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 + + color : string + matplotlib color + """ + 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, color=color, + # ) + ax.plot( + self.data[xval][cindxs], + self.data[yval][cindxs], + "o", + color=color, + markersize=2, + ) + + 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 + + @staticmethod + def regrid_one_met( + one_track, + 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 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 + ---------- + logmass_range : (float, float) + range of new mass grid + default is -1 to 2 (0.1 to 100 M_sun) + + 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(one_track["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 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(one_track["eep"], return_inverse=True) + for k, cval in enumerate(uvals): + (cindxs,) = np.where(k == indices) + cur_masses = one_track["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 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, one_track[cname][cindxs]) + vals = f(new_mass_vals[new_gindxs]) + + new_grid[cname] = np.concatenate((new_grid[cname], vals)) + + # update the grid + one_track = new_grid + + # setup the condensed grid + if condense: + new_grid = {} + for cname in one_track.keys(): + new_grid[cname] = np.array([]) + + # loop over each mass track and condense + 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(one_track["logL"][cindxs])) + delta_logT = np.absolute(np.diff(one_track["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 one_track.keys(): + new_grid[cname] = np.concatenate( + (new_grid[cname], one_track[cname][cindxs][nindxs]) + ) + + # update the grid + one_track = new_grid + + return one_track + + +class ETMist(EvolTracks): + """ + MIST evolutionary tracks + """ + + def __init__(self): + super().__init__() + self.name = "MIST EvolTracks" + self.source = "MIST" + + # 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}_vvcrit0.4.fits" + for cstr in self.orig_FeH + ] + + 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 + + Returns + ------- + orig_tracks : astropy Table + Table with evolutionary track info as columns for all metallicities + """ + if filename is None: + if isinstance(self.orig_files, list): + files = self.orig_files + else: + files = [self.orig_files] + else: + files = filename + + itables = [Table.read(cfile) for cfile in files] + + return itables + + 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 + ---------- + 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. + """ + # 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): + """ + 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) 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..0b38a17ca --- /dev/null +++ b/beast/plotting/plot_evoltracks.py @@ -0,0 +1,99 @@ +""" +Make a plot of the evolutionary tracks +""" +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" + ) + 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": + et = ETParsec() + else: + et = ETMist() + + # read in the tracks + et.data = et.load_orig_tables(filename=args.filename) + + orig_metrics = et.grid_metrics() + + # plot, plot, plot + 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"])) + + logmass_range = [-1.0, 2.47] + 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:", orig_metrics[ckey]) + print("regrid:", regrid_metrics[ckey]) + + # get the new grid metrics + # et.grid_metrics() + + # plot, plot, plot + 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" + 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() 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(