Merge pull request #74 from hover2pi/v1.2.3

Fixed mass and radius inference

env/environment-3.7.yml

@@ -23,4 +23,4 @@ dependencies:
   - pytest==7.1.1
   - sphinx==4.5.0
   - scipy==1.7.3
-  - svo-filters==0.4.1
+  - svo-filters==0.4.4

env/environment-3.8.yml

@@ -23,4 +23,4 @@ dependencies:
   - pytest==7.1.1
   - sphinx==4.5.0
   - scipy==1.8.0
-  - svo-filters==0.4.1
+  - svo-filters==0.4.4

sedkit/__init__.py

@@ -17,7 +17,7 @@
 __version_commit__ = ''
 _regex_git_hash = re.compile(r'.*\+g(\w+)')
 
-__version__ = '1.2.2'
+__version__ = '1.2.3'
 
 # from pkg_resources import get_distribution, DistributionNotFound
 # try:

sedkit/relations.py

@@ -240,12 +240,12 @@ def evaluate(self, rel_name, x_val, xunits=None, yunits=None, fit_local=False, p
                 out_yunits = full_rel['yunit'].to(yunits) * yunits if yunits is not None else full_rel['yunit'] or 1
 
                 # Use local points for relation
-                if isinstance(fit_local, int):
+                if isinstance(fit_local, int) and fit_local is not False:
 
                     # Trim relation data to nearby points, refit with low order polynomial, and evaluate
                     idx = np.argmin(np.abs(full_rel['x'] - (x_val[0] if isinstance(x_val, (list, tuple)) else x_val)))
                     x_min, x_max = full_rel['x'][max(0, idx - fit_local)], full_rel['x'][min(idx + fit_local, len(full_rel['x'])-1)]
-                    self.add_relation('mass(Lbol)', 2, xrange=[x_min, x_max], yunit=full_rel['yunit'], reject_outliers=True, plot=False)
+                    self.add_relation(rel_name, 2, xrange=[x_min, x_max], yunit=full_rel['yunit'], reject_outliers=True, plot=False)
                     rel = self.relations[rel_name]
 
                 # ... or use the full relation
@@ -285,7 +285,7 @@ def evaluate(self, rel_name, x_val, xunits=None, yunits=None, fit_local=False, p
                 else:
                     return y_val, self.ref
 
-            except ValueError as exc:
+            except IndexError as exc:
 
                 print(exc)
                 print("Could not evaluate the {} relation at {}".format(rel_name, x_val))

sedkit/sed.py

@@ -108,7 +108,7 @@ class SED:
     wave_units: astropy.units.quantity.Quantity
         The desired wavelength units
     """
-    def __init__(self, name='My Target', verbose=True, method_list=None, **kwargs):
+    def __init__(self, name='My Target', verbose=True, method_list=None, substellar=False, **kwargs):
         """
         Initialize an SED object
 
@@ -147,12 +147,13 @@ def __init__(self, name='My Target', verbose=True, method_list=None, **kwargs):
         self._evo_model = None
 
         # Static attributes
-        self.evo_model = 'DUSTY00'
+        self.evo_model = None
         self.frame = 'icrs'
         self.reddening = 0
         self.SpT = None
         self.multiple = False
         self.mainsequence = rel.DwarfSequence()
+        self.substellar = substellar
 
         # Dictionary to keep track of references
         self._refs = {}
@@ -281,7 +282,7 @@ def add_photometry(self, band, mag, mag_unc=None, system='Vega', ref=None, **kwa
             # Get the bandpass
             if isinstance(band, str):
                 bp = svo.Filter(band)
-            elif isinstance(band, svo.Filter):
+            elif isinstance(band, svo.Filter) or hasattr(band, 'wave_peak'):
                 bp, band = band, band.name
             else:
                 self.message('Not a recognized bandpass: {}'.format(band))
@@ -1048,6 +1049,8 @@ def evo_model(self, model):
 
             self._evo_model = iso.Isochrone(model, verbose=self.verbose)
 
+            self.message("Setting evolutionary model to {}".format(model))
+
         # Set as uncalculated
         self.calculated = False
 
@@ -1783,10 +1786,10 @@ def infer_logg(self, plot=False):
         """
         Estimate the surface gravity from model isochrones given an age and Lbol
         """
-        self.logg = None
+        self._logg = None
 
         # Try model isochrones
-        if self.age is not None and self.Lbol_sun is not None:
+        if self.evo_model is not None and self.age is not None and self.Lbol_sun is not None:
 
             # Default
             logg = None
@@ -1808,21 +1811,24 @@ def infer_logg(self, plot=False):
         else:
             self.message('Could not calculate logg without Lbol and age')
 
-    def infer_mass(self, mass_units=q.Msun, isochrone=False, plot=False):
+    def infer_mass(self, mass_units=q.Msun, plot=False):
         """
         Estimate the mass from model isochrones given an age and Lbol
 
         Parameters
         ----------
         mass_units: astropy.units.quantity.Quantity
             The units for the mass
-        isochrone: bool
-            Use the model isochrones
+        plot: bool
+            Plot the relation and the evaluated radius
         """
-        self.mass = None
+        self._mass = None
+
+        if self.substellar:
+            mass_units = q.Mjup
 
         # Try model isochrones
-        if isochrone and self.age is not None and self.Lbol_sun is not None:
+        if self.evo_model is not None and self.age is not None and self.Lbol_sun is not None:
 
             # Default
             mass = None
@@ -1862,21 +1868,24 @@ def infer_mass(self, mass_units=q.Msun, isochrone=False, plot=False):
                 self.message('Could not calculate mass without Lbol, M_2MASS.J, or M_2MASS.Ks')
 
 
-    def infer_radius(self, radius_units=q.Rsun, isochrone=False, plot=False):
+    def infer_radius(self, radius_units=q.Rsun, plot=False):
         """
         Estimate the radius from model isochrones given an age and Lbol
 
         Parameters
         ----------
         radius_units: astropy.units.quantity.Quantity
             The radius units
-        isochrone: bool
-            Use the model isochrones
+        plot: bool
+            Plot the relation and the evaluated radius
         """
-        self.radius = None
+        self._radius = None
+
+        if self.substellar:
+            radius_units = q.Rjup
 
         # Try model isochrones
-        if isochrone and (self.age is not None) and (self.Lbol_sun is not None):
+        if (self.evo_model is not None) and (self.age is not None) and (self.Lbol_sun is not None):
 
             # Default
             radius = None
@@ -1897,43 +1906,45 @@ def infer_radius(self, radius_units=q.Rsun, isochrone=False, plot=False):
             if self.radius is None and self.spectral_type is not None:
 
                 # Infer from Dwarf Sequence
-                self.radius = self.mainsequence.evaluate('radius(spt)', self.spectral_type, plot=plot)
+                self.radius = self.mainsequence.evaluate('radius(spt)', self.spectral_type, yunits=radius_units, fit_local=10, plot=plot)
 
             # Try radius(Lbol) relation
             elif self.radius is None and self.Lbol_sun is not None:
 
                 # Infer from Dwarf Sequence
-                self.radius = self.mainsequence.evaluate('radius(Lbol)', self.Lbol_sun, plot=plot)
+                self.radius = self.mainsequence.evaluate('radius(Lbol)', self.Lbol_sun, yunits=radius_units, fit_local=10, plot=plot)
 
             # Try radius(M_J) relation
             elif self.radius is None and self.get_mag('2MASS.J', 'abs') is not None:
 
                 # Infer from Dwarf Sequence
-                self.radius = self.mainsequence.evaluate('radius(M_J)', self.get_mag('2MASS.J'), plot=plot)
+                self.radius = self.mainsequence.evaluate('radius(M_J)', self.get_mag('2MASS.J'), yunits=radius_units, fit_local=10, plot=plot)
 
             # Try radius(M_Ks) relation
             elif self.radius is None and self.get_mag('2MASS.Ks', 'abs') is not None:
 
                 # Infer from Dwarf Sequence
-                self.radius = self.mainsequence.evaluate('radius(M_J)', self.get_mag('2MASS.Ks'), plot=plot)
+                self.radius = self.mainsequence.evaluate('radius(M_J)', self.get_mag('2MASS.Ks'), yunits=radius_units, fit_local=10, plot=plot)
 
             # No dice
             else:
                 self.message('Could not calculate radius without spectral_type, Lbol, M_2MASS.J, or M_2MASS.Ks')
 
-    def infer_Teff(self, teff_units=q.K, isochrone=False, plot=False):
+    def infer_Teff(self, teff_units=q.K, plot=False):
         """
         Infer the effective temperature
 
         Parameters
         ----------
         teff_units: astropy.units.quantity.Quantity
             The temperature units to use
+        plot: bool
+            Plot the relation and the evaluated radius
         """
-        self.Teff = None
+        self._Teff = None
 
         # Try model isochrones
-        if isochrone and self.age is not None and self.Lbol_sun is not None:
+        if self.evo_model is not None and self.age is not None and self.Lbol_sun is not None:
 
             # Default
             teff = None

setup.py

@@ -68,7 +68,7 @@ def run(self):
         'Topic :: Software Development :: Libraries :: Python Modules',
     ],
     packages=find_packages(exclude=["examples"]),
-    version='1.2.2',
+    version='1.2.3',
     setup_requires=['setuptools_scm'],
     install_requires=['numpy', 'astropy', 'bokeh', 'emcee', 'pysynphot', 'scipy', 'astroquery', 'dustmaps', 'pandas','svo_filters', 'healpy'],
     include_package_data=True,