diff --git a/.travis.yml b/.travis.yml index 9a9ff10e..c3048310 100755 --- a/.travis.yml +++ b/.travis.yml @@ -17,8 +17,9 @@ addons: packages: - graphviz - -stages: +script: +- echo "skipping tests. This is a draft of code I want Professor Lu and Matt to take a look at. + stages: # Do the style check and a single test job, don't proceed if it fails - name: Initial tests # Test docs, astropy dev, and without optional dependencies diff --git a/IsochroneMakerReformatted.py b/IsochroneMakerReformatted.py new file mode 100644 index 00000000..e58687fb --- /dev/null +++ b/IsochroneMakerReformatted.py @@ -0,0 +1,299 @@ +import numpy as np +import pandas as pd +from astropy.io import fits +from astropy.table import Table +import os +import hoki +from hoki import load +import glob +import re + +"""Source: https://stackoverflow.com/questions/44369504/ +how-to-convert-entire-dataframe-values-to-float-in-pandas""" + +colsToElim = [x for x in range(96) if x not in [hoki.dummy_dict[y] + for y in hoki.dummy_dict.keys()]] +def merger_related(f): + # Code that eventually returns whether + # The given model (in form of dataframe) + # has the conditions necessary + # to be a merger model. + return False + +def find_mergers(f): + # returns when the merger occurs (first time) + # assuming that f is a merger model. + return False + + +def reformatter(bevo_dir): + """I will expedite the process of writing and saving files by + making the BPASS stellar input files binary fits tables This method + creates a reformatted FITS version of + all of the available BPASS data so that data can be loaded quickly + and processed in the extractor method.""" + + # The list of possible BPASS metallicities. + # HMG stands for Qasi-Chemically Homogeneous Evolution Models. + + metallicities=["zem5", "zem4", "z001", "z002", "z003", "z004", "z006", "z008", + "z010", "z014","z020","z030", "z040", "zem5_hmg", + "zem4_hmg", "z001_hmg", "z002_hmg", "z003_hmg", "z004_hmg"] + # I will use the Python HashSet since (assuming that strings + # have a good HashCode, adding to the HashSet should be fast + + setofAll = set() + for x in metallicity: + a = set(glob.glob(hoki.MODELS_PATH + "/NEWSINMODS/" + x + "/*")) + b = set(glob.glob(hoki.MODELS_PATH + "/NEWBINMODS/NEWBINMODS/" + x + "/*")) + c = set(glob.glob(hoki.MODELS_PATH + "/NEWBINMODS/NEWSECMODS/" + x + "_2/*")) + setofAll = setofAll.union(a, b, c) + + for x in setofAll: + + f = pd.read_csv(hoki.MODELS_PATH + x, sep = "\s+", + header = None) + f.apply(pd.to_numeric, errors = 'coerce') + f.transpose() + + # Drop all of the columns that are not defined by the + # HOKI package provided mapping of column name to + # column number. + f.drop(columns = [f.columns[x] for x in colsToElim], + inplace = True) + # Source: https://stackoverflow.com/questions/483666/ + # reverse-invert-a-dictionary-mapping + # I create a mapping from column NUMBER to + # column name + invmap = {u: v for v, u in hoki.dummy_dict.items()} + f.rename(columns = invmap, inplace = True) + f["model_type"] = [typ] * len(f.index) + astroTable = Table.from_pandas(f) + # We will be recalculating the photometry related stuff + # when we create the isochrone. + # Will remove photometry related columns + astroTable.remove_columns(['V-I', 'U', + 'B', 'V', 'R', 'I', 'J', 'H', 'K','u', + 'g', 'r', 'i', 'z', 'f300w', 'f336w', + 'f435w', 'f450w', 'f555w', 'f606w','f814w', + 'U2', 'B2', 'V2','R2', 'I2','J2','H2','K2', + 'u2','g2','r2','i2','z2','f300w2','f336w2', + 'f435w2','f450w2','f555w2','f606w2','f814w2', + 'Halpha','FUV','NUV']) + # Merger Related indicates whether the star is + # a model with a merger in it. + astroTable['Merger_Related'] = [False] * len(astroTable) + # Creates container directory for stellar models grouped by metallicity + # and whether they are single or binary related + # the latter included mergers, QHE Models + if not os.path.isdir(bevo_dir + "iso/" + metallicity + '/' + + "sin" + + "/"): + os.makedirs(bevo_dir + + 'iso/' + metallicity + + "/" + SinOrBin + "/") + if not os.path.isdir(bevo_dir + "iso/" + metallicity + '/' + + "bin" + + "/"): + os.makedirs(bevo_dir + + 'iso/' + metallicity + + "/" + "bin" + "/") + # Single star models that are not the hmg models go into the /sin + # directory. + if (x[0:10] == 'NEWSINMODS' and + not re.matches(r"NEWSINMODS/z[0-9]+_hmg/.*")): + astroTable.write(bevo_dir + 'iso/' + met + '/' + + "sin" + "/" + "FitsModels" + + x[16:] + ".fits", format='fits', + overwrite=True) + print(x[16:] + ".fits") + else: + # Models that count as secondary star with compact remnants + # go into the /bin subdirectory of iso. + if (x[11:21] == 'NEWSECMODS'): + astroTable.write(bevo_dir + 'iso/' + met + '/' + + "bin" + "/" + "FitsModels" + + x[29:] + "sec.fits", format='fits', + overwrite = True) + print(bevo_dir + 'iso/' + metallicity + '/' + + SinOrBin + "/" + "FitsModels" + + x[29:] +"sec.fits") + else: + # Check for merger related (models_type=0) models + # and mark them. + if (find_mergers(astroTable)): + astroTable['Merger_Related'] = [True] * len(astroTable) + # Models that count as secondary star with compact remnants + # go into the /bin subdirectory of iso. + astroTable.write(bevo_dir + 'iso/' + metallicity + '/' + + "bin" + "/" + "FitsModels" + + x[27:] + "bin.fits", format='fits', + overwrite = True) + print(x[27:] + ".fits") + + except FileNotFoundError: + print("File not Readable/Found") + notFound.append(x) + #The following is commented out since I have yet to delete this version yet. +""" still_to_do=set() + + + if (x=="zem5" or x=="zem4" or x=="z001" or x=="z002" or x=="z003" or x=="z004"): + c=set(c).union(set(glob.glob(hoki.MODELS_PATH+"/NEWSINMODS/"+x+"hmg"+"/*"))) + still_to_do=still_to_do.union(a, b, c) + still_to_do=still_to_do.difference(included) + lenModPath=len(HOKI.MODELS_PATH) + for x in still_to_do: + try: + f = pd.read_csv(x, sep="\s+", + header=None) + f.apply(pd.to_numeric, errors='coerce') + f.transpose() + f.drop(columns=[f.columns[x] for x in colsToElim], inplace=True) + included.add(hoki.MODELS_PATH+x) + # Source: https://stackoverflow.com/questions/483666/ + # reverse-invert-a-dictionary-mapping + invmap = {u: v for v, u in hoki.dummy_dict.items()} + f.rename(columns=invmap, inplace=True) + f["model_type"]=[typ]*len(f.index) + astroTable = Table.from_pandas(f) + astroTable.remove_columns(['?','modelimf','mixedimf', 'V-I', 'U', + 'B', 'V', 'R', 'I', 'J', 'H', 'K','u', + 'g', 'r', 'i', 'z', 'f300w', 'f336w', + 'f435w', 'f450w', 'f555w', 'f606w','f814w', + 'U2', 'B2', 'V2','R2', 'I2','J2','H2','K2', + 'u2','g2','r2','i2','z2','f300w2','f336w2', + 'f435w2','f450w2','f555w2','f606w2','f814w2', + 'Halpha','FUV','NUV']) + if not os.path.isdir(bevo_dir + "iso/" + metallicity + '/' + + SinOrBin + "/"): + os.makedirs(bevo_dir + 'iso/' + metallicity + + "/" + SinOrBin + "/") + if (x[lenModPath+0:lenModPath+10]=='NEWSINMODS'): + nameMatcher=re.match(re".*/NEWSINMODS/[a-z0-9]+hmg/.*") + if (bool(nameMatcher)): + astroTable["model_type"]=[4]*len(astroTable) + else: + astroTable["model_type"]=[-12]*len(astroTable) + astroTable.write(bevo_dir + 'iso/' + met + '/' + + SinOrBin + "/" + "FitsModels" + + x[lenModPath+16:] + ".fits", format='fits', + overwrite=True) + print(x[lenModPath+16:]+".fits") + else: + print(x[lenModPath+11:lenModPath+21]) + + if (x[lenModPath+11:lenModPath+21]=='NEWSECMODS'): + astroTable.write(bevo_dir + 'iso/' + met + '/' + + SinOrBin + "/" + "FitsModels" + + x[lenModPath+29:] + "sec.fits", format='fits', + overwrite=True) + print(bevo_dir + 'iso/' + metallicity + '/' + + SinOrBin + "/" + "FitsModels" + + x[lenModPath+29:] +"sec.fits") + # Added designation for secondary models that can be + # either regular or single star + astroTable["model_type"]=[-15]*len(astroTable) + + else: + astroTable.write(bevo_dir + 'iso/' + metallicity + '/' + SinOrBin + + "/" + "FitsModels" + + x[lenModPath+27:] +"bin.fits", format='fits', + overwrite=True) + # Added designation for secondary models that can be + # either regular or single star + astroTable["model_type"]=[-10]*len(f.index) + print(x[lenModPath+27:] + ".fits") + """ + + except FileNotFoundError: + print("File not Readable/Found") + notFound.append(x) + + """ The follwing method takes in the reformatted stellar model files listed in + the corresponding model-metallicity's cluster input files originally meant + for color magnitude diagram makers. This is suppoesed to create the isochrone + files out of which BPASS creates its isochrone objects. For each stellar model + file, we try to find the row in the stellar model representing a model with the + age closest to the inputted age and is within the specified margin error for + log(age in years). Then that row is concattenated with other similar rows from + different stellar model files to create the isochrone table. The age closest + method may get rid of some data but is meant to improve accuracy of the age of the + stellar model and improve computational efficiency to some extent.""" +def extractor(SinOrBin, model, age, metallicity, input_dir, bpass_evo_dir, + margin): + entries = glob.glob(input_dir + "/*") + print(len(entries)) + stopLen = 4 + bigOne = None; + initlMass = np.nan + initlMass2 = np.nan + for x in entries: + if x[0] != '.' and not os.path.isdir(x): + # To Obtain initial mass I will use Regex + # to extract the initial mass of a model. + try: + org = Table.read(x, format = 'fits') + org = org.to_pandas() + starIndex = "Not Applicable" + if (x[-8:-5] == 'bin' and org['model_type'][0] == 0): + starIndex = findLastRelMax(org) + f = org[np.where(np.abs(np.log10(org['age'])-age) <= margin)[0]] + if (len(f) != 0): + f=f[np.where(np.abs(f['age'] - 10 ** age) == + np.min(np.abs(f['age'] - 10 ** age)))[0]] + if len(f) != 0: + index = f.index + indexlen = len(index) + if initial: + initial = False + bigOne = f + # mass and mass2 stand for the initial masses + # of the primary and secondary respectively. + bigOne['mass'] = [initlMass] * indexlen + bigOne['mass2'] = [initlMass2] * indexlen + bigOne['single'] = [True] * indexlen + bigOne['mergered?'] = [False] * indexlen + if (x[-8:-5] == 'bin' or x[-8:-5] == 'sec'): + bigOne['single'] = [False] * indexlen + # I only need to check one row since + # whether a model is a merger model + # is the same for all rows of the model + if bigOne['Merger_Related'][0]: + # Find the row in which the model merges. + # Note that I pass in the original grid. + merge_pt = find_mergers(org) + bigOne['mergered?'] = [x>merge for x in bigOne.index] + print(bigOne.shape) + else: + f['mass'] = [initlMass] * indexlen + f['mass2'] = [initlMass2] * indexlen + f['single'] = [True] * indexlen + f['Merger Related']=[False] * indexlen + if (x[-8:-5] == 'bin'): + f['single'] = [False] * indexlen + if bigOne['Merger_Related'][0]: + merge_pt = find_mergers(org) + bigOne['mergered?'] = [x > merge for x + in bigOne.index] + bigOne = pd.concat([f, bigOne], axis = 0) + print(bigOne.shape) + except FileNotFoundError: + print('File Not Found/Readable') + + try: + if not isinstance(bigOne, type(None)) and not (bigOne['age'].empty): + bigOne = bigOne.apply(pd.to_numeric, errors = 'coerce') + reduced = Table.from_pandas(bigOne) + if not os.path.isdir(bpass_evo_dir + + 'iso/' + model + '/' + metallicity + '/'): + os.makedirs(bpass_evo_dir + + 'iso/'+ model + '/' + metallicity + '/') + print(reduced.columns) + reduced.write(bpass_evo_dir + 'iso/' + model+"/" + metallicity + + '/' + 'iso' + str(age) + SinOrBin + + '.fits', format = 'fits', overwrite = True) + print(bigOne.shape) + except IndexError: + print ('It looks like there are no stars in out input file' + + ' that have the specified age...') \ No newline at end of file diff --git a/README.md b/README.md index cbb42b51..7bcc2b70 100755 --- a/README.md +++ b/README.md @@ -1,3 +1,6 @@ +# Ryota Inagaki's fork of the SPISEA repository. +This is my repository for storing my SPISEA code additions, which I hope to add by the end of this academic year to the dev branch. It is subject to modification throughout the next several months as I complete the code. +This code will eventually be merged with the main SPISEA repository. **For the most recent version, please check the Interpolation branch.** # SPISEA: Stellar Population Interface for Stellar Evolution and Atmospheres SPISEA is an python package that generates single-age, single-metallicity diff --git a/spisea/evolution.py b/spisea/evolution.py index e3ddf23c..426a26f2 100755 --- a/spisea/evolution.py +++ b/spisea/evolution.py @@ -7,6 +7,8 @@ import pdb import warnings from astropy.table import Table, vstack, Column +import astropy.constants as cs +import astropy.units as un from scipy import interpolate import pylab as py from spisea.utils import objects @@ -1615,6 +1617,182 @@ def get_orig_pisa_isochrones(metallicity=0.015): return data + +#==============================# +# BPASS v2.2 models +#==============================# +#These models were published by Stanway and Elridge et. al in 2018 +#The code is based from the Merged Evolution model methods outlined in the . After all, I want some degree of consistency with regards to +#Our original isochrone files are the stellar models (single, binary, and secondary) listed in the BPASS input files from one of the stellar models for ages 10**6.0, 10**6.1, ... , 10**11.0 years. + +class BPASS(StellarEvolution): + + """ These models were published by Stanway and Elridge et. al in 2018 + Some of the code is based from the Merged Evolution model methods + by Hosek et. al. After all, I want some degree of consistency with + regards to how isochrones are made.Also, I do not want to reinvent + the wheel twice. + If one has not checked the IsochroneMakerReformatted.py file, each + isochrone file is made of star from stellar model that are closest + to specified input age and are within the specified margin of error + for log(Age in years). This is to be doubly sure that I am getting + the same age.""" + + def __init__(self): + + """ + Note: All of the possible mass list items include possible mass for + single star, binary primary and secondary. Recall that BPASS age bins + are 10**6.0, 10**6.1, ... , 10**11.0 yrs """ + + + self.age_list = [round(6.0 + x * 0.1, 1) for x in range(0, 51)] + self.z_list = [ + 10 ** -5, + 10 ** -4, + 10 ** -3, + 0.002, + 0.003, + 0.004, + 0.006, + 0.008, + 0.010, + 0.014, + 0.020, + 0.030, + 0.040, + ] + self.models_dir = models_dir + "/BPASS/v2.2/" + self.mass_list = [round(0.1, 1)] + [round(0.12 + x * 0.020, 2) + for x in range(95)] + self.mass_list = self.mass_list + [round(2.05 + x * 0.05, 2) + for x in range(20)] + self.mass_list = self.mass_list + [round(3.1 + 0.1 * x, 1) + for x in range(70)] + self.mass_list = self.mass_list + [11 + x for x in range(90)] \ + + [125 + 25 * x for x in range(8)] + [400] + [500] + self.z_solar = 0.020 + #The possible Metallicities of a BPASS Stellar model + self.z_file_map = { + 10 ** -5: 'zem5/', + 10 ** -4: 'zem4/', + 0.00300: 'z003/', + 0.00200: 'z002/', + 0.00300: 'z003/', + 0.00400: 'z004/', + 0.00600: 'z006/', + 0.00800: 'z008/', + 0.01000: 'z010/', + 0.01400: 'z014/', + 0.02000: 'z020/', + 0.03000: 'z030/', + 0.04000: 'z040/', + } + + def isochrone(self, filePath='', age=1 * 10 ** 8.1, metallicity=0.0, + binary='bin', model='imf135_300'): + """ + This function adds several necessary columns (e.g. isWR, logg, and + whether a star is WR) to the existing isochrone files' tables. + + Important note on convention: I will specify the phase of the star + as follows: + White dwarf -> 101 + All other stars (no neutron stars or black holes) -> -1 + Just to make sure, I do have word that neutron stars + and black holes are NOT INCLUDED in the BPASS + evolution grids. + If you are REALLY curious about black holes, try using TUI + with BPASS. + """ + print(filePath) + + oldmetallicity = metallicity + metallicity = self.z_solar * 10 ** metallicity + log_age = math.log10(age) + print(log_age) + if log_age < np.min(self.age_list) or log_age > np.max(self.age_list): + logger.error('Requested age {0} is out of bounds.'.format(log_age)) + if (metallicity < np.min(self.z_list) or + metallicity > np.max(self.z_list)): + logger.error('Requested metallicity {0} is out of bounds.'. + format(metallicity)) + + # Borrow from matt hosek's code: + # Look for the closest valid age to the input age + + iso_file = 'iso' + str(round(log_age, 1)) + binary\ + + '.fits' + # Find the closest valid metallicity to the given metallicity + closest_metallicity=min([x for x in self.z_file_map], key=lambda x: abs(metallicity-x)) + z_dir=self.z_file_map[closest_metallicity] + print(z_dir) + + # Use the full path to the desired isochrone file. + + full_iso_file = self.models_dir + 'iso/' + model+'/' + z_dir + iso_file + if filePath: + iso = Table.read(filePath, format='fits') + else: + iso = Table.read(full_iso_file, format='fits') + print(iso.columns) + # Create columns for isWR, logg, + # Phase. I will need to find how + # phase for individual stars will be added + # then I will fill it in with appropriate values. + + isWR = Column([False] * len(iso), name='isWR') + isWR2 = Column([False] * len(iso), name='isWR2') + colG = Column([0.0] * len(iso), name='logg') + colP = Column([-1] * len(iso), name='phase') + colP2 = Column([-1] * len(iso), name='phase2') + iso.add_column(colG) + iso.add_column(isWR) + iso.add_column(currentMass) + iso.add_column(colP) + iso.add_column(colTWR) + iso.add_column(isWR2) + # We may as well delete a for loop here + + iso['age']=np.log10(iso['age']) + iso['logg'] = np.log10(( iso['M1'] * cs.GM_sun / ((10**iso['log(R1)']*cs.R_sun) ** 2))*un.s*un.s/un.m) + iso['logg2'] = np.log10(( iso['M2'] * cs.GM_sun / ((10**iso['log(R2)']*cs.R_sun) ** 2))*un.s*un.s/un.m) + iso.rename_column('M1', 'mass_current') + for x in iso: + x['age'] = np.log10(x['age']) + + # Using Stanway and Elridge Criterion for calculating whether + # a star is a WR star or not for using the PotsDam Atmospheres + + x['isWR'] = x['X'] < 0.40 and x['log(T1)'] > 4.45 + x['isWR']=x['isWR'] and x['secondary'] + x['isWR2']=x['isWR'] + + # Calculated logg=log10(Mass*G/(Radius of the star)^2) + # I made one solar radius equal to 695508 * 10 ** 3 meters + # Source: + # I also make sure that the computer understands whether a star is a + #white dwarf is not a white dwarf. This may help when running the + #atmosphere models for the photometry. + if (x['logg']>6.9 and x['log(L1)'] < -1 and x['log(T1)']<1.5): + x['phase']=101 + if (x['logg2']>6.9 and x['log(L2)'] < -1 and x['log(T2)']<1.5): + x['phase2']=101 + + #Changing column name to + #Making sure that names of the columns are consistent with + # general format of isochrone. + + iso.rename_column('log(T1)', 'logT') + iso.rename_column('M2', 'mass_current2') + iso.rename_column('log(L1)', 'logL') + iso.rename_column('age', 'logAge') + iso.meta['log_age'] = age + iso.meta['metallicity_in'] = oldmetallicity + iso.meta['metallicity_act'] = np.log10(closest_metallicity / + self.z_solar) + return iso + class Isochrone(object): def __init__(self, log_age): self.log_age = log_age diff --git a/spisea/synthetic.py b/spisea/synthetic.py index 30a4bb52..daa5c9b3 100755 --- a/spisea/synthetic.py +++ b/spisea/synthetic.py @@ -859,6 +859,305 @@ def plot_mass_luminosity(self, savefile=None): return +class Isochrone_Binary(Isochrone): + """ + Base Isochrone class. + + Parameters + ---------- + logAge : float + The age of the isochrone, in log(years) + AKs : float + The total extinction in Ks filter, in magnitudes + distance : float + The distance of the isochrone, in pc + metallicity : float, optional + The metallicity of the isochrone, in [M/H]. + Default is 0. + evo_model: model evolution class, optional + Set the stellar evolution model class. + Default is evolution.MISTv1(). + atm_func: model atmosphere function, optional + Set the stellar atmosphere models for the stars. + Default is get_merged_atmosphere. + wd_atm_func: white dwarf model atmosphere function, optional + Set the stellar atmosphere models for the white dwafs. + Default is get_wd_atmosphere + mass_sampling : int, optional + Sample the raw isochrone every `mass_sampling` steps. The default + is mass_sampling = 0, which is the native isochrone mass sampling + of the evolution model. + wave_range : list, optional + length=2 list with the wavelength min/max of the final spectra. + Units are Angstroms. Default is [3000, 52000]. + min_mass : float or None, optional + If float, defines the minimum mass in the isochrone. + Unit is solar masses. Default is None + max_mass : float or None, optional + If float, defines the maxmimum mass in the isochrone. + Units is solar masses. Default is None. + rebin : boolean, optional + If true, rebins the atmospheres so that they are the same + resolution as the Castelli+04 atmospheres. Default is False, + which is often sufficient synthetic photometry in most cases. + """ + + def __init__(self, logAge, AKs, distance, metallicity=0.0, + evo_model=default_evo_model, atm_func=default_atm_func, + wd_atm_func = default_wd_atm_func, + red_law=default_red_law, mass_sampling=1, + wave_range=[3000, 52000], min_mass=None, max_mass=None, + filters=['ubv,U', 'ubv,V', 'ubv,B', 'ubv,R', 'ubv,I'], + rebin=True): + + + c = constants + # Changes by Ryota: make atm_func and wd_atm_func instance vars + self.atm_func = atm_func + self.wd_atm_func = wd_atm_func + self.distance = distance + self.wave_range=wave_range + self.AKs=AKs + self.red_law=red_law + self.filters=filters + # Assert that the wavelength ranges are within the limits of the + # VEGA model (0.1 - 10 microns) + try: + assert wave_range[0] > 1000 + assert wave_range[1] < 100000 + except: + print('Desired wavelength range invalid. Limit to 1000 - 10000 A') + return + + # Get solar metallicity models for a population at a specific age. + # Takes about 0.1 seconds. + evol = evo_model.isochrone(age=10**logAge, + metallicity=metallicity) + + # Eliminate cases where log g is less than 0 + idx = np.where(evol['logg'] > 0) + evol = evol[idx] + + # Trim to desired mass range + if min_mass != None: + idx = np.where(evol['mass'] >= min_mass) + evol = evol[idx] + if max_mass != None: + idx = np.where(evol['mass'] <= max_mass) + evol = evol[idx] + + # Trim down the table by selecting every Nth point where + # N = mass sampling factor. + evol = evol[::mass_sampling] + + # Give luminosity, temperature, mass, radius units (astropy units). + L_all = 10 ** evol['logL'] * c.L_sun # luminsoity in W + T_all = 10 ** evol['logT'] * units.K + # TO DO: Conditionals to make sure I am using valid values + + R_all = np.sqrt(L_all / (4.0 * math.pi * c.sigma_sb * T_all**4)) + mass_all = evol['mass'] * units.Msun # masses in solar masses + logg_all = evol['logg'] # in cgs + mass_curr_all = evol['mass_current'] * units.Msun + phase_all = evol['phase'] + phase_all2 = evol['phase2'] + isWR_all = evol['isWR'] + isWR2=evol['isWR2] + Tef2 = 10 ** evol['log(T2)'] * units.K + R2 = 10 ** evol['log(R2)'] * c.R_sun + L2 = 10 ** evol['log(L2)'] * c.L_sun + singles = Table([mass_all, L_all, T_all, R_all, logg_all, isWR_all, + mass_curr_all, phase_all, evol['single']], + names = ['mass','log_a','L', 'Teff', 'R', 'gravity', + 'isWR', 'mass_current', 'phase', 'single']) + primaries = Table([mass_all, np.log10(c.R_sun * 10 ** evol['log(a)'] * + 1 / c.au), + L_all, T_all, R_all, logg_all, isWR_all, + mass_curr_all, phase_all2, evol['single']], + names = ['mass', 'log_a', 'L', 'Teff', 'R', 'gravity', + 'isWR', 'mass_current', 'phase', 'single']) + + secondaries = Table([evol['mass2'] * units.Msun, + np.log10(c.R_sun * + 10 ** evol['log(a)'] * + 1 / c.au), + L2, Tef2, R2, + evol['logg2'], isWR2, + evol['mass_current2'], + phase_all2, + evol['single'], + evol['mergered']], + names = ['mass','log_a', 'L', 'Teff', 'R', + 'gravity', 'isWR', 'mass_current2', + 'phase','single', 'merged']) + secondaries = secondaries[np.where(~secondaries['single'])[0]] + secondaries.remove_column('single') + singles = singles[np.where(singles['single'])[0]] + singles.remove_column('single') + primaries = p rimaries[np.where(~primaries['single'])[0]] + primaries.remove_column('single') + + + # I try to make sure that we have a queue + # of atm_function results to process + # If we have null values + self.spec_list_si = [] # For single Stars + self.spec_list2_pri = [] # For primary stars + self.spec_list3_sec = [] # For secondary stars + # Turns into an attribute since we will access this in another function + self.pairings2 = {"Singles": self.spec_list_si, + "Primaries": self.spec_list2_pri, + "Secondaries": self.spec_list3_sec} + pairings={"Singles": singles, "Primaries": primaries + , "Secondaries": secondaries} + self.pairings=pairings + self.pairings2=pairings2 + # For each temperature extract the synthetic photometry. + for x in pairings + tab = pairings[x] + atm_list = pairings2[x] + L_all = tab['L'] + T_all = tab['Teff'] + R_all = tab['R'] + logg_all = tab['gravity'] + gravity = tab['gravity'] + # For each temperature extract the synthetic photometry. + for ii in range(len(tab['Teff'])): + # Loop is currently taking about 0.11 s per iteration + + # Get the atmosphere model now. + # Wavelength is in Angstroms + # This is the time-intensive call... everything else is negligable. + # If source is a star, pull from star atmospheres. If it is a WD, + # pull from WD atmospheres + + if (tab[ii]['mergered?']==True and 'Secondaries'==x): + #For merged models, I make sure that if I end up reading the secondary + # of a model that has already been merged + # TO DO: Add info regarding this to the spec + tab[ii]['mass_current2'] = np.nan #The star is all gobbled up + tab[ii]['Teff'] = np.nan + tab[ii]['R'] = np.nan # No more secondary star with atm left + tab[ii]['L'] = np.nan # No more secondary star with atm left + tab[ii]['gravity']=np.nan + tab[ii]['isWR']=False # Star no longer exits + gravity = float( logg_all[ii] ) + L = float( L_all[ii].cgs / (units.erg / units.s)) # in erg/s + T = float( T_all[ii] / units.K) # in Kelvin + R = float( R_all[ii].to('pc') / units.pc) # in pc + phase = tab['phase'] + if (np.isfinite(gravity) and np.isfinite(L) and + np.isfinite(T) and + np.isfinite(R)): + if phase == 101: + star = wd_atm_func(temperature=T, gravity=gravity, + metallicity=metallicity, + verbose=False) + else: + star = atm_func(temperature=T, gravity=gravity, + metallicity=metallicity, + rebin=rebin) + + # Trim wavelength range down to JHKL range (0.5 - 5.2 microns) + star = spectrum.trimSpectrum(star, wave_range[0], + wave_range[1]) + + # Convert into flux observed at Earth (unreddened) + star *= (R / distance)**2 # in erg s^-1 cm^-2 A^-1 + + # Redden the spectrum. This doesn't take much time at all. + red = red_law.reddening(AKs).resample(star.wave) + star *= red + + # Save the final spectrum to our spec_list for later use. + atm_list.append(star) + else: + atm_list.append(None) + + self.singles = singles + self.primaries=primaries + self.secondaries = secondaries + + self.make_photometry() + # Append all the meta data to the summary table. + for tab in (singles, primaries, secondaries): + tab.meta['REDLAW'] = red_law.name + tab.meta['ATMFUNC'] = atm_func.__name__ + tab.meta['EVOMODEL'] = 'BPASS v2.2' + tab.meta['LOGAGE'] = logAge + tab.meta['AKS'] = AKs + tab.meta['DISTANCE'] = distance + tab.meta['METAL_IN'] = evol.meta['metallicity_in'] + tab.meta['METAL_ACT'] = evol.meta['metallicity_act'] + tab.meta['WAVEMIN'] = wave_range[0] + tab.meta['WAVEMAX'] = wave_range[1] + + + t2 = time.time() + print( 'Isochrone generation took {0:f} s.'.format(t2-t1)) + return + def make_photometry(self, rebin=True, vega=vega): + """ + Make synthetic photometry for the specified filters. This function + udpates the self.points table to include new columns with the + photometry. + + """ + startTime = time.time() + + meta = self.singles.meta + + print( 'Making photometry for isochrone: log(t) = %.2f AKs = %.2f dist = %d' % \ + (meta['LOGAGE'], meta['AKS'], meta['DISTANCE'])) + print( ' Starting at: ', datetime.datetime.now(), ' Usually takes ~5 minutes') + + # npoints = len(self.points) + verbose_fmt = 'M = {0:7.3f} Msun T = {1:5.0f} K m_{2:s} = {3:4.2f}' + + # Loop through the filters, get filter info, make photometry for + # all stars in this filter. + for ii in self.filters: + prt_fmt = 'Starting filter: {0:s} Elapsed time: {1:.2f} seconds' + print( prt_fmt.format(ii, time.time() - startTime)) + + filt = get_filter_info(ii, rebin=rebin, vega=vega) + filt_name = get_filter_col_name(ii) + + # Make the column to hold magnitudes in this filter. Add to points table. + col_name = 'm_' + filt_name + mag_col = Column(np.zeros(len(self.singles), dtype=float), name=col_name) + self.singles.add_column(mag_col) + mag_col = Column(np.zeros(len(self.secondaries), dtype=float), name=col_name) + self.secondaries.add_column(mag_col) + mag_col = Column(np.zeros(len(self.primaries), dtype=float), name=col_name) + self.primaries.add_column(mag_col) + + # Loop through each star in the isochrone and do the filter integration + print('Starting synthetic photometry') + for x in self.pairings: + listofStars=self.pairings2[x] + table=self.pairings[x] + length_of_list=len(listofStars) + for ss in range(listofStars): + star = listofStars[ss] + if (star!=None): + # These are already extincted, observed spectra. + star_mag = mag_in_filter(star, filt) + else: + star_mag = np.nan + + table[col_name][ss] = star_mag + + if (self.verbose and (ss % 100) == 0): + print( verbose_fmt.format(self.points['mass'][ss], + self.points['Teff'][ss], + filt_name, star_mag)) + + endTime = time.time() + print( ' Time taken: {0:.2f} seconds'.format(endTime - startTime)) + return + + class IsochronePhot(Isochrone): """ Make an isochrone with synthetic photometry in various filters. @@ -871,7 +1170,6 @@ class IsochronePhot(Isochrone): AKs : float The total extinction in Ks filter, in magnitudes - distance : float The distance of the isochrone, in pc