From 254c79898440a74c76f1d0d9174df2e7a94d9558 Mon Sep 17 00:00:00 2001 From: Ryota Inagaki <66146094+RyotaInagaki1@users.noreply.github.com> Date: Thu, 30 Jul 2020 14:53:14 -0700 Subject: [PATCH 01/15] BPASS Model Wrote out method for extracting BPASS evolutionary model isochrone and object. Still in development/debug phase. --- spisea/evolution.py | 83 +++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 83 insertions(+) diff --git a/spisea/evolution.py b/spisea/evolution.py index e3ddf23c..65d4cdb4 100755 --- a/spisea/evolution.py +++ b/spisea/evolution.py @@ -1615,6 +1615,89 @@ 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): + def __init__(self): + #Note: All of the possible mass list items include possible mass for single star, binary primary and secondary. + self.age_list=[round(6.0+x*0.1,1) for x in range(0, 51)] + #Possible metallicities for BPASS models + 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='/g/lu/models/evolution/BPASS/v2.2/' + #Possible masses + self.mass_list=[round(0.1, 1)]+[round(0.12+x*0.02, 2) for x in range(95)]+[round(2.05+x*0.05, 2) for x in range(20)]+[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.02 + self.z_file_map={10**-5: 'zem5/', 10**-4: 'zem4/', 0.003:'z003/',0.002: 'z002/', 0.003: 'z003/', 0.004: 'z004/', 0.006: + 'z006/',0.008: 'z008/', 0.010:'z010/',0.014:'z014/', 0.020: 'z020/', 0.030: 'z030/',0.040:'z040/' } + def isochrone(self, filePath="", age=1*10**8, metallicity=0.0, binary=True): + oldmetallicity=metallicity + #Convert log metallicity to metallicity with respect to solar masses + metallicity=self.z_solar*10**metallicity + log_age=math.log10(age) + if (log_agenp.max(self.age_list)): + logger.error('Requested age {0} is out of bounds.'.format(log_age)) + if (metallicitynp.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 + if log_age != self.age_list[0]: + age_idx = searchsorted(self.age_list, log_age, side='right') + else: + age_idx = searchsorted(self.age_list, log_age, side='left') + iso_file = "iso"+str(round(self.age_list[age_idx], 1))+".fits" + #Find the closest valid age to the given age + zindex=searchsorted(self.z_list, metallicity, side='left') + if (zindex==len(self.z_list)): + zindex=zindex-1 + z_dir=self.z_file_map[self.z_list[zindex]] + #Use the full path to the desired isochrone file. + full_iso_file=self.models_dir+'iso/'+z_dir+iso_file + print(z_dir) + print(full_iso_file) + if (filePath): + iso=Table.read(filePath, format='ascii') + else: + iso=Table.read(full_iso_file, format='ascii') + #Create columns for isWR, logg, Phase then fill it in with appropriate values. + isWR=Column([False]*len(iso), name='isWR') + colG=Column([0.0]*len(iso), name='logg') + currentMass=Column([0.0]*len(iso), name='mass_current') + colZ=Column([self.z_list[zindex]]*len(iso), name='Z') + colP=Column([-1]*len(iso), name='phase') + iso.add_column(colG) + iso.add_column(colZ) + iso.add_column(isWR) + iso.add_column(currentMass) + iso.add_column(colP) + for x in iso: + x['age']=math.log10(x['age']) + x['mass_current']=x['M1'] + #Using Stanway and Elridge Criterion for calculating whether a star is a WR star or not. + x['isWR']=(x['X']<0.40 and x['log(T1)']>4.0) + #Calculated logg=log10(Mass*G/(Radius of the star)^2) + x['logg']=math.log10(((6.67*10**-11)*x["M1"]*1.9891*10**30)/(((695508*10**3)*10**x['log(R1)'])**2)) + #Making sure that names of the columns are consistent with general format of isochrone. + iso.rename_column('log(T1)', 'logT') + iso.rename_column('log(L1)', 'logL') + iso.rename_column('M1', 'mass') + + iso.rename_column('age', 'logAge') + iso.meta['log_age'] = age + iso.meta['metallicity_in'] = oldmetallicity + iso.meta['metallicity_act'] = np.log10(self.z_list[zindex] / self.z_solar) + return iso + + + + + + class Isochrone(object): def __init__(self, log_age): self.log_age = log_age From c1a48d13bf177fc26bd71ccd289ef00e8cd59f1a Mon Sep 17 00:00:00 2001 From: Ryota Inagaki <66146094+RyotaInagaki1@users.noreply.github.com> Date: Thu, 30 Jul 2020 15:12:54 -0700 Subject: [PATCH 02/15] Version as of 3:12 PM 7/30/2020 --- spisea/evolution.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/spisea/evolution.py b/spisea/evolution.py index 65d4cdb4..06a85aa8 100755 --- a/spisea/evolution.py +++ b/spisea/evolution.py @@ -1677,12 +1677,14 @@ def isochrone(self, filePath="", age=1*10**8, metallicity=0.0, binary=True): iso.add_column(colP) for x in iso: x['age']=math.log10(x['age']) + #Would like to confirm whether M1 is the current mass of the stellar model at the given time. Just to be very very very sure. x['mass_current']=x['M1'] #Using Stanway and Elridge Criterion for calculating whether a star is a WR star or not. + #That means hydrogen mass fraction is less than 0.4 and T>10000 Kelvin x['isWR']=(x['X']<0.40 and x['log(T1)']>4.0) #Calculated logg=log10(Mass*G/(Radius of the star)^2) x['logg']=math.log10(((6.67*10**-11)*x["M1"]*1.9891*10**30)/(((695508*10**3)*10**x['log(R1)'])**2)) - #Making sure that names of the columns are consistent with general format of isochrone. + #Making sure that names of the columns are consistent with general format of isochrone. This is just so that the Synthetic.py works and can access necessary columns. iso.rename_column('log(T1)', 'logT') iso.rename_column('log(L1)', 'logL') iso.rename_column('M1', 'mass') From 80d45714c12d67319c3a41369da24f059d55dc19 Mon Sep 17 00:00:00 2001 From: RyotaInagaki1 Date: Fri, 21 Aug 2020 09:52:40 -0700 Subject: [PATCH 03/15] Evolution --- spisea/evolution.py | 214 +++++++++++++++++++++++++++++++------------- 1 file changed, 153 insertions(+), 61 deletions(-) diff --git a/spisea/evolution.py b/spisea/evolution.py index 06a85aa8..379cdce7 100755 --- a/spisea/evolution.py +++ b/spisea/evolution.py @@ -1623,82 +1623,174 @@ def get_orig_pisa_isochrones(metallicity=0.015): #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): +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. - self.age_list=[round(6.0+x*0.1,1) for x in range(0, 51)] - #Possible metallicities for BPASS models - 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='/g/lu/models/evolution/BPASS/v2.2/' - #Possible masses - self.mass_list=[round(0.1, 1)]+[round(0.12+x*0.02, 2) for x in range(95)]+[round(2.05+x*0.05, 2) for x in range(20)]+[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.02 - self.z_file_map={10**-5: 'zem5/', 10**-4: 'zem4/', 0.003:'z003/',0.002: 'z002/', 0.003: 'z003/', 0.004: 'z004/', 0.006: - 'z006/',0.008: 'z008/', 0.010:'z010/',0.014:'z014/', 0.020: 'z020/', 0.030: 'z030/',0.040:'z040/' } - def isochrone(self, filePath="", age=1*10**8, metallicity=0.0, binary=True): - oldmetallicity=metallicity - #Convert log metallicity to metallicity with respect to solar masses - metallicity=self.z_solar*10**metallicity - log_age=math.log10(age) - if (log_agenp.max(self.age_list)): + + """ + 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 + 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 (metallicitynp.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 - if log_age != self.age_list[0]: - age_idx = searchsorted(self.age_list, log_age, side='right') - else: - age_idx = searchsorted(self.age_list, log_age, side='left') - iso_file = "iso"+str(round(self.age_list[age_idx], 1))+".fits" - #Find the closest valid age to the given age - zindex=searchsorted(self.z_list, metallicity, side='left') - if (zindex==len(self.z_list)): - zindex=zindex-1 - z_dir=self.z_file_map[self.z_list[zindex]] - #Use the full path to the desired isochrone file. - full_iso_file=self.models_dir+'iso/'+z_dir+iso_file + 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) - print(full_iso_file) - if (filePath): - iso=Table.read(filePath, format='ascii') + + # 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='ascii') - #Create columns for isWR, logg, Phase then fill it in with appropriate values. - isWR=Column([False]*len(iso), name='isWR') - colG=Column([0.0]*len(iso), name='logg') - currentMass=Column([0.0]*len(iso), name='mass_current') - colZ=Column([self.z_list[zindex]]*len(iso), name='Z') - colP=Column([-1]*len(iso), name='phase') + 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') + currentMass = Column([0.0] * len(iso), name='mass_current') + colP = Column([-1] * len(iso), name='phase') + colP2 = Column([-1] * len(iso), name='phase2') + colTWR= Column([np.nan]*len(iso), name='logT_WR') iso.add_column(colG) - iso.add_column(colZ) 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 / ((cs.R_sun) ** 2))*un.s*un.s/un.m) + iso['logg2'] = np.log10(( iso['M2'] * cs.GM_sun / ((cs.R_sun) ** 2))*un.s*un.s/un.m) for x in iso: - x['age']=math.log10(x['age']) - #Would like to confirm whether M1 is the current mass of the stellar model at the given time. Just to be very very very sure. - x['mass_current']=x['M1'] - #Using Stanway and Elridge Criterion for calculating whether a star is a WR star or not. - #That means hydrogen mass fraction is less than 0.4 and T>10000 Kelvin - x['isWR']=(x['X']<0.40 and x['log(T1)']>4.0) - #Calculated logg=log10(Mass*G/(Radius of the star)^2) - x['logg']=math.log10(((6.67*10**-11)*x["M1"]*1.9891*10**30)/(((695508*10**3)*10**x['log(R1)'])**2)) - #Making sure that names of the columns are consistent with general format of isochrone. This is just so that the Synthetic.py works and can access necessary columns. + x['age'] = np.log10(x['age']) + x['mass_current'] = x['M1'] + + # 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'] + if x['isWR']: + x['logT_WR']=x['log(T1)'] + + # 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('M1', 'mass') - iso.rename_column('age', 'logAge') iso.meta['log_age'] = age iso.meta['metallicity_in'] = oldmetallicity - iso.meta['metallicity_act'] = np.log10(self.z_list[zindex] / self.z_solar) + iso.meta['metallicity_act'] = np.log10(closest_metallicity / + self.z_solar) return iso - - - - - class Isochrone(object): def __init__(self, log_age): From 5bf619f5823036abafd4445c5a7b0c5e414f9823 Mon Sep 17 00:00:00 2001 From: RyotaInagaki1 Date: Fri, 21 Aug 2020 09:57:21 -0700 Subject: [PATCH 04/15] Binary Isochrone Add --- spisea/synthetic.py | 286 ++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 286 insertions(+) diff --git a/spisea/synthetic.py b/spisea/synthetic.py index 30a4bb52..acae9663 100755 --- a/spisea/synthetic.py +++ b/spisea/synthetic.py @@ -859,6 +859,292 @@ 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): + + + t1 = time.time() + + 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 + # 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'] + 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 + tab = Table([mass_all, evol['mass2']*units.Msun,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_all, + evol['secondary'], evol['mass_current2'], + Tef2, R2, L2, evol['logg2'], evol['single']], + names=['mass','mass2','log_a','L', 'Teff', 'R', 'logg', 'isWR', 'mass_current', 'phase', 'secondary','mass_current2', 'Teff2', 'R2', 'L2', 'logg2', 'single']) + primaries = Table([mass_all, L_all, T_all, R_all, logg_all, isWR_all, mass_curr_all, phase_all, evol['single'], evol['secondary']], + names=['mass', 'L', 'Teff', 'R', 'logg', 'isWR', 'mass_current', 'phase', 'single', 'secondary']) + primaries=primaries[np.where(~primaries['secondary'])[0]] + primaries.remove_column('secondary') + + secondaries=Table([evol['mass2']*units.Msun,np.log10(c.R_sun*10**evol['log(a)']*1/c.au), + evol['secondary'], evol['mass_current2'], + Tef2, R2, L2, evol['logg2'], isWR2, evol['single']], + names=['mass2','log_a', 'secondary','mass_current2', 'Teff2', 'R2', 'L2', 'logg2', 'isWR2', 'single']) + secondaries=secondaries[np.where(~secondaries['single'])[0]] + secondaries.remove_column('single') + + + # Initialize output for stellar spectra + self.spec_list = [] + row_blacklist=[] + # For each temperature extract the synthetic photometry. + for ii in range(len(tab['Teff'])): + # Loop is currently taking about 0.11 s per iteration + 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 = phase_all[ii] + + # 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 (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. + self.spec_list.append(star) + else: + row_blacklist.append(ii) + row_blacklist2=[] + for ii in range(len(primaries['Teff'])): + # Loop is currently taking about 0.11 s per iteration + gravity = float( logg_all[ii] ) + L = float( primaires['L'][ii].cgs / (units.erg / units.s)) # in erg/s + T = float( primaries['Teff'][ii] / units.K) # in Kelvin + R = float( primaries['R'][ii].to('pc') / units.pc) # in pc + + # 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 not (np.isfinite(gravity) and np.isfinite(L) and np.isfinite(T) + and np.isfinite(R)): + row_blacklist2.append(ii) + row_blacklist3=[] + for ii in range(len(primaries['Teff'])): + # Loop is currently taking about 0.11 s per iteration + gravity = float( logg_all[ii] ) + L = float( secondaries['L2'][ii].cgs / (units.erg / units.s)) # in erg/s + T = float( secondaries['Teff2'][ii] / units.K) # in Kelvin + R = float( secondaries['R2'][ii].to('pc') / units.pc) # in pc + + # 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 not (np.isfinite(gravity) and np.isfinite(L) and np.isfinite(T) + and np.isfinite(R)): + row_blacklist3.append(ii) + + + tab.remove_rows(row_blacklist) + primaries.remove_rows(row_blacklist2) + secondaries.remove_rows(row_blacklist3) + self.points = tab + self.primaries=primaries + self.secondaries=secondaries + + self.make_photometry() + # Append all the meta data to the summary table. + tab.meta['REDLAW'] = red_law.name + tab.meta['ATMFUNC'] = atm_func.__name__ + tab.meta['EVOMODEL'] = type(evo_model).__name__ + tab.meta['LOGAGE'] = logAge + tab.meta['AKS'] = AKs + tab.meta['DISTANCE'] = distance + print(evol.meta) + 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.points.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(npoints, dtype=float), name=col_name) + self.points.add_column(mag_col) + + # Loop through each star in the isochrone and do the filter integration + print('Starting synthetic photometry') + for ss in range(npoints): + star = self.spec_list[ss] # These are already extincted, observed spectra. + star_mag = mag_in_filter(star, filt) + + self.points[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)) + self.primaries[col_name]=self.points[np.where(~self.points['secondary'])[0]][col_name] + self.secondaries[col_name]=self.points[np.where(~self.points['single'])[0]][col_name] + + endTime = time.time() + print( ' Time taken: {0:.2f} seconds'.format(endTime - startTime)) + + if self.save_file != None: + with warnings.catch_warnings(): + warnings.simplefilter('ignore') + self.points.write(self.save_file, overwrite=True) + return + + class IsochronePhot(Isochrone): """ Make an isochrone with synthetic photometry in various filters. From 7b3aeacf48261ee2a0b7fb1afe0e8c0526b7a754 Mon Sep 17 00:00:00 2001 From: RyotaInagaki1 Date: Thu, 10 Sep 2020 15:38:57 -0700 Subject: [PATCH 05/15] Reformatting --- .travis.yml | 5 +- IsochroneMakerReformatted.py | 299 +++++++++++++++++++++++++++++++++++ spisea/synthetic.py | 100 +++++++----- 3 files changed, 359 insertions(+), 45 deletions(-) create mode 100644 IsochroneMakerReformatted.py 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/spisea/synthetic.py b/spisea/synthetic.py index 5cd61303..daa5c9b3 100755 --- a/spisea/synthetic.py +++ b/spisea/synthetic.py @@ -910,8 +910,6 @@ def __init__(self, logAge, AKs, distance, metallicity=0.0, filters=['ubv,U', 'ubv,V', 'ubv,B', 'ubv,R', 'ubv,I'], rebin=True): - - t1 = time.time() c = constants # Changes by Ryota: make atm_func and wd_atm_func instance vars @@ -953,8 +951,8 @@ def __init__(self, logAge, AKs, distance, metallicity=0.0, 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 + 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)) @@ -962,61 +960,73 @@ def __init__(self, logAge, AKs, distance, metallicity=0.0, logg_all = evol['logg'] # in cgs mass_curr_all = evol['mass_current'] * units.Msun phase_all = evol['phase'] - phase_all = evol['phase2'] + 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 + 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', + 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), + 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', + 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) - , evol['mass_current2'], - Tef2, R2, L2, evol['logg2'], isWR2, evol['single'], evol['mergered']], - names=['mass','log_a', 'mass_current2', 'Teff', 'R', 'L', 'gravity', 'isWR', 'single']) - secondaries=secondaries[np.where(~secondaries['single'])[0]] + 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 = singles[np.where(singles['single'])[0]] singles.remove_column('single') - primaries=primaries[np.where(~primaries['single'])[0]] + 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 + # 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 + 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} + 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'] + 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 + # 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 @@ -1025,28 +1035,32 @@ def __init__(self, logAge, AKs, distance, metallicity=0.0, #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 left - tab[ii]['L']=np.nan # No more secondary star left + 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 = phase_all[ii] (Make sure phase can be added in) - if (np.isfinite(gravity) and np.isfinite(L) and np.isfinite(T) and + 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, + star = wd_atm_func(temperature=T, gravity=gravity, + metallicity=metallicity, verbose=False) else: - star = atm_func(temperature=T, gravity=gravity, metallicity=metallicity, + 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]) + 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 From 22e99bd098b0fcb08728d0099e22f628df6b96c6 Mon Sep 17 00:00:00 2001 From: Ryota Inagaki <66146094+RyotaInagaki1@users.noreply.github.com> Date: Wed, 30 Dec 2020 11:57:29 -0800 Subject: [PATCH 06/15] Update ReadMe, reflect that repo is in development --- README.md | 3 +++ 1 file changed, 3 insertions(+) diff --git a/README.md b/README.md index cbb42b51..ca6fd404 100755 --- a/README.md +++ b/README.md @@ -1,3 +1,6 @@ +# Ryota Inagaki's fork to the SPISEA repository. +This is my repository for storing my SPISEA code additions. 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. # SPISEA: Stellar Population Interface for Stellar Evolution and Atmospheres SPISEA is an python package that generates single-age, single-metallicity From 5d1c4e41ea8fdb946a41d30d3d90310c1b290346 Mon Sep 17 00:00:00 2001 From: Ryota Inagaki <66146094+RyotaInagaki1@users.noreply.github.com> Date: Sun, 3 Jan 2021 17:47:35 -0800 Subject: [PATCH 07/15] Update README.md --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index ca6fd404..f8dbb9ec 100755 --- a/README.md +++ b/README.md @@ -1,4 +1,4 @@ -# Ryota Inagaki's fork to the SPISEA repository. +# Ryota Inagaki's fork of the SPISEA repository. This is my repository for storing my SPISEA code additions. 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. # SPISEA: Stellar Population Interface for Stellar Evolution and Atmospheres From bc087ac838a803b2202e7dd67398e2057320b0e9 Mon Sep 17 00:00:00 2001 From: Ryota Inagaki <66146094+RyotaInagaki1@users.noreply.github.com> Date: Tue, 5 Jan 2021 18:34:24 -0800 Subject: [PATCH 08/15] Update README.md --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index f8dbb9ec..e8e75542 100755 --- a/README.md +++ b/README.md @@ -1,6 +1,6 @@ # Ryota Inagaki's fork of the SPISEA repository. This is my repository for storing my SPISEA code additions. 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. +This code will eventually be merged with the main SPISEA repository. **For Fall 2020 documentation and tutorials, please consult Fall 2020 Branch.** # SPISEA: Stellar Population Interface for Stellar Evolution and Atmospheres SPISEA is an python package that generates single-age, single-metallicity From 0aa7f3c939353d20999e74e71025d989c1e25050 Mon Sep 17 00:00:00 2001 From: Ryota Inagaki <66146094+RyotaInagaki1@users.noreply.github.com> Date: Wed, 13 Jan 2021 10:05:22 -0800 Subject: [PATCH 09/15] Update README.md --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index e8e75542..cde61aa8 100755 --- a/README.md +++ b/README.md @@ -1,6 +1,6 @@ # Ryota Inagaki's fork of the SPISEA repository. This is my repository for storing my SPISEA code additions. 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 Fall 2020 documentation and tutorials, please consult Fall 2020 Branch.** +This code will eventually be merged with the main SPISEA repository. **For Fall 2020 documentation and tutorials, please consult Fall-2020 Branch.** # SPISEA: Stellar Population Interface for Stellar Evolution and Atmospheres SPISEA is an python package that generates single-age, single-metallicity From 1bfab95774ed8c00e73ef87b155a40f22c94feca Mon Sep 17 00:00:00 2001 From: Ryota Inagaki <66146094+RyotaInagaki1@users.noreply.github.com> Date: Wed, 13 Jan 2021 10:06:28 -0800 Subject: [PATCH 10/15] Update README.md --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index cde61aa8..3dc0d855 100755 --- a/README.md +++ b/README.md @@ -1,6 +1,6 @@ # Ryota Inagaki's fork of the SPISEA repository. This is my repository for storing my SPISEA code additions. 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 Fall 2020 documentation and tutorials, please consult Fall-2020 Branch.** +This code will eventually be merged with the main SPISEA repository. **For Fall 2020 documentation and tutorials, please consult Fall-2020 Branch.** For the in-development Spring 2021 version, please consult the dev branch. # SPISEA: Stellar Population Interface for Stellar Evolution and Atmospheres SPISEA is an python package that generates single-age, single-metallicity From 801762bc5c836a3af7c5694b5e962c10601668f7 Mon Sep 17 00:00:00 2001 From: Ryota Inagaki <66146094+RyotaInagaki1@users.noreply.github.com> Date: Sat, 23 Jan 2021 19:40:28 -0800 Subject: [PATCH 11/15] Update README.md --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index 3dc0d855..b5baef59 100755 --- a/README.md +++ b/README.md @@ -1,5 +1,5 @@ # Ryota Inagaki's fork of the SPISEA repository. -This is my repository for storing my SPISEA code additions. It is subject to modification throughout the next several months as I complete the code. +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 Fall 2020 documentation and tutorials, please consult Fall-2020 Branch.** For the in-development Spring 2021 version, please consult the dev branch. # SPISEA: Stellar Population Interface for Stellar Evolution and Atmospheres From 57ee6302ff77caf2b547277ddda44534fbe28c40 Mon Sep 17 00:00:00 2001 From: Ryota Inagaki <66146094+RyotaInagaki1@users.noreply.github.com> Date: Tue, 17 Aug 2021 15:40:26 -0700 Subject: [PATCH 12/15] Update README.md --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index b5baef59..9754a579 100755 --- a/README.md +++ b/README.md @@ -1,6 +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 Fall 2020 documentation and tutorials, please consult Fall-2020 Branch.** For the in-development Spring 2021 version, please consult the dev branch. +This code will eventually be merged with the main SPISEA repository. **For Fall 2020 documentation and tutorials, please consult Fall-2020 Branch.** 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 From 154f75b5a453bdc4efacafb6eae74d0b75f2a6b7 Mon Sep 17 00:00:00 2001 From: Ryota Inagaki <66146094+RyotaInagaki1@users.noreply.github.com> Date: Tue, 17 Aug 2021 15:40:39 -0700 Subject: [PATCH 13/15] Update README.md --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index 9754a579..2ea0b669 100755 --- a/README.md +++ b/README.md @@ -1,6 +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 Fall 2020 documentation and tutorials, please consult Fall-2020 Branch.** For the most recent version, please check the Interpolation branch. +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 From 7c89b33af08840c6506c5c469385cefc8cb09dd2 Mon Sep 17 00:00:00 2001 From: Ryota Inagaki <66146094+RyotaInagaki1@users.noreply.github.com> Date: Tue, 17 Aug 2021 15:40:56 -0700 Subject: [PATCH 14/15] Update README.md --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index 2ea0b669..8c047ed1 100755 --- a/README.md +++ b/README.md @@ -1,6 +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. ** +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 From c35da72224387f9a8ccb5f260519d564d1b4f8c2 Mon Sep 17 00:00:00 2001 From: Ryota Inagaki <66146094+RyotaInagaki1@users.noreply.github.com> Date: Tue, 17 Aug 2021 15:41:06 -0700 Subject: [PATCH 15/15] Update README.md --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index 8c047ed1..7bcc2b70 100755 --- a/README.md +++ b/README.md @@ -1,6 +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.* +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