From 8617c031cf46b03ab3d3a6fa778ae1d065d4d05f Mon Sep 17 00:00:00 2001 From: Joe Filippazzo Date: Fri, 12 Aug 2022 14:33:32 -0400 Subject: [PATCH] Fixed mass and radius inference --- env/environment-3.7.yml | 2 +- env/environment-3.8.yml | 2 +- sedkit/__init__.py | 2 +- sedkit/relations.py | 6 ++--- sedkit/sed.py | 55 ++++++++++++++++++++++++----------------- setup.py | 2 +- 6 files changed, 40 insertions(+), 29 deletions(-) diff --git a/env/environment-3.7.yml b/env/environment-3.7.yml index 2441274e..a7172a9e 100644 --- a/env/environment-3.7.yml +++ b/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 diff --git a/env/environment-3.8.yml b/env/environment-3.8.yml index 26dd390f..b52bf2d5 100644 --- a/env/environment-3.8.yml +++ b/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 diff --git a/sedkit/__init__.py b/sedkit/__init__.py index f469225e..7cfd4443 100644 --- a/sedkit/__init__.py +++ b/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: diff --git a/sedkit/relations.py b/sedkit/relations.py index 4645de4f..0db98829 100644 --- a/sedkit/relations.py +++ b/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)) diff --git a/sedkit/sed.py b/sedkit/sed.py index f9f3d2c8..0f6ec5ee 100755 --- a/sedkit/sed.py +++ b/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,7 +1811,7 @@ 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 @@ -1816,13 +1819,16 @@ def infer_mass(self, mass_units=q.Msun, isochrone=False, plot=False): ---------- 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,7 +1868,7 @@ 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 @@ -1870,13 +1876,16 @@ def infer_radius(self, radius_units=q.Rsun, isochrone=False, plot=False): ---------- 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,31 +1906,31 @@ 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 @@ -1929,11 +1938,13 @@ def infer_Teff(self, teff_units=q.K, isochrone=False, plot=False): ---------- 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 diff --git a/setup.py b/setup.py index 3cf66757..f6c27115 100755 --- a/setup.py +++ b/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,