From aef54ce63c8ca607886d3b09e8a77d42f6da3e2a Mon Sep 17 00:00:00 2001 From: Tomas Stolker Date: Mon, 20 Jan 2020 11:09:32 +0100 Subject: [PATCH] Support of FitPlanck for multiple Planck functions, improved documentation --- docs/running.rst | 2 +- species/analysis/fit_model.py | 66 +++---------- species/analysis/fit_planck.py | 153 +++++++++++++++++++------------ species/analysis/fit_spectrum.py | 16 ++-- species/core/box.py | 4 +- species/data/companions.py | 4 +- species/data/database.py | 48 +++++----- species/data/vlm_plx.py | 14 +-- species/read/read_model.py | 6 +- species/read/read_planck.py | 51 ++++++++--- species/util/phot_util.py | 21 +++-- species/util/read_util.py | 16 ++-- test/test_read/test_color.py | 63 +++++++++++++ test/test_read/test_planck.py | 10 +- 14 files changed, 283 insertions(+), 191 deletions(-) create mode 100644 test/test_read/test_color.py diff --git a/docs/running.rst b/docs/running.rst index 5d00361a..43698b7a 100644 --- a/docs/running.rst +++ b/docs/running.rst @@ -5,7 +5,7 @@ Running species To get a first impression, this page shows what a typical workflow with `species` looks like. -First, the atmospheric model spectra will be downloaded and added to the database. Available photometry of PZ Tel B will be read from :class:`~species.data.companions` dictionary and added to the database. Alternatively, the :func:`~species.data.database.Database.add_object` function can be used for manually creating an planet or other type of object in the database. Then, the model spectra are read from the database, interpolated at the provided parameter values, and stored in a :class:`~species.core.box.ModelBox`. The data of PZ Tel B is read from the database and stored in a :class:`~species.core.box.ObjectBox`. The :class:`~species.core.box.Box` objects are given to the :func:`~species.plot.plot_spectrum.plot_spectrum` function, together with the filter names for which the photometry is available. +First, the atmospheric model spectra will be downloaded and added to the database. Available photometry of PZ Tel B will be read from the :class:`~species.data.companions` dictionary and added to the database. Alternatively, the :func:`~species.data.database.Database.add_object` function can be used for manually creating an planet or other type of object in the database. Then, the model spectra are read from the database, interpolated at the provided parameter values, and stored in a :class:`~species.core.box.ModelBox`. The data of PZ Tel B is read from the database and stored in a :class:`~species.core.box.ObjectBox`. The :class:`~species.core.box.Box` objects are given to the :func:`~species.plot.plot_spectrum.plot_spectrum` function, together with the filter names for which the photometry is available. The following code can be executed from the command line or within a `Jupyter notebook `_. diff --git a/species/analysis/fit_model.py b/species/analysis/fit_model.py index 26f0926e..d7c94371 100644 --- a/species/analysis/fit_model.py +++ b/species/analysis/fit_model.py @@ -1,5 +1,5 @@ """ -Module for fitting atmospheric models. +Module with functionalities for fitting atmospheric model specra. """ import sys @@ -22,7 +22,7 @@ def lnprior(param, modelpar, prior=None): """ - Function for the prior probability. + Internal function for the prior probability. Parameters ---------- @@ -76,11 +76,10 @@ def lnlike(param, synphot, distance, spectrum, - instrument, modelspec, weighting): """ - Function for the likelihood probability. + Internal function for the likelihood probability. Parameters ---------- @@ -95,8 +94,6 @@ def lnlike(param, Distance (pc). spectrum : numpy.ndarray Wavelength (micron), apparent flux (W m-2 micron-1), and flux error (W m-2 micron-1). - instrument : str - Instrument that was used for the spectrum (currently only 'gpi' possible). modelspec : species.read.read_model.ReadModel weighting : float, None Weighting applied to the spectrum when calculating the likelihood function in order @@ -151,11 +148,10 @@ def lnprob(param, distance, prior, spectrum, - instrument, modelspec, weighting): """ - Function for the posterior probability. + Internal function for the posterior probability. Parameters ---------- @@ -176,8 +172,6 @@ def lnprob(param, used if set to None. spectrum : numpy.ndarray Wavelength (micron), apparent flux (W m-2 micron-1), and flux error (W m-2 micron-1). - instrument : str - Instrument that was used for the spectrum (currently only 'gpi' possible). modelspec : species.read.read_model.ReadModel weighting : float, None Weighting applied to the spectrum when calculating the likelihood function in order @@ -205,7 +199,6 @@ def lnprob(param, synphot, distance, spectrum, - instrument, modelspec, weighting) @@ -217,11 +210,11 @@ def lnprob(param, class FitModel: """ - Fit atmospheric model spectra to photometric and spectral data. + Class for fitting atmospheric model spectra to photometric data. """ def __init__(self, - objname, + object_name, filters, model, bounds, @@ -230,7 +223,7 @@ def __init__(self, """ Parameters ---------- - objname : str + object_name : str Object name in the database. filters : tuple(str, ) Filter IDs for which the photometry is selected. All available photometry of the @@ -251,7 +244,7 @@ def __init__(self, None """ - self.object = read_object.ReadObject(objname) + self.object = read_object.ReadObject(object_name) self.distance = self.object.get_distance() self.model = model @@ -260,13 +253,8 @@ def __init__(self, if not inc_phot and not inc_spec: raise ValueError('No photometric or spectral data has been selected.') - if self.bounds is not None and 'teff' in self.bounds: - teff_bound = self.bounds['teff'] - else: - teff_bound = None - if self.bounds is not None: - readmodel = read_model.ReadModel(self.model, None, teff_bound) + readmodel = read_model.ReadModel(self.model) bounds_grid = readmodel.get_bounds() for item in bounds_grid: @@ -287,11 +275,11 @@ def __init__(self, if not filters: species_db = database.Database() - objectbox = species_db.get_object(objname, None) - filters = objectbox.filter + objectbox = species_db.get_object(object_name, None) + filters = objectbox.filters for item in filters: - readmodel = read_model.ReadModel(self.model, item, teff_bound) + readmodel = read_model.ReadModel(self.model, filter_name=item) readmodel.interpolate_model() self.modelphot.append(readmodel) @@ -309,7 +297,7 @@ def __init__(self, if inc_spec: self.spectrum = self.object.get_spectrum() self.instrument = self.object.get_instrument() - self.modelspec = read_model.ReadModel(self.model, (0.9, 2.5), teff_bound) + self.modelspec = read_model.ReadModel(self.model, wavel_range=(0.9, 2.5)) else: self.spectrum = None @@ -374,7 +362,7 @@ def run_mcmc(self, high=self.bounds[item][1], size=nwalkers) - with Pool(processes=cpu_count()) as pool: + with Pool(processes=cpu_count()): sampler = emcee.EnsembleSampler(nwalkers, ndim, lnprob, @@ -386,37 +374,11 @@ def run_mcmc(self, self.distance, prior, self.spectrum, - self.instrument, self.modelspec, weighting])) sampler.run_mcmc(initial, nsteps, progress=True) - # sampler = emcee.EnsembleSampler(nwalkers=nwalkers, - # dim=ndim, - # lnpostfn=lnprob, - # a=2., - # args=([self.bounds, - # self.modelpar, - # self.modelphot, - # self.objphot, - # self.synphot, - # self.distance, - # prior, - # self.spectrum, - # self.instrument, - # self.modelspec, - # weighting])) - - # progbar = progress.bar.Bar('\rRunning MCMC...', - # max=nsteps, - # suffix='%(percent)d%%') - - # for i, _ in enumerate(sampler.sample(initial, iterations=nsteps)): - # progbar.next() - - # progbar.finish() - species_db = database.Database() species_db.add_samples(sampler=sampler, diff --git a/species/analysis/fit_planck.py b/species/analysis/fit_planck.py index 2e12e03b..dce3b8d8 100644 --- a/species/analysis/fit_planck.py +++ b/species/analysis/fit_planck.py @@ -1,5 +1,5 @@ """ -Module for fitting atmospheric models. +Module with functionalities for fitting a Planck spectrum. """ import sys @@ -17,19 +17,18 @@ def lnprior(param, - bounds, - modelpar): + bounds): """ - Function for the prior probability. + Internal function for the prior probability. Parameters ---------- param : numpy.ndarray Parameter values. bounds : dict - Parameter boundaries. - modelpar : list(str, ) - Parameter names. + Parameter boundaries for 'teff' and 'radius'. The values should be provided in a list + such that multiple Planck functions can be combined, e.g. ``{'teff': [(1000., 2000.), + (500., 1500.)], 'radius': [(0.5, 1.5), (1.5, 2.0)]}``. Returns ------- @@ -39,8 +38,7 @@ def lnprior(param, ln_prior = 0 - for i, item in enumerate(modelpar): - + for i, item in enumerate(bounds): if bounds[item][0] <= param[i] <= bounds[item][1]: ln_prior += 0. @@ -52,30 +50,33 @@ def lnprior(param, def lnlike(param, - modelpar, + bounds, objphot, synphot, distance, spectrum, - instrument, weighting): """ - Function for the likelihood probability. + Internal function for the likelihood probability. Parameters ---------- param : numpy.ndarray Parameter values. - modelpar : list(str, ) - Parameter names. + bounds : dict + Parameter boundaries for 'teff' and 'radius'. The values should be provided in a list + such that multiple Planck functions can be combined, e.g. ``{'teff': [(1000., 2000.), + (500., 1500.)], 'radius': [(0.5, 1.5), (1.5, 2.0)]}``. objphot : list(tuple(float, float), ) + List with the photometric fluxes and uncertainties. synphot : list(species.analysis.photometry.SyntheticPhotometry, ) + List with the :class:`~species.analysis.photometry.SyntheticPhotometry` objects for + calculation of synthetic photometry from the model spectra. distance : float Distance (pc). - spectrum : numpy.ndarray - Wavelength (micron), apparent flux (W m-2 micron-1), and flux error (W m-2 micron-1). - instrument : str - Instrument that was used for the spectrum (currently only 'gpi' possible). + spectrum : numpy.ndarray, None + Spectrum array with the wavelength (micron), flux (W m-2 micron-1), and error + (W m-2 micron-1). Not used if set to None. weighting : float, None Weighting applied to the spectrum when calculating the likelihood function in order to not have a spectrum dominate the chi-squared value. For example, with `weighting=3` @@ -90,7 +91,7 @@ def lnlike(param, """ paramdict = {} - for i, item in enumerate(modelpar): + for i, item in enumerate(bounds.keys()): paramdict[item] = param[i] paramdict['distance'] = distance @@ -98,15 +99,15 @@ def lnlike(param, chisq = 0. if objphot is not None: - for i, item in enumerate(objphot): - readplanck = read_planck.ReadPlanck(synphot[i].filter_name) - flux = readplanck.get_photometry(paramdict, synphot=synphot[i]) + for i, obj_item in enumerate(objphot): + readplanck = read_planck.ReadPlanck(filter_name=synphot[i].filter_name) + flux = readplanck.get_flux(paramdict, synphot=synphot[i]) - chisq += (item[0]-flux)**2 / item[1]**2 + chisq += (obj_item[0]-flux)**2 / obj_item[1]**2 if spectrum is not None: - # TODO check if the wavelength range of get_planck is broad enought for spectres - readplanck = read_planck.ReadPlanck((spectrum[0, 0], spectrum[-1, 0])) + readplanck = read_planck.ReadPlanck((0.9*spectrum[0, 0], 1.1*spectrum[-1, 0])) + model = readplanck.get_spectrum(paramdict, 100.) flux_new = spectres.spectres(new_spec_wavs=spectrum[:, 0], @@ -126,32 +127,32 @@ def lnlike(param, def lnprob(param, bounds, - modelpar, objphot, synphot, distance, spectrum, - instrument, weighting): """ - Function for the posterior probability. + Internal function for the posterior probability. Parameters ---------- param : numpy.ndarray Parameter values. bounds : dict - Parameter boundaries. - modelpar : list(str, ) - Parameter names. + Parameter boundaries for 'teff' and 'radius'. The values should be provided in a list + such that multiple Planck functions can be combined, e.g. ``{'teff': [(1000., 2000.), + (500., 1500.)], 'radius': [(0.5, 1.5), (1.5, 2.0)]}``. objphot : list(tuple(float, float), ) + List with the photometric fluxes and uncertainties. synphot : list(species.analysis.photometry.SyntheticPhotometry, ) + List with the :class:`~species.analysis.photometry.SyntheticPhotometry` objects for + calculation of synthetic photometry from the model spectra. distance : float Distance (pc). - spectrum : numpy.ndarray - Wavelength (micron), apparent flux (W m-2 micron-1), and flux error (W m-2 micron-1). - instrument : str - Instrument that was used for the spectrum (currently only 'gpi' possible). + spectrum : numpy.ndarray, None + Spectrum array with the wavelength (micron), flux (W m-2 micron-1), and error + (W m-2 micron-1). Not used if set to None. weighting : float, None Weighting applied to the spectrum when calculating the likelihood function in order to not have a spectrum dominate the chi-squared value. For example, with `weighting=3` @@ -165,19 +166,18 @@ def lnprob(param, Log posterior probability. """ - ln_prior = lnprior(param, bounds, modelpar) + ln_prior = lnprior(param, bounds) if math.isinf(ln_prior): ln_prob = -np.inf else: ln_prob = ln_prior + lnlike(param, - modelpar, + bounds, objphot, synphot, distance, spectrum, - instrument, weighting) if np.isnan(ln_prob): @@ -188,11 +188,11 @@ def lnprob(param, class FitPlanck: """ - Fit Planck spectrum to photometric and spectral data. + Class for fitting Planck spectra to spectral and photometric data. """ def __init__(self, - objname, + object_name, filters, bounds, inc_phot=True, @@ -200,13 +200,15 @@ def __init__(self, """ Parameters ---------- - objname : str + object_name : str Object name in the database. filters : tuple(str, ) - Filter IDs for which the photometry is selected. All available photometry of the - object is selected if set to None. + Filter names for which the photometry is selected. All available photometry of the + object are used if set to None. bounds : dict - Parameter boundaries for 'teff' and 'radius'. + Parameter boundaries for 'teff' and 'radius'. The values should be provided either as + float or as list of floats such that multiple Planck functions can be combined, + e.g. ``{'teff': [(1000., 2000.), (500., 1500.)], 'radius': [(0.5, 1.5), (1.5, 2.0)]}``. inc_phot : bool Include photometry data with the fit. inc_spec : bool @@ -218,27 +220,40 @@ def __init__(self, None """ - self.object = read_object.ReadObject(objname) - self.distance = self.object.get_distance() - - self.model = 'planck' - self.bounds = bounds - self.modelpar = ['teff', 'radius'] - if not inc_phot and not inc_spec: raise ValueError('No photometric or spectral data has been selected.') - if 'teff' not in self.bounds or 'radius' not in self.bounds: + if 'teff' not in bounds or 'radius' not in bounds: raise ValueError('The \'bounds\' dictionary should contain \'teff\' and \'radius\'.') + self.model = 'planck' + + self.object = read_object.ReadObject(object_name) + self.distance = self.object.get_distance() + + if isinstance(bounds['teff'], list) and isinstance(bounds['radius'], list): + self.modelpar = [] + self.bounds = {} + + for i, item in enumerate(bounds['teff']): + self.modelpar.append(f'teff_{i}') + self.modelpar.append(f'radius_{i}') + + self.bounds[f'teff_{i}'] = bounds['teff'][i] + self.bounds[f'radius_{i}'] = bounds['radius'][i] + + else: + self.modelpar = ['teff', 'radius'] + self.bounds = bounds + if inc_phot: self.synphot = [] self.objphot = [] if not filters: species_db = database.Database() - objectbox = species_db.get_object(objname, None) - filters = objectbox.filter + objectbox = species_db.get_object(object_name, None) + filters = objectbox.filters for item in filters: sphot = photometry.SyntheticPhotometry(item) @@ -276,7 +291,9 @@ def run_mcmc(self, Number of steps per walker. guess : dict, None Guess for the 'teff' and 'radius'. Random values between the boundary values are used - if a dictionary value is set to None. + if a value is set to None. The values should be provided either as float or in a list + of floats such that multiple Planck functions can be combined, e.g. + ``{'teff': [1500., 1000.], 'radius': [1., 2.]``. tag : str Database tag where the MCMC samples are stored. weighting : float, None @@ -292,14 +309,30 @@ def run_mcmc(self, None """ - sigma = {'teff': 5., 'radius': 0.01} - sys.stdout.write('Running MCMC...') sys.stdout.flush() - ndim = 2 + ndim = len(self.bounds) + + if ndim == 2: + sigma = {'teff': 5., 'radius': 0.01} + + else: + sigma = {} + + for i, item in enumerate(guess['teff']): + sigma[f'teff_{i}'] = 5. + guess[f'teff_{i}'] = guess['teff'][i] + + for i, item in enumerate(guess['radius']): + sigma[f'radius_{i}'] = 0.01 + guess[f'radius_{i}'] = guess['radius'][i] + + del guess['teff'] + del guess['radius'] initial = np.zeros((nwalkers, ndim)) + for i, item in enumerate(self.modelpar): if guess[item] is not None: initial[:, i] = guess[item] + np.random.normal(0, sigma[item], nwalkers) @@ -309,17 +342,15 @@ def run_mcmc(self, high=self.bounds[item][1], size=nwalkers) - with Pool(processes=cpu_count()) as pool: + with Pool(processes=cpu_count()): sampler = emcee.EnsembleSampler(nwalkers, ndim, lnprob, args=([self.bounds, - self.modelpar, self.objphot, self.synphot, self.distance, self.spectrum, - self.instrument, weighting])) sampler.run_mcmc(initial, nsteps, progress=True) diff --git a/species/analysis/fit_spectrum.py b/species/analysis/fit_spectrum.py index f6f6569b..b95ced6c 100644 --- a/species/analysis/fit_spectrum.py +++ b/species/analysis/fit_spectrum.py @@ -1,5 +1,5 @@ """ -Module for fitting a calibration spectrum. +Module with functionalities for fitting a calibration spectrum. """ import sys @@ -23,6 +23,8 @@ def lnprob(param, specphot, bands): """ + Internal function for the posterior probability. + Parameters ---------- param : numpy.ndarray @@ -75,18 +77,18 @@ def lnprob(param, class FitSpectrum: """ - Fit a calibration spectrum to photometric data. + Class for fitting a calibration spectrum to photometric data. """ def __init__(self, - objname, + object_name, filters, spectrum, bounds): """ Parameters ---------- - objname : str + object_name : str Object name in the database. filters : tuple(str, ) Filter IDs for which the photometry is selected. All available photometry of the @@ -101,7 +103,7 @@ def __init__(self, None """ - self.object = read_object.ReadObject(objname) + self.object = read_object.ReadObject(object_name) self.spectrum = spectrum self.bounds = bounds @@ -111,8 +113,8 @@ def __init__(self, if filters is None: species_db = database.Database() - objectbox = species_db.get_object(objname, None) - filters = objectbox.filter + objectbox = species_db.get_object(object_name, None) + filters = objectbox.filters for item in filters: readcalib = read_calibration.ReadCalibration(self.spectrum, item) diff --git a/species/core/box.py b/species/core/box.py index a5cb1dae..9c7c461d 100644 --- a/species/core/box.py +++ b/species/core/box.py @@ -54,7 +54,7 @@ def create_box(boxtype, elif boxtype == 'object': box = ObjectBox() box.name = kwargs['name'] - box.filter = kwargs['filter'] + box.filters = kwargs['filters'] box.magnitude = kwargs['magnitude'] box.flux = kwargs['flux'] box.distance = kwargs['distance'] @@ -222,7 +222,7 @@ def __init__(self): """ self.name = None - self.filter = None + self.filters = None self.magnitude = None self.flux = None self.distance = None diff --git a/species/data/companions.py b/species/data/companions.py index d68d354d..c2565b91 100644 --- a/species/data/companions.py +++ b/species/data/companions.py @@ -97,12 +97,12 @@ def get_data(): 'Paranal/SPHERE.IRDIS_D_H23_3': (17.95, 0.17), # Keppler et al. 2018 'Paranal/SPHERE.IRDIS_D_K12_1': (16.65, 0.06), # Müller et al. 2018 'Paranal/SPHERE.IRDIS_D_K12_2': (16.44, 0.05), # Müller et al. 2018 - 'Paranal/NACO.Lp': (14.61, 0.78), # Stolker et al. in prep. + 'Paranal/NACO.Lp': (14.46, 0.50), # Stolker et al. in prep. 'Paranal/NACO.NB405': (14.15, 0.21), # Stolker et al. in prep. 'Paranal/NACO.Mp': (13.52, 0.26)}}, # Stolker et al. in prep. 'PDS 70 c': {'distance': (113.43, 0.52), - 'app_mag': {'Paranal/NACO.NB405': (14.93, 0.57)}}, # Stolker et al. in prep. + 'app_mag': {'Paranal/NACO.NB405': (15.05, 0.59)}}, # Stolker et al. in prep. '2M1207 b': {'distance': (64.42, 0.65), 'app_mag': {'HST/NICMOS1.F090M': (22.58, 0.35), # Song et al. 2006 diff --git a/species/data/database.py b/species/data/database.py index cc6e4210..503b1a2f 100644 --- a/species/data/database.py +++ b/species/data/database.py @@ -1,5 +1,5 @@ """ -Database module. +Module with functionalities for reading and writing of data. """ import os @@ -23,7 +23,7 @@ class Database: """ - Text. + Class for fitting atmospheric model spectra to photometric data. """ def __init__(self): @@ -618,6 +618,10 @@ def add_samples(self, except emcee.autocorr.AutocorrError: int_auto = None + sys.stdout.write('The chain is shorter than 50 times the integrated autocorrelation ' + 'time. [WARNING]\n') + sys.stdout.flush() + if int_auto is not None: for i, item in enumerate(int_auto): dset.attrs['autocorrelation'+str(i)] = float(item) @@ -764,8 +768,8 @@ def get_mcmc_spectra(self, spectrum_name = dset.attrs['spectrum'] if spec_res is not None and spectrum_type == 'calibration': - warnings.warn("Smoothing of the spectral resolution is not implemented for calibration " - "spectra.") + warnings.warn('Smoothing of the spectral resolution is not implemented for calibration ' + 'spectra.') if dset.attrs.__contains__('distance'): distance = dset.attrs['distance'] @@ -799,21 +803,21 @@ def get_mcmc_spectra(self, suffix='%(percent)d%%') for i in range(samples.shape[0]): - model_par = {} + model_param = {} for j in range(samples.shape[1]): - model_par[param[j]] = samples[i, j] + model_param[param[j]] = samples[i, j] if distance: - model_par['distance'] = distance + model_param['distance'] = distance if spectrum_type == 'model': if spectrum_name == 'planck': - specbox = readmodel.get_spectrum(model_par, spec_res) + specbox = readmodel.get_spectrum(model_param, spec_res) else: - specbox = readmodel.get_model(model_par, spec_res) + specbox = readmodel.get_model(model_param, spec_res) elif spectrum_type == 'calibration': - specbox = readcalib.get_spectrum(model_par) + specbox = readcalib.get_spectrum(model_param) box.type = 'mcmc' @@ -881,17 +885,17 @@ def get_mcmc_photometry(self, suffix='%(percent)d%%') for i in range(samples.shape[0]): - model_par = {} + model_param = {} for j in range(nparam): - model_par[param[j]] = samples[i, j] + model_param[param[j]] = samples[i, j] if distance: - model_par['distance'] = distance + model_param['distance'] = distance if spectrum_type == 'model': - mcmc_phot[i, 0], _ = readmodel.get_magnitude(model_par) + mcmc_phot[i, 0], _ = readmodel.get_magnitude(model_param) # elif spectrum_type == 'calibration': - # specbox = readcalib.get_spectrum(model_par) + # specbox = readcalib.get_spectrum(model_param) progbar.next() @@ -901,7 +905,7 @@ def get_mcmc_photometry(self, def get_object(self, object_name, - filter_id=None, + filters=None, inc_phot=True, inc_spec=True): """ @@ -909,7 +913,7 @@ def get_object(self, ---------- object_name : str Object name in the database. - filter_id : tuple(str, ) + filters : tuple(str, ) Filter IDs for which the photometry is selected. All available photometry of the object is selected if set to None. inc_phot : bool @@ -936,8 +940,8 @@ def get_object(self, magnitude = {} flux = {} - if filter_id: - for item in filter_id: + if filters: + for item in filters: data = dset[item] magnitude[item] = np.asarray(data[0:2]) @@ -952,13 +956,13 @@ def get_object(self, magnitude[name] = np.asarray(dset[name][0:2]) flux[name] = np.asarray(dset[name][2:4]) - filterids = tuple(magnitude.keys()) + filters = tuple(magnitude.keys()) else: magnitude = None flux = None - filterids = None + filters = None if inc_spec and 'objects/'+object_name+'/spectrum' in h5_file: spectrum = np.asarray(h5_file['objects/'+object_name+'/spectrum']) @@ -972,7 +976,7 @@ def get_object(self, return box.create_box('object', name=object_name, - filter=filterids, + filters=filters, magnitude=magnitude, flux=flux, distance=distance, diff --git a/species/data/vlm_plx.py b/species/data/vlm_plx.py index 1b2210dd..7c67a05c 100644 --- a/species/data/vlm_plx.py +++ b/species/data/vlm_plx.py @@ -101,15 +101,15 @@ def add_vlm_plx(input_path, sys.stdout.write(' [DONE]\n') sys.stdout.flush() - sys.stdout.write('Querying SIMBAD...') - sys.stdout.flush() + # sys.stdout.write('Querying SIMBAD...') + # sys.stdout.flush() - simbad_id = queries.get_simbad(name) + # simbad_id = queries.get_simbad(name) - dset = database.create_dataset(group+'/simbad', (np.size(simbad_id), ), dtype=dtype) - dset[...] = simbad_id + # dset = database.create_dataset(group+'/simbad', (np.size(simbad_id), ), dtype=dtype) + # dset[...] = simbad_id - sys.stdout.write(' [DONE]\n') - sys.stdout.flush() + # sys.stdout.write(' [DONE]\n') + # sys.stdout.flush() database.close() diff --git a/species/read/read_model.py b/species/read/read_model.py index 128c2904..bac7574c 100644 --- a/species/read/read_model.py +++ b/species/read/read_model.py @@ -33,9 +33,9 @@ def __init__(self, ---------- model : str Name of the atmospheric model. - wavel_range : tuple(float, float), str, None - Wavelength range (micron) or filter ID. Full spectrum is selected if set to None. Not - used if ``filter_name`` is not None. + wavel_range : tuple(float, float), None + Wavelength range (micron). Full spectrum is selected if set to None. Not used if + ``filter_name`` is not None. filter_name : str, None Filter ID that is used for the wavelength range. The ``wavel_range`` is used if set to None. diff --git a/species/read/read_planck.py b/species/read/read_planck.py index 181e789c..dee0993d 100644 --- a/species/read/read_planck.py +++ b/species/read/read_planck.py @@ -19,12 +19,17 @@ class ReadPlanck: """ def __init__(self, - filter_name): + wavel_range=None, + filter_name=None): """ Parameters ---------- + wavel_range : tuple(float, float), None + Wavelength range (micron). A wavelength range of 0.1-1000 micron is used if set to + None. Not used if ``filter_name`` is not None. filter_name : str, None - Filter ID that is used for the wavelength range. Full spectrum is used if set to None. + Filter name that is used for the wavelength range. The ``wavel_range`` is used if set + to None. Returns ------- @@ -36,14 +41,15 @@ def __init__(self, self.wl_points = None self.wl_index = None - if isinstance(filter_name, str): - self.filter_name = filter_name - transmission = read_filter.ReadFilter(filter_name) + self.filter_name = filter_name + self.wavel_range = wavel_range + + if self.filter_name is not None: + transmission = read_filter.ReadFilter(self.filter_name) self.wavel_range = transmission.wavelength_range() - else: - self.filter_name = None - self.wavel_range = filter_name + elif self.wavel_range is None: + self.wavel_range = (0.1, 1000.) config_file = os.path.join(os.getcwd(), 'species_config.ini') @@ -85,12 +91,15 @@ def get_spectrum(self, model_param, spec_res): """ - Function for calculating a Planck spectrum. + Function for calculating a Planck spectrum or a combination of multiple Planck spectra. Parameters ---------- model_param : dict - Dictionary with the 'teff' (K), 'radius' (Rjup), and 'distance' (pc). + Dictionary with the 'teff' (K), 'radius' (Rjup), and 'distance' (pc). The values of + 'teff' and 'radius' can be a single float, or a list with floats for a combination of + multiple Planck functions, e.g. ``{'teff': [1500., 1000.], 'radius': [1., 2.], + 'distance': 10.}``. spec_res : float Spectral resolution. @@ -107,10 +116,26 @@ def get_spectrum(self, wavel_points = np.asarray(wavel_points) # [micron] - scaling = ((model_param['radius']*constants.R_JUP) / - (model_param['distance']*constants.PARSEC))**2 + n_planck = (len(model_param)-1) // 2 + + if n_planck == 1: + scaling = ((model_param['radius']*constants.R_JUP) / + (model_param['distance']*constants.PARSEC))**2 + + flux = self.planck(np.copy(wavel_points), + model_param['teff'], + scaling) # [W m-2 micron-1] + + else: + flux = np.zeros(wavel_points.shape) + + for i in range(n_planck): + scaling = ((model_param[f'radius_{i}']*constants.R_JUP) / + (model_param['distance']*constants.PARSEC))**2 - flux = self.planck(np.copy(wavel_points), model_param['teff'], scaling) # [W m-2 micron-1] + flux += self.planck(np.copy(wavel_points), + model_param[f'teff_{i}'], + scaling) # [W m-2 micron-1] return box.create_box(boxtype='model', model='planck', diff --git a/species/util/phot_util.py b/species/util/phot_util.py index f1788b9f..039d60a2 100644 --- a/species/util/phot_util.py +++ b/species/util/phot_util.py @@ -41,16 +41,16 @@ def multi_photometry(datatype, if datatype == 'model': for item in filters: if spectrum == 'planck': - readmodel = read_planck.ReadPlanck(item) + readmodel = read_planck.ReadPlanck(filter_name=item) else: - readmodel = read_model.ReadModel(spectrum, item) + readmodel = read_model.ReadModel(spectrum, filter_name=item) - flux[item] = readmodel.get_photometry(parameters) + flux[item] = readmodel.get_flux(parameters) elif datatype == 'calibration': for item in filters: - readcalib = read_calibration.ReadCalibration(spectrum, item) - flux[item] = readcalib.get_photometry(parameters) + readcalib = read_calibration.ReadCalibration(spectrum, filter_name=item) + flux[item] = readcalib.get_flux(parameters) sys.stdout.write(' [DONE]\n') sys.stdout.flush() @@ -109,7 +109,7 @@ def get_residuals(datatype, """ if filters is None: - filters = objectbox.filter + filters = objectbox.filters if inc_phot: model_phot = multi_photometry(datatype=datatype, @@ -134,8 +134,13 @@ def get_residuals(datatype, if inc_spec: wavel_range = (0.9*objectbox.spectrum[0, 0], 1.1*objectbox.spectrum[-1, 0]) - readmodel = read_model.ReadModel(spectrum, wavel_range) - model = readmodel.get_model(parameters) + if spectrum == 'planck': + readmodel = read_planck.ReadPlanck(wavel_range=wavel_range) + model = readmodel.get_spectrum(model_param=parameters, spec_res=1000.) + + else: + readmodel = read_model.ReadModel(spectrum, wavel_range=wavel_range) + model = readmodel.get_model(model_param=parameters, spec_res=1000.) wl_new = objectbox.spectrum[:, 0] diff --git a/species/util/read_util.py b/species/util/read_util.py index bdea8738..183d61cd 100644 --- a/species/util/read_util.py +++ b/species/util/read_util.py @@ -12,11 +12,11 @@ from species.read import read_model, read_planck -def get_mass(model_par): +def get_mass(model_param): """ Parameters ---------- - model_par : dict + model_param : dict Model parameter values. Should contain the surface gravity and radius. Returns @@ -25,9 +25,9 @@ def get_mass(model_par): Mass (Mjup). """ - logg = 1e-2 * 10.**model_par['logg'] # [m s-1] + logg = 1e-2 * 10.**model_param['logg'] # [m s-1] - radius = model_par['radius'] # [Rjup] + radius = model_param['radius'] # [Rjup] radius *= constants.R_JUP # [m] mass = logg*radius**2/constants.GRAVITY # [kg] @@ -54,12 +54,12 @@ def add_luminosity(modelbox): """ if modelbox.model == 'planck': - readmodel = read_planck.ReadPlanck(wavelength=(1e-1, 1e3)) - fullspec = readmodel.get_spectrum(model_par=modelbox.parameters, spec_res=1000.) + readmodel = read_planck.ReadPlanck(wavel_range=(1e-1, 1e3)) + fullspec = readmodel.get_spectrum(model_param=modelbox.parameters, spec_res=1000.) else: - readmodel = read_model.ReadModel(model=modelbox.model, wavelength=None, teff=None) - fullspec = readmodel.get_model(model_par=modelbox.parameters) + readmodel = read_model.ReadModel(model=modelbox.model) + fullspec = readmodel.get_model(model_param=modelbox.parameters) flux = simps(fullspec.flux, fullspec.wavelength) diff --git a/test/test_read/test_color.py b/test/test_read/test_color.py new file mode 100644 index 00000000..d8695c70 --- /dev/null +++ b/test/test_read/test_color.py @@ -0,0 +1,63 @@ +import os +import shutil +import pytest + +import numpy as np + +import species +from species.util import test_util + + +class TestColor: + + def setup_class(self): + self.limit = 1e-10 + + def teardown_class(self): + os.remove('species_database.hdf5') + os.remove('species_config.ini') + shutil.rmtree('data/') + + def test_species_init(self): + test_util.create_config('./') + species.SpeciesInit('./') + + def test_read_color_magnitude(self): + database = species.Database() + database.add_photometry('vlm-plx') + database.add_photometry('leggett') + + read_colormag = species.ReadColorMagnitude(['vlm-plx', 'leggett'], + ('MKO/NSFCam.J', 'MKO/NSFCam.H'), + 'MKO/NSFCam.J') + + assert read_colormag.filters_color == ('MKO/NSFCam.J', 'MKO/NSFCam.H') + assert read_colormag.filter_mag == 'MKO/NSFCam.J' + + def test_get_color_magnitude(self): + read_colormag = species.ReadColorMagnitude(['vlm-plx', 'leggett'], + ('MKO/NSFCam.J', 'MKO/NSFCam.H'), + 'MKO/NSFCam.J') + + colormag_box = read_colormag.get_color_magnitude(object_type=None) + + assert np.nansum(colormag_box.color) == pytest.approx(234.09106) + assert np.nansum(colormag_box.magnitude) == pytest.approx(6794.118) + + def test_read_color_color(self): + read_colorcolor = species.ReadColorColor(['vlm-plx', ], + (('MKO/NSFCam.J', 'MKO/NSFCam.H'), + ('MKO/NSFCam.H', 'MKO/NSFCam.K'))) + + assert read_colorcolor.filters_colors == (('MKO/NSFCam.J', 'MKO/NSFCam.H'), + ('MKO/NSFCam.H', 'MKO/NSFCam.K')) + + def test_get_color_color(self): + read_colorcolor = species.ReadColorColor(['vlm-plx', ], + (('MKO/NSFCam.J', 'MKO/NSFCam.H'), + ('MKO/NSFCam.H', 'MKO/NSFCam.K'))) + + colorcolor_box = read_colorcolor.get_color_color(object_type=None) + + assert np.nansum(colorcolor_box.color1) == pytest.approx(162.13107) + assert np.nansum(colorcolor_box.color2) == pytest.approx(141.24994) diff --git a/test/test_read/test_planck.py b/test/test_read/test_planck.py index 732a4523..d57900da 100644 --- a/test/test_read/test_planck.py +++ b/test/test_read/test_planck.py @@ -24,14 +24,14 @@ def test_species_init(self): species.SpeciesInit('./') def test_read_planck(self): - read_planck = species.ReadPlanck('MKO/NSFCam.J') + read_planck = species.ReadPlanck(filter_name='MKO/NSFCam.J') assert read_planck.wavel_range == pytest.approx((1.1308, 1.3812)) - read_planck = species.ReadPlanck((1., 5.)) + read_planck = species.ReadPlanck(wavel_range=(1., 5.)) assert read_planck.wavel_range == (1., 5.) def test_get_spectrum(self): - read_planck = species.ReadPlanck('MKO/NSFCam.J') + read_planck = species.ReadPlanck(filter_name='MKO/NSFCam.J') modelbox = read_planck.get_spectrum({'teff': 2000., 'radius': 1., 'distance': 10.}, 100.) assert modelbox.model == 'planck' @@ -42,11 +42,11 @@ def test_get_spectrum(self): assert np.allclose(np.sum(modelbox.flux), 1.827363592502394e-12, rtol=self.limit, atol=0.) def test_get_flux(self): - read_planck = species.ReadPlanck('MKO/NSFCam.J') + read_planck = species.ReadPlanck(filter_name='MKO/NSFCam.J') flux = read_planck.get_flux({'teff': 2000., 'radius': 1., 'distance': 10.}) assert flux == 8.322445907073985e-14 - synphot = species.SyntheticPhotometry('MKO/NSFCam.J') + synphot = species.SyntheticPhotometry(filter_name='MKO/NSFCam.J') flux = read_planck.get_flux({'teff': 2000., 'radius': 1., 'distance': 10.}, synphot=synphot) assert flux == 8.322445907073985e-14