Fix Merge head

Author: RyotaInagaki1

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
@@ -1697,6 +1699,9 @@ def isochrone(self, filePath='', age=1 * 10 ** 8.1, metallicity=0.0,
         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.
         """
@@ -1739,10 +1744,8 @@ def isochrone(self, filePath='', age=1 * 10 ** 8.1, metallicity=0.0,
         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(isWR)
         iso.add_column(currentMass)
@@ -1752,20 +1755,18 @@ def isochrone(self, filePath='', age=1 * 10 ** 8.1, metallicity=0.0,
         # 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)
+        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'])
-            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

spisea/imf/imf.py

@@ -113,7 +113,7 @@ def generate_cluster(self, totalMass, seed=None):
         # Generate output arrays.
         masses = np.array([], dtype=float)
         isMultiple = np.array([], dtype=bool)
-        compMasses = []
+        compMasses = np.array([], dtype=np.object)
         systemMasses = np.array([], dtype=float)
 
         # Loop through and add stars to the cluster until we get to
@@ -139,8 +139,8 @@ def generate_cluster(self, totalMass, seed=None):
 
             # Dealing with multiplicity
             if self._multi_props != None:
-                compMasses = np.empty((len(newMasses),), dtype=np.object)
-                compMasses.fill([])
+                newCompMasses = np.empty((len(newMasses),), dtype=np.object)
+                newCompMasses.fill([])
 
                 # Determine the multiplicity of every star
                 MF = self._multi_props.multiplicity_fraction(newMasses)
@@ -152,13 +152,14 @@ def generate_cluster(self, totalMass, seed=None):
                 newSystemMasses = newMasses.copy()
 
                 # Function to calculate multiple systems more efficiently
-                compMasses, newSystemMasses, newIsMultiple = self.calc_multi(newMasses, compMasses,
-                                                                             newSystemMasses, newIsMultiple,
-                                                                             CSF, MF)
+                newCompMasses, newSystemMasses, newIsMultiple = self.calc_multi(newMasses, newCompMasses,
+                                                                                newSystemMasses, newIsMultiple,
+                                                                                CSF, MF)
 
                 newTotalMassTally = newSystemMasses.sum()
                 isMultiple = np.append(isMultiple, newIsMultiple)
                 systemMasses = np.append(systemMasses, newSystemMasses)
+                compMasses = np.append(compMasses, newCompMasses)
             else:
                 newTotalMassTally = newMasses.sum()
 

spisea/imf/tests/test_imf.py

@@ -24,8 +24,13 @@ def test_generate_cluster():
     # range of the requested mass.
     assert np.abs(M_cl - sysMass.sum()) < 120.0
 
-    return
+    # Check that enough companions were generated.
+    # Should be greater than 25% of the stars with companions.
+    mdx = np.where(isMulti)[0]
+    for mm in mdx:
+        assert len(compMass[mm]) > 0
 
+    return
 
 def test_prim_power():
     from .. import imf