From f14a46683a604f5f3b3dda8d15d7426dc9d4a25a Mon Sep 17 00:00:00 2001 From: Tomas Stolker Date: Thu, 6 May 2021 11:26:11 +0200 Subject: [PATCH 1/7] Improved tick labels for logarithmic wavelengths, bug fix after refactoring of ReadModel --- species/data/database.py | 4 ++-- species/plot/plot_spectrum.py | 21 ++++++++++++++++----- species/read/read_model.py | 22 +++++++++++----------- 3 files changed, 29 insertions(+), 18 deletions(-) diff --git a/species/data/database.py b/species/data/database.py index dacb2b1a..53e9c0ae 100644 --- a/species/data/database.py +++ b/species/data/database.py @@ -919,7 +919,7 @@ def add_object(self, for i, hdu_item in enumerate(hdulist): data = np.asarray(hdu_item.data) - corr_warn = f'The covariance matrix from {value[1]} contains ' \ + corr_warn = f'The matrix from {value[1]} contains ' \ f'ones along the diagonal. Converting this ' \ f'correlation matrix into a covariance matrix.' @@ -956,7 +956,7 @@ def add_object(self, print(' - Covariance matrix:') if np.all(np.diag(data) == 1.): - warnings.warn(f'The covariance matrix from {value[1]} contains ones on ' + warnings.warn(f'The matrix from {value[1]} contains ones on ' f'the diagonal. Converting this correlation matrix into a ' f'covariance matrix.') diff --git a/species/plot/plot_spectrum.py b/species/plot/plot_spectrum.py index 191984f4..9de03171 100644 --- a/species/plot/plot_spectrum.py +++ b/species/plot/plot_spectrum.py @@ -14,7 +14,7 @@ import matplotlib.pyplot as plt from typeguard import typechecked -from matplotlib.ticker import AutoMinorLocator, MultipleLocator +from matplotlib.ticker import AutoMinorLocator, MultipleLocator, ScalarFormatter from species.core import box, constants from species.read import read_filter @@ -207,11 +207,13 @@ def plot_spectrum(boxes: list, # ax1.set_yticks([1e-5, 1e-4, 1e-3, 1e-2, 1e-1, 1e0]) # ax3.set_yticks([-2., 0., 2.]) - if filters is not None and scale[0] == 'linear': - ax2.xaxis.set_minor_locator(AutoMinorLocator(5)) + if filters is not None: + if scale[0] == 'linear': + ax2.xaxis.set_minor_locator(AutoMinorLocator(5)) - if residuals is not None and scale[0] == 'linear': - ax3.xaxis.set_minor_locator(AutoMinorLocator(5)) + if residuals is not None: + if scale[0] == 'linear': + ax3.xaxis.set_minor_locator(AutoMinorLocator(5)) if residuals is not None and filters is not None: ax1.set_xlabel('') @@ -755,6 +757,15 @@ def plot_spectrum(boxes: list, else: ax1.legend(**legend) + if scale[0] == 'log': + ax1.xaxis.set_major_formatter(ScalarFormatter()) + + if ax2 is not None: + ax2.xaxis.set_major_formatter(ScalarFormatter()) + + if ax3 is not None: + ax3.xaxis.set_major_formatter(ScalarFormatter()) + # filters = ['Paranal/SPHERE.ZIMPOL_N_Ha', # 'MUSE/Hbeta', # 'ALMA/855'] diff --git a/species/read/read_model.py b/species/read/read_model.py index 2b25ff1e..4e0da962 100644 --- a/species/read/read_model.py +++ b/species/read/read_model.py @@ -692,25 +692,25 @@ def get_model(self, if 'lognorm_radius' in model_param and 'lognorm_sigma' in model_param and \ 'lognorm_ext' in model_param: - model_box.flux = self.apply_lognorm_ext(model_box.wavelength, - model_box.flux, - model_param['lognorm_radius'], - model_param['lognorm_sigma'], - model_param['lognorm_ext']) + model_box.flux = self.apply_ext_norm_distr(model_box.wavelength, + model_box.flux, + model_param['lognorm_radius'], + model_param['lognorm_sigma'], + model_param['lognorm_ext']) if 'powerlaw_max' in model_param and 'powerlaw_exp' in model_param and \ 'powerlaw_ext' in model_param: - model_box.flux = self.apply_powerlaw_ext(model_box.wavelength, - model_box.flux, - model_param['powerlaw_max'], - model_param['powerlaw_exp'], - model_param['powerlaw_ext']) + model_box.flux = self.apply_ext_plaw_distr(model_box.wavelength, + model_box.flux, + model_param['powerlaw_max'], + model_param['powerlaw_exp'], + model_param['powerlaw_ext']) if 'ism_ext' in model_param: ism_reddening = model_param.get('ism_red', 3.1) - model_box.flux = self.apply_ism_ext(model_box.wavelength, + model_box.flux = self.apply_ext_ism(model_box.wavelength, model_box.flux, model_param['ism_ext'], ism_reddening) From e7c5052aa4b1b9c87eb765e4b71e294c2cf7234a Mon Sep 17 00:00:00 2001 From: Tomas Stolker Date: Fri, 7 May 2021 09:08:49 +0200 Subject: [PATCH 2/7] Changed out_file parameter in Database.get_samples to json_file --- species/data/database.py | 26 +++++++++++--------------- 1 file changed, 11 insertions(+), 15 deletions(-) diff --git a/species/data/database.py b/species/data/database.py index 53e9c0ae..1c44ce8d 100644 --- a/species/data/database.py +++ b/species/data/database.py @@ -4,6 +4,7 @@ import configparser import os +import json import warnings from typing import Dict, List, Optional, Tuple, Union @@ -1828,7 +1829,7 @@ def get_samples(self, tag: str, burnin: Optional[int] = None, random: Optional[int] = None, - out_file: Optional[str] = None) -> box.SamplesBox: + json_file: Optional[str] = None) -> box.SamplesBox: """ Parameters ---------- @@ -1841,10 +1842,9 @@ def get_samples(self, random : int, None Number of random samples to select. All samples (with the burnin excluded) are selected if set to ``None``. - out_file : str, None - Output file to store the posterior samples. The data will be stored in a FITS file if - the argument of ``out_file`` ends with `.fits`. Otherwise, the data will be written to - a text file. The data will not be written to a file if the argument is set to ``None``. + json_file : str, None + JSON file to store the posterior samples. The data will not be written if the argument + is set to ``None``. Returns ------- @@ -1899,18 +1899,14 @@ def get_samples(self, median_sample = self.get_median_sample(tag, burnin) - if out_file is not None: - header = '' - for i, item in enumerate(param): - header += f'{item}' - if i != len(param) - 1: - header += ' - ' + if json_file is not None: + samples_dict = {} - if out_file.endswith('.fits'): - fits.writeto(out_file, samples, overwrite=True) + for i, item in enumerate(param): + samples_dict[item] = list(samples[:, i]) - else: - np.savetxt(out_file, samples, header=header) + with open(json_file, 'w') as out_file: + json.dump(samples_dict, out_file, indent=4) return box.create_box('samples', spectrum=spectrum, From 6157890f7bb0c50539a9f9b471b690785776912c Mon Sep 17 00:00:00 2001 From: Tomas Stolker Date: Mon, 10 May 2021 16:17:02 +0200 Subject: [PATCH 3/7] Bug fix with changing name of apply_ism_ext to apply_ext_ism --- species/read/read_model.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/species/read/read_model.py b/species/read/read_model.py index 4e0da962..aad2964e 100644 --- a/species/read/read_model.py +++ b/species/read/read_model.py @@ -509,7 +509,7 @@ def apply_powerlaw_ext(wavelength: np.ndarray, @staticmethod @typechecked - def apply_ism_ext(wavelengths: np.ndarray, + def apply_ext_ism(wavelengths: np.ndarray, flux: np.ndarray, v_band_ext: float, v_band_red: float) -> np.ndarray: @@ -963,7 +963,7 @@ def get_data(self, if 'ism_ext' in model_param: ism_reddening = model_param.get('ism_red', 3.1) - model_box.flux = self.apply_ism_ext(model_box.wavelength, + model_box.flux = self.apply_ext_ism(model_box.wavelength, model_box.flux, model_param['ism_ext'], ism_reddening) From 727983c0e2cc189b6940f7ee13fe9d42ab329502 Mon Sep 17 00:00:00 2001 From: Tomas Stolker Date: Sat, 15 May 2021 16:54:11 +0200 Subject: [PATCH 4/7] Added adapt_logg parameter to ReadIsochrone.get_color_magnitude, storing the masses and raddi in the Box from ReadIsochrone, added read_util.get_radius function for converting a log(g) and mass into a radius, updated companion data, bug fix and additional output in Database.add_comparison, support for the Keck/NIRC2.NB_4.05 filter, pylint improvements --- species/analysis/compare_spectra.py | 2 +- species/core/box.py | 22 +++++-- species/core/constants.py | 1 + species/data/companions.py | 80 ++++++++++++++++++++------ species/data/database.py | 10 +++- species/data/filters.py | 14 +++++ species/plot/plot_color.py | 39 ++++++++----- species/plot/plot_comparison.py | 12 ++-- species/plot/plot_spectrum.py | 19 ++++-- species/read/read_isochrone.py | 89 +++++++++++++++++++++-------- species/util/phot_util.py | 4 +- species/util/read_util.py | 33 +++++++++-- 12 files changed, 245 insertions(+), 80 deletions(-) diff --git a/species/analysis/compare_spectra.py b/species/analysis/compare_spectra.py index dcf8a647..69863364 100644 --- a/species/analysis/compare_spectra.py +++ b/species/analysis/compare_spectra.py @@ -310,7 +310,7 @@ def compare_model(self, """ Method for finding the best fitting spectrum from a grid of atmospheric model spectra by evaluating the goodness-of-fit statistic from Cushing et al. (2008). Currently, this method - only supports model grids with only :math:`T_\mathrm{eff}` and :math:`\log(g)` as free + only supports model grids with only :math:`T_\\mathrm{eff}` and :math:`\\log(g)` as free parameters (e.g. BT-Settl). Please create an issue on Github if support for models with more than two parameters is required. diff --git a/species/core/box.py b/species/core/box.py index 3b917054..8107197d 100644 --- a/species/core/box.py +++ b/species/core/box.py @@ -27,8 +27,13 @@ def create_box(boxtype, box.filter_mag = kwargs['filter_mag'] box.color = kwargs['color'] box.magnitude = kwargs['magnitude'] - box.sptype = kwargs['sptype'] box.names = kwargs['names'] + if 'sptype' in kwargs: + box.sptype = kwargs['sptype'] + if 'mass' in kwargs: + box.mass = kwargs['mass'] + if 'radius' in kwargs: + box.radius = kwargs['radius'] if boxtype == 'colorcolor': box = ColorColorBox() @@ -37,8 +42,13 @@ def create_box(boxtype, box.filters = kwargs['filters'] box.color1 = kwargs['color1'] box.color2 = kwargs['color2'] - box.sptype = kwargs['sptype'] box.names = kwargs['names'] + if 'sptype' in kwargs: + box.sptype = kwargs['sptype'] + if 'mass' in kwargs: + box.mass = kwargs['mass'] + if 'radius' in kwargs: + box.radius = kwargs['radius'] elif boxtype == 'isochrone': box = IsochroneBox() @@ -171,8 +181,10 @@ def __init__(self): self.filter_mag = None self.color = None self.magnitude = None - self.sptype = None self.names = None + self.sptype = None + self.mass = None + self.radius = None class ColorColorBox(Box): @@ -193,8 +205,10 @@ def __init__(self): self.filters = None self.color1 = None self.color2 = None - self.sptype = None self.names = None + self.sptype = None + self.mass = None + self.radius = None class IsochroneBox(Box): diff --git a/species/core/constants.py b/species/core/constants.py index 92241aaa..4f55a2bc 100644 --- a/species/core/constants.py +++ b/species/core/constants.py @@ -12,3 +12,4 @@ M_EARTH = 5.9722e24 # (kg) R_EARTH = 6.3781e6 # (m) SIGMA_SB = 5.670374419e-8 # (W m−2 K−4) +ATOMIC_MASS = 1.66053906660e-27 # (kg) diff --git a/species/data/companions.py b/species/data/companions.py index 45e84788..a011b4db 100644 --- a/species/data/companions.py +++ b/species/data/companions.py @@ -53,10 +53,10 @@ def get_data() -> Dict[str, Dict[str, Union[bool, Tuple[float, float], 'Paranal/NACO.Lp': (15.33, 0.12), # Stolker et al. 2020 'Paranal/NACO.NB405': (15.23, 0.22), # Stolker et al. 2020 'Paranal/NACO.Mp': (14.65, 0.29)}, # Stolker et al. 2020 - 'semi_major': (110., 45.), # Cheetham et al. 2019 - 'mass_star': (1.96, 0.04), # Chauvin et al. 2017 - 'mass_companion': (9.9, 1.8), # Marleau et al. 2019 - 'accretion': False}, + 'semi_major': (110., 45.), # Cheetham et al. 2019 + 'mass_star': (1.96, 0.04), # Chauvin et al. 2017 + 'mass_companion': (9.9, 1.8), # Marleau et al. 2019 + 'accretion': False}, '51 Eri b': {'distance': (29.78, 0.12), 'app_mag': {'MKO/NSFCam.J': (19.04, 0.40), # Rajan et al. 2017 @@ -158,7 +158,9 @@ def get_data() -> Dict[str, Dict[str, Union[bool, Tuple[float, float], 'semi_major': (20.8, 0.7), # Wang et al. 2021 'mass_star': (0.98, 0.07), # Wang et al. 2021 'mass_companion': (3.2, 1.6), # Wang et al. 2021 - 'accretion': True}, # Haffert et al. 2019 + 'accretion': True, # Haffert et al. 2019 + 'line_flux': {'h-alpha': (8.1e-19, 0.3e-19), # Hashimoto et al. 2020 + 'h-beta': (2.3e-19, 2.3e-19)}}, # Hashimoto et al. 2020 'PDS 70 c': {'distance': (113.43, 0.52), 'app_mag': {'Paranal/NACO.NB405': (14.91, 0.35), # Stolker et al. 2020 @@ -179,7 +181,7 @@ def get_data() -> Dict[str, Dict[str, Union[bool, Tuple[float, float], 'Paranal/NACO.Lp': (15.28, 0.14)}, # Chauvin et al. 2004 'semi_major': (46., 46.), # Patience et al. 2010 'mass_star': (25.*constants.M_JUP/constants.M_SUN, - 5.*constants.M_JUP/constants.M_SUN), # Mohanty et al. 2007 + 5.*constants.M_JUP/constants.M_SUN), # Mohanty et al. 2007 'mass_companion': (8., 2.), # Mohanty et al. 2007 'accretion': False}, @@ -223,6 +225,9 @@ def get_data() -> Dict[str, Dict[str, Union[bool, Tuple[float, float], 'Magellan/VisAO.ip': (18.89, 0.24), # Wu et al. 2017 'Magellan/VisAO.zp': (16.40, 0.10), # Wu et al. 2017 'Magellan/VisAO.Ys': (15.88, 0.10), # Wu et al. 2017 + 'MKO/NSFCam.H': (14.02, 0.13), # Stolker et al. in prep. + 'Paranal/NACO.NB405': (12.29, 0.07), # Stolker et al. in prep. + 'Paranal/NACO.Mp': (11.97, 0.08), # Stolker et al. in prep. 'Paranal/NACO.Ks': [(13.474, 0.031), # Ginski et al. 2014 (13.386, 0.032), # Ginski et al. 2014 (13.496, 0.050), # Ginski et al. 2014 @@ -233,7 +238,10 @@ def get_data() -> Dict[str, Dict[str, Union[bool, Tuple[float, float], 'semi_major': (150., 50.), # Schwarz et al. 2016 'mass_star': (1.03, 0.05), # MacGregor et al. 2017 'mass_companion': (25., 15.), # Wu et al. 2017 - 'accretion': True}, # Seifahrt et al. 2007 + 'radius_companion': (3.6, 0.1), # Stolker et al. in prep. + 'accretion': True, # Seifahrt et al. 2007 + 'line_flux': {'h-alpha': (3.31e-18, 0.04e-18), # Stolker et al. in prep. + 'pa-beta': (1.32e-18, 0.01e-18)}}, # Stolker et al. in prep. 'PZ Tel B': {'distance': (47.13, 0.13), 'app_mag': {'Paranal/SPHERE.ZIMPOL_R_PRIM': (17.84, 0.31), # Maire et al. 2015 @@ -250,17 +258,18 @@ def get_data() -> Dict[str, Dict[str, Union[bool, Tuple[float, float], 'Paranal/NACO.Mp': (10.93, 0.03), # Stolker et al. 2020 'Gemini/NICI.ED286': (11.68, 0.14), # Biller et al. 2010 'Gemini/NIRI.H2S1v2-1-G0220': (11.39, 0.14)}, # Biller et al. 2010 - 'semi_major': (25., 25.), # Maire et al. 2016 - 'mass_star': (1.2, 1.2), # Ginski et al. 2014 - 'mass_companion': (55., 17.), # Maire et al. 2016 - 'accretion': False}, + 'semi_major': (25., 25.), # Maire et al. 2016 + 'mass_star': (1.2, 1.2), # Ginski et al. 2014 + 'mass_companion': (55., 17.), # Maire et al. 2016 + 'accretion': False, + 'line_flux': {'h-alpha': (2.2e-18, 0.9e-18)}}, # Musso Barcucci et al. 2019 'kappa And b': {'distance': (50.06, 0.87), 'app_mag': {'Subaru/CIAO.J': (15.86, 0.21), # Bonnefoy et al. 2014 'Subaru/CIAO.H': (14.95, 0.13), # Bonnefoy et al. 2014 'Subaru/CIAO.Ks': (14.32, 0.09), # Bonnefoy et al. 2014 'Keck/NIRC2.Lp': (13.12, 0.1), # Bonnefoy et al. 2014 - # 'Keck/NIRC2.NB_4.05': (13.0, 0.2), # Bonnefoy et al. 2014 + 'Keck/NIRC2.NB_4.05': (13.0, 0.2), # Bonnefoy et al. 2014 'LBT/LMIRCam.M_77K': (13.3, 0.3)}, # Bonnefoy et al. 2014 'semi_major': (55., 55.), # Bonnefoy et al. 2014 'mass_star': (2.7, 0.1), # Jones et al. 2016 @@ -292,7 +301,7 @@ def get_data() -> Dict[str, Dict[str, Union[bool, Tuple[float, float], 'Keck/NIRC2.H': (15.88, 0.05), # Daemgen et al. 2017 'Keck/NIRC2.Ks': (15.01, 0.06), # Daemgen et al. 2017 'Keck/NIRC2.Lp': (13.97, 0.06), # Daemgen et al. 2017 - # 'Keck/NIRC2.NB_4.05': (13.90, 0.08), # Daemgen et al. 2017 + 'Keck/NIRC2.NB_4.05': (13.90, 0.08), # Daemgen et al. 2017 'Keck/NIRC2.Ms': (14.01, 0.23)}, # Daemgen et al. 2017 'semi_major': (157., 157.), # Currie et al. 2014 'mass_star': (0.89, 0.08), # Kraus et al. 2014 @@ -337,7 +346,9 @@ def get_data() -> Dict[str, Dict[str, Union[bool, Tuple[float, float], 'semi_major': (84., 84.), # Delorme et al. 2013 'mass_star': (0.19, 0.02), # Delorme et al. 2013 'mass_companion': (13., 1.), # Delorme et al. 2013 - 'accretion': True}, # Eriksson et al. 2020 + 'accretion': True, # Eriksson et al. 2020 + 'line_flux': {'h-alpha': (12.80e-19, 0.70e-19), + 'h-beta': (1.39e-19, 0.10e-19)}}, # Eriksson et al. 2020 '1RXS 1609 B': {'distance': (139.67, 1.33), 'app_mag': {'Gemini/NIRI.J-G0202w': (17.90, 0.12), # Lafreniere et al. 2008 @@ -358,7 +369,9 @@ def get_data() -> Dict[str, Dict[str, Union[bool, Tuple[float, float], 'semi_major': (320., 320.), # Bowler et al. 2011 'mass_star': (0.9, 0.1), # Bowler et al. 2011 'mass_companion': (14., 2.), # Bowler et al. 2011 - 'accretion': True}, # Bowler et al. 2011 + 'accretion': True, # Bowler et al. 2011 + 'line_flux': {'h-alpha': (7.08e-19, 2.12e-18), # Zhou et al. 2014 + 'pa-beta': (1.12e-18, 0.03e-18)}}, # Bowler et al. 2011 'HD 72946 B': {'distance': (25.87, 0.03), 'app_mag': {'Paranal/SPHERE.IRDIS_D_H23_2': (14.56, 0.07), # Maire et al. 2019 @@ -434,7 +447,9 @@ def get_data() -> Dict[str, Dict[str, Union[bool, Tuple[float, float], 'mass_star': (2.0, 0.3), # Mendigutía et al. 2014 'mass_companion': (0.13*constants.M_SUN/constants.M_JUP, 0.03*constants.M_SUN/constants.M_JUP), # Lacour et al. 2016 - 'accretion': True}, # Close et al. 2014 + 'radius_companion': (19.1, 1.0), # Christiaens et al. 2018 + 'accretion': True, # Close et al. 2014 + 'line_flux': {'h-alpha': (7.6e-17, 3.5e-17)}}, # Cugno et al. 2019 TODO extinction? 'CS Cha B': {'distance': (168.77, 1.92), 'app_mag': {'Paranal/SPHERE.IRDIS_B_J': (19.16, 0.21), # Ginski et al. 2018 @@ -444,6 +459,37 @@ def get_data() -> Dict[str, Dict[str, Union[bool, Tuple[float, float], 'mass_star': (1.0, 0.1), # Ginski et al. 2018 'mass_companion': (0.3*constants.M_SUN/constants.M_JUP, 0.1*constants.M_SUN/constants.M_JUP), # Haffert et al. 2020 - 'accretion': True}} # Haffert et al. 2020 + 'accretion': True, # Haffert et al. 2020 + 'line_flux': {'h-alpha': (17.3e-20, 2.1e-20)}}, # Haffert et al. 2020 + + 'CT Cha B': {'distance': (189.95, 0.42), + 'app_mag': {'Paranal/NACO.J': (16.61, 0.30), # Schmidt et al. 2008 + 'Paranal/NACO.Ks': [(14.95, 0.30), # Schmidt et al. 2008 + (14.89, 0.30)]}, # Schmidt et al. 2008 + 'semi_major': (430., 0.), # Wu et al. 2015 + 'mass_star': (0.55, 0.), # Hartmann et al. 1998 + 'mass_companion': (19., 5.), # Wu et al. 2015 + 'accretion': True}, # Wu et al. 2015 + + 'SR 12 C': {'distance': (125., 25.), # Bouvier & Appenzeller 1992 + 'app_mag': {'MKO/NSFCam.J': (15.93, 0.03), # Kuzuhara et al. 2011 + 'MKO/NSFCam.H': (15.18, 0.03), # Kuzuhara et al. 2011 + 'MKO/NSFCam.Ks': (14.57, 0.03), # Kuzuhara et al. 2011 + 'MKO/NSFCam.Lp': (13.10, 0.08)}, # Kuzuhara et al. 2011 + 'semi_major': (1100., 0.), # Bowler et al. 2014 + 'mass_star': (1.05, 0.05), # Bowler et al. 2014 + 'mass_companion': (13., 7.), # Kuzuhara et al. 2011 + 'accretion': True, # Santamaría-Miranda et al. 2017 + 'line_flux': {'h-alpha': (1.34e-18, 0.05e-18), # Santamaría-Miranda et al. 2017 (erratum) + 'h-beta': (2.19e-19, 0.03e-19)}}, # Santamaría-Miranda et al. 2017 (erratum) + + 'DH Tau B': {'distance': (133.45, 0.45), + 'app_mag': {'Subaru/CIAO.J': (15.71, 0.05), # Itoh et al. 2005 + 'Subaru/CIAO.H': (14.96, 0.04), # Itoh et al. 2005 + 'Subaru/CIAO.Ks': (14.19, 0.02)}, # Itoh et al. 2005 + 'semi_major': (330., 0.), # Patience et al. 2012 + 'mass_star': (0.33, 0.), # Patience et al. 2012 + 'mass_companion': (11., 3.), # Patience et al. 2012 + 'accretion': True}} # Zhou et al. 2014 return data diff --git a/species/data/database.py b/species/data/database.py index 1c44ce8d..f3a9ffda 100644 --- a/species/data/database.py +++ b/species/data/database.py @@ -1432,7 +1432,8 @@ def get_compare_sample(self, model_param['distance'] = dset.attrs['distance'] if n_spec_name == 1: - model_param['radius'] = dset.attrs[f'radius_{item}'] + spec_name = dset.attrs['spec_name0'] + model_param['radius'] = dset.attrs[f'radius_{spec_name}'] else: if spec_fix is None: @@ -2070,12 +2071,15 @@ def add_comparison(self, # Indices of the best-fit model best_index = np.unravel_index(goodness_sum.argmin(), goodness_sum.shape) + dset.attrs['best_fit'] = goodness_sum[best_index] print('Best-fit parameters:') + print(f' - Goodness-of-fit = {goodness_sum[best_index]:.2e}') for i, item in enumerate(model_param): best_param = coord_points[i][best_index[i]] dset.attrs[f'best_param{i}'] = best_param + print(f' - {item} = {best_param}') for i, item in enumerate(spec_name): scaling = flux_scaling[best_index[0], best_index[1], best_index[2], i] @@ -2084,7 +2088,7 @@ def add_comparison(self, radius /= constants.R_JUP # (Rjup) dset.attrs[f'radius_{item}'] = radius - print(f' - {item} radius (Rjup) = {radius:.2f}') + print(f' - {item} radius (Rjup) = {radius:.2f}') dset.attrs[f'scaling_{item}'] = scaling - print(f' - {item} scaling = {scaling:.2e}') + print(f' - {item} scaling = {scaling:.2e}') diff --git a/species/data/filters.py b/species/data/filters.py index f0571678..2afeb2a5 100644 --- a/species/data/filters.py +++ b/species/data/filters.py @@ -70,6 +70,20 @@ def download_filter(filter_id: str) -> Tuple[Optional[np.ndarray], os.remove('VisAO_zp_filter_curve.dat') + elif filter_id == 'Keck/NIRC2.NB_4.05': + # The filter profile of Br_alpha has been digitized from + # https://www2.keck.hawaii.edu/inst/nirc2/filters.html + + url = 'https://home.strw.leidenuniv.nl/~stolker/species/Keck_NIRC2.NB_4.05.dat' + urllib.request.urlretrieve(url, 'Keck_NIRC2.NB_4.05.dat') + + wavelength, transmission = np.loadtxt('Keck_NIRC2.NB_4.05.dat', unpack=True) + + # Not sure if energy- or photon-counting detector + det_type = 'photon' + + os.remove('Keck_NIRC2.NB_4.05.dat') + elif filter_id in ['LCO/VisAO.Ys', 'Magellan/VisAO.Ys']: url = 'https://xwcl.science/magao/visao/VisAO_Ys_filter_curve.dat' urllib.request.urlretrieve(url, 'VisAO_Ys_filter_curve.dat') diff --git a/species/plot/plot_color.py b/species/plot/plot_color.py index 2c7a16b4..517ad62e 100644 --- a/species/plot/plot_color.py +++ b/species/plot/plot_color.py @@ -75,7 +75,8 @@ def plot_color_magnitude(boxes: list, Plot accreting, directly imaged objects with a different symbol than the regular, directly imaged objects. The object names from ``objects`` will be compared with the data from :func:`~species.data.companions.get_data` to check if a companion is accreting or not. - reddening : list(tuple(tuple(str, str), tuple(str, float), str, float, tuple(float, float)), None + reddening : list(tuple(tuple(str, str), tuple(str, float), str, float, + tuple(float, float)), None Include reddening arrows by providing a list with tuples. Each tuple contains the filter names for the color, the filter name and value of the magnitude, the particle radius (um), and the start position (color, mag) of the arrow in the plot, so ``((filter_color_1, @@ -211,7 +212,8 @@ def plot_color_magnitude(boxes: list, label = plot_util.model_name(item.library) if item.library == 'zhu2015': - ax1.plot(item.color, item.magnitude, marker='x', ms=5, linestyle=model_linestyle[model_count[1]], + ax1.plot(item.color, item.magnitude, marker='x', ms=5, + linestyle=model_linestyle[model_count[1]], linewidth=0.6, color='gray', label=label, zorder=0) xlim = ax1.get_xlim() @@ -223,11 +225,12 @@ def plot_color_magnitude(boxes: list, if item.magnitude[i] > ylim[1]: ax1.annotate(teff_label, (item.color[i], item.magnitude[i]), color='gray', fontsize=8, ha='left', va='center', - xytext=(item.color[i]+0.1, item.magnitude[i]+0.05), zorder=3) + xytext=(item.color[i]+0.1, item.magnitude[i]+0.05), + zorder=3) else: ax1.plot(item.color, item.magnitude, linestyle=model_linestyle[model_count[1]], - linewidth=1., color=model_color[model_count[0]], label=label, zorder=0) + lw=1., color=model_color[model_count[0]], label=label, zorder=0) if mass_labels is not None: interp_magnitude = interp1d(item.sptype, item.magnitude) @@ -262,12 +265,13 @@ def plot_color_magnitude(boxes: list, if xlim[0]+0.2 < pos_color < xlim[1]-0.2 and \ ylim[1]+0.2 < pos_mag < ylim[0]-0.2: - ax1.scatter(pos_color, pos_mag, c=model_color[model_count[0]], s=15, - edgecolor='none', zorder=0) + ax1.scatter(pos_color, pos_mag, c=model_color[model_count[0]], + s=15, edgecolor='none', zorder=0) ax1.annotate(mass_label, (pos_color, pos_mag), color=model_color[model_count[0]], fontsize=9, - xytext=mass_xytext, zorder=3, ha=mass_ha, va='center') + xytext=mass_xytext, zorder=3, ha=mass_ha, + va='center') else: ax1.plot(item.color, item.magnitude, linestyle=model_linestyle[model_count[1]], @@ -515,8 +519,8 @@ def plot_color_magnitude(boxes: list, @typechecked def plot_color_color(boxes: list, objects: Optional[Union[List[Tuple[str, Tuple[str, str], Tuple[str, str]]], - List[Tuple[str, Tuple[str, str], Tuple[str, str], Optional[dict], - Optional[dict]]]]] = None, + List[Tuple[str, Tuple[str, str], Tuple[str, str], + Optional[dict], Optional[dict]]]]] = None, mass_labels: Optional[Union[List[float], List[Tuple[float, str]]]] = None, teff_labels: Optional[Union[List[float], List[Tuple[float, str]]]] = None, companion_labels: bool = False, @@ -562,7 +566,8 @@ def plot_color_color(boxes: list, to None. companion_labels : bool Plot labels with the names of the directly imaged companions. - reddening : list(tuple(tuple(str, str), tuple(str, str), tuple(str, float), str, float, tuple(float, float)), None + reddening : list(tuple(tuple(str, str), tuple(str, str), tuple(str, float), str, float, + tuple(float, float)), None Include reddening arrows by providing a list with tuples. Each tuple contains the filter names for the color, the filter name for the magnitude, the particle radius (um), and the start position (color, mag) of the arrow in the plot, so (filter_color_1, filter_color_2, @@ -685,7 +690,8 @@ def plot_color_color(boxes: list, label = plot_util.model_name(item.library) if item.library == 'zhu2015': - ax1.plot(item.color1, item.color2, marker='x', ms=5, linestyle=model_linestyle[model_count[1]], + ax1.plot(item.color1, item.color2, marker='x', ms=5, + linestyle=model_linestyle[model_count[1]], linewidth=0.6, color='gray', label=label, zorder=0) xlim = ax1.get_xlim() @@ -697,11 +703,12 @@ def plot_color_color(boxes: list, if item.color2[i] < ylim[1]: ax1.annotate(teff_label, (item.color1[i], item.color2[i]), color='gray', fontsize=8, ha='left', va='center', - xytext=(item.color1[i]+0.1, item.color2[i]-0.05), zorder=3) + xytext=(item.color1[i]+0.1, item.color2[i]-0.05), + zorder=3) else: ax1.plot(item.color1, item.color2, linestyle=model_linestyle[model_count[1]], - linewidth=1., color=model_color[model_count[0]], label=label, zorder=0) + lw=1., color=model_color[model_count[0]], label=label, zorder=0) if mass_labels is not None: interp_color1 = interp1d(item.sptype, item.color1) @@ -737,12 +744,14 @@ def plot_color_color(boxes: list, if xlim[0]+0.2 < pos_color1 < xlim[1]-0.2 and \ ylim[0]+0.2 < pos_color2 < ylim[1]-0.2: - ax1.scatter(pos_color1, pos_color2, c=model_color[model_count[0]], + ax1.scatter(pos_color1, pos_color2, + c=model_color[model_count[0]], s=15, edgecolor='none', zorder=0) ax1.annotate(mass_label, (pos_color1, pos_color2), color=model_color[model_count[0]], fontsize=9, - xytext=mass_xytext, ha=mass_ha, va='center', zorder=3) + xytext=mass_xytext, ha=mass_ha, va='center', + zorder=3) else: ax1.plot(item.color1, item.color2, linestyle=model_linestyle[model_count[1]], diff --git a/species/plot/plot_comparison.py b/species/plot/plot_comparison.py index 8b3f48a5..80b2725e 100644 --- a/species/plot/plot_comparison.py +++ b/species/plot/plot_comparison.py @@ -18,7 +18,7 @@ from species.core import constants from species.read import read_object -from species.util import dust_util, read_util, plot_util +from species.util import dust_util, read_util @typechecked @@ -353,8 +353,8 @@ def plot_grid_statistic(tag: str, read_obj = read_object.ReadObject(dset.attrs['object_name']) n_wavel = 0 - for key, value in read_obj.get_spectrum().items(): - n_wavel += value[0].shape[0] + for item in read_obj.get_spectrum().values(): + n_wavel += item[0].shape[0] goodness_fit = np.array(dset) @@ -436,14 +436,16 @@ def plot_grid_statistic(tag: str, x_grid, y_grid = np.meshgrid(x_new, y_new) goodness_fit = fit_interp((y_grid, x_grid)) + goodness_fit = np.log10(goodness_fit) + goodness_fit -= np.amin(goodness_fit) - c = ax.contourf(x_grid, y_grid, np.log10(goodness_fit)) + c = ax.contourf(x_grid, y_grid, goodness_fit) cb = mpl.colorbar.Colorbar(ax=ax_cb, mappable=c, orientation='vertical', ticklocation='right', format='%.1f') cb.ax.tick_params(width=0.8, length=5, labelsize=12, direction='in', color='black') - cb.ax.set_ylabel(r'$\mathregular{log}\,G_k$', rotation=270, labelpad=22, fontsize=13.) + cb.ax.set_ylabel(r'$\Delta\mathregular{log}\,G_k$', rotation=270, labelpad=22, fontsize=13.) if len(coord_points[2]) != 1: extra_interp = RegularGridInterpolator((coord_points[1], coord_points[0]), extra_map) diff --git a/species/plot/plot_spectrum.py b/species/plot/plot_spectrum.py index 9de03171..2e864410 100644 --- a/species/plot/plot_spectrum.py +++ b/species/plot/plot_spectrum.py @@ -89,7 +89,7 @@ def plot_spectrum(boxes: list, Alternatively, a list with two values can be provided to separate the model and data handles in two legends. Each of these two elements can be set to ``None``. For example, ``[None, {'loc': 'upper left', 'fontsize: 12.}]``, if only the data points should be - included in a legend. + included in a legend. figsize : tuple(float, float) Figure size. object_type : str @@ -551,15 +551,24 @@ def plot_spectrum(boxes: list, if quantity == 'flux': flux_scaling = wavelength + scale_tmp = flux_scaling / scaling + if isinstance(boxitem.flux[item][0], np.ndarray): for i in range(boxitem.flux[item].shape[1]): - plot_obj = ax1.errorbar(wavelength, flux_scaling*boxitem.flux[item][0, i]/scaling, xerr=fwhm/2., - yerr=flux_scaling*boxitem.flux[item][1, i]/scaling, marker='s', ms=5, zorder=3, color='black') + plot_obj = ax1.errorbar(wavelength, + scale_tmp*boxitem.flux[item][0, i], + xerr=fwhm/2., + yerr=scale_tmp*boxitem.flux[item][1, i], + marker='s', ms=5, zorder=3, color='black') else: - plot_obj = ax1.errorbar(wavelength, flux_scaling*boxitem.flux[item][0]/scaling, xerr=fwhm/2., - yerr=flux_scaling*boxitem.flux[item][1]/scaling, marker='s', ms=5, zorder=3, color='black') + + plot_obj = ax1.errorbar(wavelength, + scale_tmp*boxitem.flux[item][0], + xerr=fwhm/2., + yerr=scale_tmp*boxitem.flux[item][1], + marker='s', ms=5, zorder=3, color='black') plot_kwargs[j][item] = {'marker': 's', 'ms': 5., 'color': plot_obj[0].get_color()} diff --git a/species/read/read_isochrone.py b/species/read/read_isochrone.py index 598c24bc..001ebf12 100644 --- a/species/read/read_isochrone.py +++ b/species/read/read_isochrone.py @@ -11,11 +11,12 @@ import h5py import numpy as np -from typeguard import typechecked from scipy.interpolate import griddata +from typeguard import typechecked from species.core import box from species.read import read_model +from species.util import read_util class ReadIsochrone: @@ -161,7 +162,8 @@ def get_color_magnitude(self, masses: np.ndarray, model: str, filters_color: Tuple[str, str], - filter_mag: str) -> box.ColorMagBox: + filter_mag: str, + adapt_logg: bool = False) -> box.ColorMagBox: """ Function for calculating color-magnitude combinations from a selected isochrone. @@ -179,6 +181,11 @@ def get_color_magnitude(self, filter_mag : str Filter name for the absolute magnitude as listed in the file with the isochrone data. The value should be equal to one of the ``filters_color`` values. + adapt_logg : bool + Adapt :math:`\\log(g)` to the upper or lower boundary of the atmospheric model grid + whenever the :math:`\\log(g)` that has been calculated from the isochrone mass and + radius lies outside the available range of the synthetic spectra. Typically + :math:`\\log(g)` has only a minor impact on the broadband magnitudes and colors. Returns ------- @@ -194,6 +201,8 @@ def get_color_magnitude(self, model1 = read_model.ReadModel(model=model, filter_name=filters_color[0]) model2 = read_model.ReadModel(model=model, filter_name=filters_color[1]) + param_bounds = model1.get_bounds() + if model1.get_parameters() != ['teff', 'logg']: raise ValueError('Creating synthetic colors and magnitudes from isochrones is ' 'currently only implemented for models with only Teff and log(g) ' @@ -202,6 +211,7 @@ def get_color_magnitude(self, mag1 = np.zeros(isochrone.masses.shape[0]) mag2 = np.zeros(isochrone.masses.shape[0]) + radius = np.zeros(isochrone.masses.shape[0]) for i, mass_item in enumerate(isochrone.masses): model_param = {'teff': isochrone.teff[i], @@ -209,6 +219,8 @@ def get_color_magnitude(self, 'mass': mass_item, 'distance': 10.} + radius[i] = read_util.get_radius(model_param['logg'], model_param['mass']) # (Rjup) + if np.isnan(isochrone.teff[i]): mag1[i] = np.nan mag2[i] = np.nan @@ -217,26 +229,48 @@ def get_color_magnitude(self, f'{model_param}. Setting the magnitudes to NaN.') else: - for item_bounds in model1.get_bounds(): - if model_param[item_bounds] <= model1.get_bounds()[item_bounds][0]: - mag1[i] = np.nan - mag2[i] = np.nan - - warnings.warn(f'The value of {item_bounds} is {model_param[item_bounds]}, ' - f'which is below the lower bound of the model grid ' - f'({model1.get_bounds()[item_bounds][0]}). Setting the ' - f'magnitudes to NaN for the following isochrone sample: ' - f'{model_param}.') - - elif model_param[item_bounds] >= model1.get_bounds()[item_bounds][1]: - mag1[i] = np.nan - mag2[i] = np.nan - - warnings.warn(f'The value of {item_bounds} is {model_param[item_bounds]}, ' - f'which is above the upper bound of the model grid ' - f'({model1.get_bounds()[item_bounds][1]}). Setting the ' - f'magnitudes to NaN for the following isochrone sample: ' - f'{model_param}.') + for item_bounds in param_bounds: + if model_param[item_bounds] <= param_bounds[item_bounds][0]: + if adapt_logg and item_bounds == 'logg': + warnings.warn(f'The log(g) is {model_param[item_bounds]} but the ' + f'lower boundary of the model grid is ' + f'{param_bounds[item_bounds][0]}. Adapting ' + f'log(g) to {param_bounds[item_bounds][0]} since ' + f'adapt_logg=True.') + + model_param['logg'] = param_bounds['logg'][0] + + else: + mag1[i] = np.nan + mag2[i] = np.nan + + warnings.warn(f'The value of {item_bounds} is ' + f'{model_param[item_bounds]}, which is below the lower ' + f'bound of the model grid ' + f'({param_bounds[item_bounds][0]}). Setting the ' + f'magnitudes to NaN for the following isochrone sample: ' + f'{model_param}.') + + elif model_param[item_bounds] >= param_bounds[item_bounds][1]: + if adapt_logg and item_bounds == 'logg': + warnings.warn(f'The log(g) is {model_param[item_bounds]} but the ' + f'upper boundary of the model grid is ' + f'{param_bounds[item_bounds][1]}. Adapting ' + f'log(g) to {param_bounds[item_bounds][1]} since ' + f'adapt_logg=True.') + + model_param['logg'] = param_bounds['logg'][1] + + else: + mag1[i] = np.nan + mag2[i] = np.nan + + warnings.warn(f'The value of {item_bounds} is ' + f'{model_param[item_bounds]}, which is above the upper ' + f'bound of the model grid ' + f'({param_bounds[item_bounds][1]}). Setting the ' + f'magnitudes to NaN for the following isochrone sample: ' + f'{model_param}.') if not np.isnan(mag1[i]): mag1[i], _ = model1.get_magnitude(model_param) @@ -259,8 +293,10 @@ def get_color_magnitude(self, filter_mag=filter_mag, color=mag1-mag2, magnitude=abs_mag, + names=None, sptype=masses, - names=None) + mass=masses, + radius=radius) @typechecked def get_color_color(self, @@ -310,6 +346,7 @@ def get_color_color(self, mag2 = np.zeros(isochrone.masses.shape[0]) mag3 = np.zeros(isochrone.masses.shape[0]) mag4 = np.zeros(isochrone.masses.shape[0]) + radius = np.zeros(isochrone.masses.shape[0]) for i, mass_item in enumerate(isochrone.masses): model_param = {'teff': isochrone.teff[i], @@ -317,6 +354,8 @@ def get_color_color(self, 'mass': mass_item, 'distance': 10.} + radius[i] = read_util.get_radius(model_param['logg'], model_param['mass']) # (Rjup) + if np.isnan(isochrone.teff[i]): mag1[i] = np.nan mag2[i] = np.nan @@ -365,5 +404,7 @@ def get_color_color(self, filters=filters_colors, color1=mag1-mag2, color2=mag3-mag4, + names=None, sptype=masses, - names=None) + mass=masses, + radius=radius) diff --git a/species/util/phot_util.py b/species/util/phot_util.py index 1d4ccede..e710c3f0 100644 --- a/species/util/phot_util.py +++ b/species/util/phot_util.py @@ -256,9 +256,9 @@ def get_residuals(datatype: str, verbose=True) else: - readmodel = read_model.ReadModel(spectrum, wavel_range=wavel_range) + # Resampling to the new wavelength points is done by the get_model method - # resampling to the new wavelength points is done in teh get_model function + readmodel = read_model.ReadModel(spectrum, wavel_range=wavel_range) model_spec = readmodel.get_model(parameters, spec_res=spec_res, diff --git a/species/util/read_util.py b/species/util/read_util.py index a77908b4..3a481eb0 100644 --- a/species/util/read_util.py +++ b/species/util/read_util.py @@ -21,6 +21,8 @@ def get_mass(logg: Union[float, np.ndarray], radius: Union[float, np.ndarray]) -> Union[float, np.ndarray]: """ + Function for converting a :math:`\\log(g)` and a radius into a mass. + Parameters ---------- logg : float, np.ndarray @@ -35,13 +37,36 @@ def get_mass(logg: Union[float, np.ndarray], """ surface_grav = 1e-2 * 10.**logg # (m s-2) - radius *= constants.R_JUP # (m) - mass = surface_grav*radius**2/constants.GRAVITY # (kg) - mass /= constants.M_JUP # (Mjup) - return mass + return mass/constants.M_JUP + + +@typechecked +def get_radius(logg: Union[float, np.ndarray], + mass: Union[float, np.ndarray]) -> Union[float, np.ndarray]: + """ + Function for converting a :math:`\\log(g)` and a mass into a radius. + + Parameters + ---------- + logg : float, np.ndarray + Log10 of the surface gravity (cgs). + mass : float, np.ndarray + Mass (Mjup). + + Returns + ------- + float, np.ndarray + Radius (Rjup). + """ + + surface_grav = 1e-2 * 10.**logg # (m s-2) + mass *= constants.M_JUP # (kg) + radius = np.sqrt(mass*constants.GRAVITY/surface_grav) # (m) + + return radius/constants.R_JUP def add_luminosity(modelbox): From 492ca08214707cc5a61c22e347c00dac11c0684a Mon Sep 17 00:00:00 2001 From: Tomas Stolker Date: Wed, 26 May 2021 14:24:15 +0200 Subject: [PATCH 5/7] Added try-except with the import of UltraNest, support for spectra of companions in data.companions --- species/__init__.py | 1 + species/analysis/emission_line.py | 6 +- species/analysis/fit_model.py | 6 +- species/data/companions.py | 118 +++++++++++++++++++++++++++++- species/data/database.py | 20 +++-- species/plot/plot_comparison.py | 6 +- 6 files changed, 143 insertions(+), 14 deletions(-) diff --git a/species/__init__.py b/species/__init__.py index da7f5ebd..465b56f2 100644 --- a/species/__init__.py +++ b/species/__init__.py @@ -65,6 +65,7 @@ from species.util.read_util import add_luminosity, \ get_mass, \ + get_radius, \ powerlaw_spectrum, \ gaussian_spectrum, \ update_spectra diff --git a/species/analysis/emission_line.py b/species/analysis/emission_line.py index ae5bcd58..18e95b44 100644 --- a/species/analysis/emission_line.py +++ b/species/analysis/emission_line.py @@ -10,7 +10,11 @@ import matplotlib as mpl import matplotlib.pyplot as plt import numpy as np -import ultranest + +try: + import ultranest +except: + warnings.warn('UltraNest could not be imported.') from astropy import units as u from astropy.modeling.fitting import LinearLSQFitter diff --git a/species/analysis/fit_model.py b/species/analysis/fit_model.py index 865a0a8d..69d56de7 100644 --- a/species/analysis/fit_model.py +++ b/species/analysis/fit_model.py @@ -12,7 +12,11 @@ import emcee import numpy as np import spectres -import ultranest + +try: + import ultranest +except: + warnings.warn('UltraNest could not be imported.') # Installation of MultiNest is not possible on readthedocs try: diff --git a/species/data/companions.py b/species/data/companions.py index a011b4db..26894f61 100644 --- a/species/data/companions.py +++ b/species/data/companions.py @@ -2,7 +2,10 @@ Module for extracting data of directly imaged planets and brown dwarfs. """ -from typing import Dict, List, Tuple, Union +import os +import urllib.request + +from typing import Dict, List, Optional, Tuple, Union from typeguard import typechecked @@ -216,7 +219,8 @@ def get_data() -> Dict[str, Dict[str, Union[bool, Tuple[float, float], 0.01*constants.M_SUN/constants.M_JUP), # Kennedy et al. 2020 'accretion': False}, - 'GQ Lup B': {'distance': (151.82, 1.10), + 'GQ Lup B': {'distance': (154.10, 0.69), # Gaia Data Release 3 + 'parallax': (6.49, 0.03), # Gaia Data Release 3 'app_mag': {'HST/WFPC2-PC.F606W': (19.19, 0.07), # Marois et al. 2007 'HST/WFPC2-PC.F814W': (17.67, 0.05), # Marois et al. 2007 'HST/NICMOS2.F171M': (13.84, 0.13), # Marois et al. 2007 @@ -225,7 +229,7 @@ def get_data() -> Dict[str, Dict[str, Union[bool, Tuple[float, float], 'Magellan/VisAO.ip': (18.89, 0.24), # Wu et al. 2017 'Magellan/VisAO.zp': (16.40, 0.10), # Wu et al. 2017 'Magellan/VisAO.Ys': (15.88, 0.10), # Wu et al. 2017 - 'MKO/NSFCam.H': (14.02, 0.13), # Stolker et al. in prep. + # 'MKO/NSFCam.H': (14.02, 0.13), # Stolker et al. in prep. 'Paranal/NACO.NB405': (12.29, 0.07), # Stolker et al. in prep. 'Paranal/NACO.Mp': (11.97, 0.08), # Stolker et al. in prep. 'Paranal/NACO.Ks': [(13.474, 0.031), # Ginski et al. 2014 @@ -449,7 +453,7 @@ def get_data() -> Dict[str, Dict[str, Union[bool, Tuple[float, float], 0.03*constants.M_SUN/constants.M_JUP), # Lacour et al. 2016 'radius_companion': (19.1, 1.0), # Christiaens et al. 2018 'accretion': True, # Close et al. 2014 - 'line_flux': {'h-alpha': (7.6e-17, 3.5e-17)}}, # Cugno et al. 2019 TODO extinction? + 'line_flux': {'h-alpha': (7.6e-17, 3.5e-17)}}, # Cugno et al. 2019 'CS Cha B': {'distance': (168.77, 1.92), 'app_mag': {'Paranal/SPHERE.IRDIS_B_J': (19.16, 0.21), # Ginski et al. 2018 @@ -493,3 +497,109 @@ def get_data() -> Dict[str, Dict[str, Union[bool, Tuple[float, float], 'accretion': True}} # Zhou et al. 2014 return data + + +@typechecked +def get_spec_data() -> Dict[str, Dict[str, Tuple[str, Optional[str], float, str]]]: + """ + Function for extracting a dictionary with the spectra of directly imaged planets. + + Returns + ------- + dict + Dictionary with the spectrum, optional covariances, spectral resolution, and filename. + """ + + spec_data = {'beta Pic b': {'GPI_YJHK': ('betapicb_gpi_yjhk.dat', + None, + 40., + 'Chilcote et al. 2017, AJ, 153, 182'), + + 'GRAVITY': ('BetaPictorisb_2018-09-22.fits', + 'BetaPictorisb_2018-09-22.fits', + 500., + 'Gravity Collaboration et al. 2020, A&A, 633, 110')}, + + '51 Eri b': {'SPHERE_YJH': ('51erib_sphere_yjh.dat', + None, + 25., + 'Samland et al. 2017, A&A, 603, 57')}, + + 'HD 206893 B': {'SPHERE_YJH': ('hd206893b_sphere_yjh.dat', + None, + 25., + 'Delorme et al. 2017, A&A, 608, 79')}, + + 'HIP 65426 B': {'SPHERE_YJH': ('hip65426b_sphere_yjh.dat', + None, + 25., + 'Cheetham et al. 2019, A&A, 622, 80')}, + + 'HR 8799 e': {'SPHERE_YJH': ('hr8799e_sphere_yjh.dat', + None, + 25., + 'Zurlo et al. 2016, A&A, 587, 57')}, + + 'PDS 70 b': {'SPHERE_YJH': ('pds70b_sphere_yjh.dat', + None, + 25., + 'Müller et al. 2018, A&A, 617, 2')}} + + return spec_data + + +@typechecked +def companion_spectra(input_path: str, + comp_name: str) -> Optional[Dict[str, Tuple[str, Optional[str], float]]]: + """ + Function for getting available spectra of directly imaged planets and brown dwarfs. + + Parameters + ---------- + input_path : str + Path of the data folder. + comp_name : str + Companion name for which the spectra will be returned. + + Returns + ------- + dict, None + Dictionary with the spectra of ``comp_name``. A ``None`` will be returned if there are not + any spectra available. + """ + + data_folder = os.path.join(input_path, 'companion_data/') + + if not os.path.exists(data_folder): + os.makedirs(data_folder) + + spec_data = get_spec_data() + + if comp_name in spec_data: + spec_dict = {} + + for key, value in spec_data[comp_name].items(): + print(f'Getting {key} spectrum of {comp_name}...', end='', flush=True) + + spec_url = f'https://home.strw.leidenuniv.nl/~stolker/species/spectra/{value[0]}' + spec_file = os.path.join(data_folder, value[0]) + + if value[1] is None: + cov_file = None + else: + cov_file = os.path.join(data_folder, value[1]) + + if not os.path.isfile(spec_file): + urllib.request.urlretrieve(spec_url, spec_file) + + spec_dict[key] = (spec_file, cov_file, value[2]) + + print(' [DONE]') + + print(f'IMPORTANT: Please cite {value[3]}') + print(' when making use of this spectrum in a publication') + + else: + spec_dict = None + + return spec_dict diff --git a/species/data/database.py b/species/data/database.py index f3a9ffda..13bd30f6 100644 --- a/species/data/database.py +++ b/species/data/database.py @@ -103,6 +103,8 @@ def list_companions() -> None: None """ + spec_data = companions.get_spec_data() + for planet_name, planet_dict in companions.get_data().items(): distance = planet_dict['distance'] app_mag = planet_dict['app_mag'] @@ -113,6 +115,10 @@ def list_companions() -> None: for mag_name, mag_dict in app_mag.items(): print(f'{mag_name} (mag) = {mag_dict[0]} +/- {mag_dict[1]}') + if planet_name in spec_data: + for key, value in spec_data[planet_name].items(): + print(f'{key} spectrum from {value[3]}') + print() @typechecked @@ -145,8 +151,9 @@ def delete_data(self, def add_companion(self, name: Union[Optional[str], Optional[List[str]]] = None) -> None: """ - Function for adding the magnitudes of directly imaged planets and brown dwarfs from - :class:`~species.data.companions.get_data` to the database. + Function for adding the magnitudes and spectra of directly imaged planets and brown dwarfs + from :class:`~species.data.companions.get_data` and + :class:`~species.data.companions.get_comp_spec`to the database. Parameters ---------- @@ -170,9 +177,12 @@ def add_companion(self, name = data.keys() for item in name: + spec_dict = companions.companion_spectra(self.input_path, item) + self.add_object(object_name=item, distance=data[item]['distance'], - app_mag=data[item]['app_mag']) + app_mag=data[item]['app_mag'], + spectrum=spec_dict) @typechecked def add_dust(self) -> None: @@ -886,8 +896,8 @@ def add_object(self, print(f' - Filename: {value[0]}') print(f' - Data shape: {read_spec[key].shape}') print(f' - Wavelength range (um): {wavelength[0]:.2f} - {wavelength[-1]:.2f}') - print(f' - Mean flux (W m-2 um-1): {np.mean(flux):.2e}') - print(f' - Mean error (W m-2 um-1): {np.mean(error):.2e}') + print(f' - Mean flux (W m-2 um-1): {np.nanmean(flux):.2e}') + print(f' - Mean error (W m-2 um-1): {np.nanmean(error):.2e}') if isinstance(deredden, float): print(f' - Dereddening A_V: {deredden}') diff --git a/species/plot/plot_comparison.py b/species/plot/plot_comparison.py index 80b2725e..28af6468 100644 --- a/species/plot/plot_comparison.py +++ b/species/plot/plot_comparison.py @@ -389,8 +389,8 @@ def plot_grid_statistic(tag: str, ax.xaxis.set_minor_locator(AutoMinorLocator(5)) ax.yaxis.set_minor_locator(AutoMinorLocator(5)) - ax.set_xlabel(r'T$_\mathrm{eff}$ (K)', fontsize=13.) - ax.set_ylabel(r'$\mathregular{log}\,g$', fontsize=13.) + ax.set_xlabel(r'T$_\mathregular{eff}$ (K)', fontsize=13.) + ax.set_ylabel(r'$\mathregular{log}\,\mathregular{g}$', fontsize=13.) if xlim is not None: ax.set_xlim(xlim[0], xlim[1]) @@ -445,7 +445,7 @@ def plot_grid_statistic(tag: str, ticklocation='right', format='%.1f') cb.ax.tick_params(width=0.8, length=5, labelsize=12, direction='in', color='black') - cb.ax.set_ylabel(r'$\Delta\mathregular{log}\,G_k$', rotation=270, labelpad=22, fontsize=13.) + cb.ax.set_ylabel(r'$\Delta\mathregular{log}\,\mathregular{G}_\mathregular{k}$', rotation=270, labelpad=22, fontsize=13.) if len(coord_points[2]) != 1: extra_interp = RegularGridInterpolator((coord_points[1], coord_points[0]), extra_map) From 1dada33c9547f5bb2721687caf3c7d5d3d171856 Mon Sep 17 00:00:00 2001 From: Tomas Stolker Date: Wed, 26 May 2021 17:32:59 +0200 Subject: [PATCH 6/7] Bug fix with missing import of warnings --- species/analysis/emission_line.py | 4 +++- species/analysis/fit_model.py | 3 ++- 2 files changed, 5 insertions(+), 2 deletions(-) diff --git a/species/analysis/emission_line.py b/species/analysis/emission_line.py index 18e95b44..b432542c 100644 --- a/species/analysis/emission_line.py +++ b/species/analysis/emission_line.py @@ -4,6 +4,7 @@ import configparser import os +import warnings from typing import Dict, Optional, Tuple, Union @@ -14,7 +15,8 @@ try: import ultranest except: - warnings.warn('UltraNest could not be imported.') + warnings.warn('UltraNest could not be imported. ' + 'Perhaps because cython was not correctly compiled?') from astropy import units as u from astropy.modeling.fitting import LinearLSQFitter diff --git a/species/analysis/fit_model.py b/species/analysis/fit_model.py index 69d56de7..355abdea 100644 --- a/species/analysis/fit_model.py +++ b/species/analysis/fit_model.py @@ -16,7 +16,8 @@ try: import ultranest except: - warnings.warn('UltraNest could not be imported.') + warnings.warn('UltraNest could not be imported. ' + 'Perhaps because cython was not correctly compiled?') # Installation of MultiNest is not possible on readthedocs try: From 4aa085b1e4cac7bbf72d869af9b05e1ed16654fe Mon Sep 17 00:00:00 2001 From: Tomas Stolker Date: Mon, 28 Jun 2021 11:02:18 +0200 Subject: [PATCH 7/7] Added the mean_wavel attribute to ObjectBox --- species/core/box.py | 2 ++ species/data/database.py | 6 ++++++ 2 files changed, 8 insertions(+) diff --git a/species/core/box.py b/species/core/box.py index 8107197d..53140ccd 100644 --- a/species/core/box.py +++ b/species/core/box.py @@ -73,6 +73,7 @@ def create_box(boxtype, box = ObjectBox() box.name = kwargs['name'] box.filters = kwargs['filters'] + box.mean_wavel = kwargs['mean_wavel'] box.magnitude = kwargs['magnitude'] box.flux = kwargs['flux'] box.distance = kwargs['distance'] @@ -317,6 +318,7 @@ def __init__(self): self.name = None self.filters = None + self.mean_wavel = None self.magnitude = None self.flux = None self.distance = None diff --git a/species/data/database.py b/species/data/database.py index 13bd30f6..8cd1905e 100644 --- a/species/data/database.py +++ b/species/data/database.py @@ -1782,6 +1782,7 @@ def get_object(self, magnitude = {} flux = {} + mean_wavel = {} for observatory in dset.keys(): if observatory not in ['distance', 'spectrum']: @@ -1792,6 +1793,9 @@ def get_object(self, magnitude[name] = dset[name][0:2] flux[name] = dset[name][2:4] + filter_trans = read_filter.ReadFilter(name) + mean_wavel[name] = filter_trans.mean_wavelength() + phot_filters = list(magnitude.keys()) else: @@ -1799,6 +1803,7 @@ def get_object(self, magnitude = None flux = None phot_filters = None + mean_wavel = None if inc_spec and f'objects/{object_name}/spectrum' in h5_file: spectrum = {} @@ -1830,6 +1835,7 @@ def get_object(self, return box.create_box('object', name=object_name, filters=phot_filters, + mean_wavel=mean_wavel, magnitude=magnitude, flux=flux, distance=distance,