SynthPop compatibility updates

spisea/evolution.py

@@ -85,6 +85,7 @@ def __init__(self, model_dir, age_list, mass_list, z_list):
         self.z_list = z_list
         self.mass_list = mass_list
         self.age_list = age_list
+        self.model_version_name = "None"
 
         return
 
@@ -93,6 +94,7 @@ def __init__(self):
         r"""
         Define intrinsic properties for Geneva stellar models.
         """
+        self.model_version_name = "Geneva"
         # populate list of model masses (in solar masses)
         mass_list = [(0.1 + i*0.005) for i in range(181)]
 
@@ -171,6 +173,10 @@ class Ekstrom12(StellarEvolution):
         If true, then use rotating Ekstrom models. Default is true.
     """
     def __init__(self, rot=True):
+        if rot:
+            self.model_version_name = "Ekstrom12-rot"
+        else:
+            self.model_version_name = "Ekstrom12-norot"
         # define metallicity parameters for Ekstrom+12 models
         self.z_list = [0.014]
 
@@ -401,6 +407,7 @@ def __init__(self):
         models.
         """
         # populate list of model masses (in solar masses)
+        self.model_version_name = "Parsec1.2s"
         #mass_list = [(0.1 + i*0.005) for i in range(181)]
 
         # define metallicity parameters for Parsec models
@@ -551,6 +558,7 @@ def __init__(self):
         Define intrinsic properties for the Pisa (Tognelli+11) stellar
         models.
         """
+        self.model_version_name = "Pisa"
         # define metallicity parameters for Pisa models
         self.z_list = [0.015]
 
@@ -710,6 +718,7 @@ class Baraffe15(StellarEvolution):
     Downloaded from `BHAC15 site <http://perso.ens-lyon.fr/isabelle.baraffe/BHAC15dir/BHAC15_tracks>`_.
     """
     def __init__(self):
+        self.model_version_name = "Baraffe15"
         # define metallicity parameters for Baraffe models
         self.z_list = [0.015]
 
@@ -1009,7 +1018,7 @@ class MISTv1(StellarEvolution):
         was downloaded on 8/2018 (solar metallicity)
         and 4/2019 (other metallicities). Default is 1.2.
     """
-    def __init__(self, version=1.2):
+    def __init__(self, version=1.2, synthpop_extension=False):
         # define metallicity parameters for MIST models
         self.z_list = [0.0000014,   # [Fe/H] = -4.00
                        0.0000045,   # [Fe/H] = -3.50
@@ -1032,15 +1041,28 @@ def __init__(self, version=1.2):
 
         # Set version directory
         self.version = version
-        if self.version == 1.0:
+        self.synthpop_extension = synthpop_extension
+        if (self.version == 1.0) and (not synthpop_extension):
+            self.model_version_name = 'MISTv1.0'
             version_dir = 'v1.0/'
+        elif (self.version == 1.0) and synthpop_extension:
+            raise ValueError('Synthpop isochrone extension not supported for MISTv1.0 isochrones')
         elif self.version == 1.2:
+            self.model_version_name = 'MISTv1.2'
             version_dir = 'v1.2/'
         else:
             raise ValueError('Version {0} not supported for MIST isochrones'.format(version))
 
         # Specify location of model files
         self.model_dir = models_dir+'MISTv1/' + version_dir
+        if self.synthpop_extension:
+            self.model_version_name = self.model_version_name + '-synthpop'
+            self.model_extension_dir = models_dir+'MISTv1/' + version_dir[:-1] + '-synthpop/'
+            if not os.path.exists(self.model_extension_dir):
+                raise ValueError(f'Missing {self.model_extension_dir}. Please download the latest SPISEA data at https://w.astro.berkeley.edu/~jlu/spisea/spisea_models.tar.gz.')
+
+        else:
+            self.model_extension_dir = None
 
         # Specifying metallicity
         self.z_solar = 0.0142
@@ -1096,10 +1118,18 @@ def isochrone(self, age=1.e8, metallicity=0.0):
 
         # generate isochrone file string
         full_iso_file = self.model_dir + 'iso/' + z_dir + iso_file
+        if self.synthpop_extension:
+            addl_iso_file = self.model_extension_dir + 'iso/' + z_dir + iso_file
 
         # return isochrone data. Column locations depend on
         # version
         iso = Table.read(full_iso_file, format='fits')
+        if self.synthpop_extension:
+            addl_iso = Table.read(addl_iso_file, format='fits')
+            for row in addl_iso:
+                iso.add_row(row)
+            iso.sort('col3')
+            iso.meta['comments2'] = addl_iso.meta['comments']
         if self.version == 1.0:
             iso.rename_column('col7', 'Z')
             iso.rename_column('col2', 'logAge')
@@ -1232,6 +1262,10 @@ class MergedBaraffePisaEkstromParsec(StellarEvolution):
         If true, then use rotating Ekstrom models. Default is true.
     """
     def __init__(self, rot=True):
+        if rot:
+            self.model_version_name = "MergedBaraffePisaEkstromParsec-rot"
+        else:
+            self.model_version_name = "MergedBaraffePisaEkstromParsec-norot"
         # populate list of model masses (in solar masses)
         mass_list = [(0.1 + i*0.005) for i in range(181)]
 
@@ -1325,6 +1359,10 @@ class MergedPisaEkstromParsec(StellarEvolution):
         If true, then use rotating Ekstrom models. Default is true.
     """
     def __init__(self, rot=True):
+        if rot:
+            self.model_version_name = "MergedPisaEkstromParsec-rot"
+        else:
+            self.model_version_name = "MergedPisaEkstromParsec-norot"
         # populate list of model masses (in solar masses)
         mass_list = [(0.1 + i*0.005) for i in range(181)]
 
@@ -1422,6 +1460,7 @@ def __init__(self):
         Define intrinsic properties for merged Siess-meynetMaeder-Padova 
         stellar models.
         """
+        self.model_version_name = "MergedSiessGenevaPadova"
         # populate list of model masses (in solar masses)
         mass_list = [(0.1 + i*0.005) for i in range(181)]
 

spisea/synthetic.py

@@ -129,7 +129,7 @@ class ResolvedCluster(Cluster):
         True for verbose output.
     """
     def __init__(self, iso, imf, cluster_mass, ifmr=None, verbose=True,
-                     seed=None):
+                     seed=None, keep_low_mass_stars=False):
         Cluster.__init__(self, iso, imf, cluster_mass, ifmr=ifmr, verbose=verbose,
                              seed=seed)
         # Provide a user warning is random seed is set
@@ -163,7 +163,7 @@ def __init__(self, iso, imf, cluster_mass, ifmr=None, verbose=True,
 
         # Trim out bad systems; specifically, stars with masses outside those provided
         # by the model isochrone (except for compact objects).
-        star_systems, compMass = self._remove_bad_systems(star_systems, compMass)
+        star_systems, compMass = self._remove_bad_systems(star_systems, compMass, keep_low_mass_stars)
 
         ##### 
         # Make a table to contain all the information about companions.
@@ -433,13 +433,19 @@ def _make_companions_table(self, star_systems, compMass):
         if len(idx) != N_comp_tot and self.verbose:
             print( 'Found {0:d} companions out of stellar mass range'.format(N_comp_tot - len(idx)))
 
+        # For low-mass stars and substellar objects below isochrone, assume no mass loss and set phase to 98
+        low_mass_idxs = (companions['mass']<np.min(self.iso.points['mass']))
+        companions['mass_current'][low_mass_idxs] = companions['mass'][low_mass_idxs]
+        companions['phase'][low_mass_idxs] = 98
+
         # Double check that everything behaved properly.
-        assert companions['mass'][idx].min() > 0
+        if len(idx)>0:
+            assert companions['mass'][idx].min() > 0
 
         return companions
 
 
-    def _remove_bad_systems(self, star_systems, compMass):
+    def _remove_bad_systems(self, star_systems, compMass, keep_low_mass_stars):
         """
         Helper function to remove stars with masses outside the isochrone
         mass range from the cluster. These stars are identified by having 
@@ -454,16 +460,34 @@ def _remove_bad_systems(self, star_systems, compMass):
         # Convert nan_to_num to avoid errors on greater than, less than comparisons
         star_systems_teff_non_nan = np.nan_to_num(star_systems['Teff'], nan=-99)
         star_systems_phase_non_nan = np.nan_to_num(star_systems['phase'], nan=-99)
-        if self.ifmr == None:
+        if (self.ifmr == None) and (not keep_low_mass_stars):
+            print('Remove low mass stars below grid and compact objects')
             # Keep only those stars with Teff assigned.
             idx = np.where(star_systems_teff_non_nan > 0)[0]
-        else:
+        elif not keep_low_mass_stars:
+            print('Remove low mass stars, keep compact objects')
             # Keep stars (with Teff) and any other compact objects (with phase info). 
             idx = np.where( (star_systems_teff_non_nan > 0) | (star_systems_phase_non_nan >= 0) )[0]
+        elif self.ifmr == None:
+            print('Remove compact objects, keep low mass stars below grid')
+            # Keep stars (with Teff) and objects below mass grid
+            idx = np.where( (star_systems_teff_non_nan > 0) | ((star_systems['mass']<np.min(self.iso.points['mass']))) )[0]
+        else:
+            print('Keep low mass stars below grid and compact objects')
+            # Keep all
+            idx = np.where( (star_systems_teff_non_nan > 0) | (star_systems_phase_non_nan >= 0) |
+                            ((star_systems['mass']<np.min(self.iso.points['mass']))) )[0]
 
         if len(idx) != N_systems and self.verbose:
             print( 'Found {0:d} stars out of mass range'.format(N_systems - len(idx)))
 
+        if keep_low_mass_stars:
+            lm_idx = np.where(star_systems['mass']<np.min(self.iso.points['mass']))[0]
+            # Adjust the properties as needed
+            star_systems['mass_current'][lm_idx] = star_systems['mass'][lm_idx]
+            star_systems['phase'][lm_idx] = 98
+            #pdb.set_trace()
+
         star_systems = star_systems[idx]
         N_systems = len(star_systems)
 
@@ -511,10 +535,10 @@ class ResolvedClusterDiffRedden(ResolvedCluster):
         True for verbose output.
     """
     def __init__(self, iso, imf, cluster_mass, deltaAKs,
-                 ifmr=None, verbose=False, seed=None):
+                 ifmr=None, verbose=False, seed=None, keep_low_mass_stars=False):
 
         ResolvedCluster.__init__(self, iso, imf, cluster_mass, ifmr=ifmr, verbose=verbose,
-                                     seed=seed)
+                                     seed=seed, keep_low_mass_stars=keep_low_mass_stars)
 
         # Set random seed, if desired
         if seed is not None:
@@ -809,6 +833,7 @@ def __init__(self, logAge, AKs, distance, metallicity=0.0,
         tab.meta['REDLAW'] = red_law.name
         tab.meta['ATMFUNC'] = atm_func.__name__
         tab.meta['EVOMODEL'] = type(evo_model).__name__
+        tab.meta['EVOMODELVERSION'] = evo_model.model_version_name
         tab.meta['LOGAGE'] = logAge
         tab.meta['AKS'] = AKs
         tab.meta['DISTANCE'] = distance
@@ -1141,6 +1166,13 @@ def check_save_file(self, evo_model, atm_func, red_law):
                 (tmp.meta['ATMFUNC'] == atm_func.__name__) &
                  (tmp.meta['REDLAW'] == red_law.name) ):
                 out_bool = True
+
+            # Check model version if it was logged
+            if 'EVOMODELVERSION' in tmp.meta:
+                if tmp.meta['EVOMODELVERSION']!=evo_model.model_version_name:
+                    out_bool=False
+            else:
+                print(f"No version information found for evolution model {type(evo_model).__name__}.")
 
         return out_bool
 
@@ -1325,6 +1357,7 @@ def __init__(self, logAge, distance, evo_model=default_evo_model,
 
         tab.meta['ATMFUNC'] = atm_func.__name__
         tab.meta['EVOMODEL'] = type(evo_model).__name__
+        tab.meta['EVOMODELVERSION'] = evo_model.model_version_name
         tab.meta['LOGAGE'] = logAge
         tab.meta['DISTANCE'] = distance
         tab.meta['WAVEMIN'] = wave_range[0]

spisea/tests/test_synthetic.py

@@ -612,6 +612,15 @@ def test_cluster_mass():
     assert np.abs(cluster_mass_out - cluster_mass) < 200.0   # within 200 Msun of desired mass.
     print('Cluster Mass: IN = ', cluster_mass, " OUT = ", cluster_mass_out)
 
+    cluster3 = syn.ResolvedCluster(iso, my_imf1, cluster_mass, ifmr=my_ifmr, keep_low_mass_stars=True)
+    clust3 = cluster3.star_systems
+    print('Constructed cluster w/ keep_low_mass_stars: %d seconds' % (time.time() - startTime))
+
+    # Check that the total mass is within tolerance of input mass
+    cluster_mass_out = clust3['systemMass'].sum()
+    assert np.abs(cluster_mass_out - cluster_mass) < 200.0   # within 200 Msun of desired mass.
+    print('Cluster Mass: IN = ', cluster_mass, " OUT = ", cluster_mass_out)
+
     ##########
     # Test with multiplicity
     ##########