# species.analysis package#

## species.analysis.compare_spectra module#

Module with functionalities for comparing a spectrum with a library of empirical or model spectra.

class species.analysis.compare_spectra.CompareSpectra(object_name: str, spec_name: Union[str, List[str]], spec_library: Optional[str] = None)[source]#

Bases: object

Class for comparing a spectrum of an object with a library of empirical or model spectra.

Parameters
• object_name (str) – Object name as stored in the database with add_object() or add_companion().

• spec_name (str, list(str)) – Name of the spectrum or list with spectrum names that are stored at the object data of object_name. The argument can be either a string or a list of strings.

• spec_library (str, None) – DEPRECATED: Name of the spectral library (‘irtf’, ‘spex’, or ‘kesseli+2017).

Returns

None

Return type

NoneType

compare_model(tag: str, model: str, av_points: Optional[Union[List[float], numpy.array]] = None, fix_logg: Optional[float] = None, scale_spec: Optional[List[str]] = None, weights: bool = True, inc_phot: Optional[List[str]] = None) None[source]#

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 $$T_\mathrm{eff}$$ and $$\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.

Parameters
• tag (str) – Database tag where for each spectrum from the spectral library the best-fit parameters will be stored. So when testing a range of values for av_ext and rad_vel, only the parameters that minimize the goodness-of-fit statistic will be stored.

• model (str) – Name of the atmospheric model grid with synthetic spectra.

• av_points (list(float), np.array, None) – List of $$A_V$$ extinction values for which the goodness-of-fit statistic will be tested. The extinction is calculated with the relation from Cardelli et al. (1989).

• fix_logg (float, None) – Fix the value of $$\log(g)$$, for example if estimated from gravity-sensitive spectral features. Typically, $$\log(g)$$ can not be accurately determined when comparing the spectra over a broad wavelength range.

• scale_spec (list(str), None) – List with names of observed spectra to which a flux scaling is applied to best match the spectral templates.

• weights (bool) – Apply a weighting based on the widths of the wavelengths bins.

• inc_phot (list(str), None) – Filter names of the photometry to include in the comparison. Photometry points are weighted by the FWHM of the filter profile. No photometric fluxes will be used if the argument is set to None.

Returns

None

Return type

NoneType

spectral_type(tag: str, spec_library, wavel_range: Optional[Tuple[Optional[float], Optional[float]]] = None, sptypes: Optional[List[str]] = None, av_ext: Optional[Union[List[float], numpy.array]] = None, rad_vel: Optional[Union[List[float], numpy.array]] = None) None[source]#

Method for finding the best fitting empirical spectra from a selected library by evaluating the goodness-of-fit statistic from Cushing et al. (2008).

Parameters
• tag (str) – Database tag where for each spectrum from the spectral library the best-fit parameters will be stored. So when testing a range of values for av_ext and rad_vel, only the parameters that minimize the goodness-of-fit statistic will be stored.

• spec_library (str) – Name of the spectral library (‘irtf’, ‘spex’, ‘kesseli+2017’, ‘bonnefoy+2014’).

• wavel_range (tuple(float, float), None) – Wavelength range (um) that is used for the empirical comparison.

• sptypes (list(str), None) – List with spectral types to compare with. The list should only contains types, for example sptypes=['M', 'L']. All available spectral types in the spec_library are compared with if set to None.

• av_ext (list(float), np.array, None) – List of $$A_V$$ for which the goodness-of-fit statistic is tested. The extinction is calculated with the empirical relation from Cardelli et al. (1989).

• rad_vel (list(float), np.array, None) – List of radial velocities (km s-1) for which the goodness-of-fit statistic is tested.

Returns

None

Return type

NoneType

## species.analysis.emission_line module#

Module with functionalities for the analysis of emission lines.

class species.analysis.emission_line.EmissionLine(object_name: str, spec_name: str, hydrogen_line: Optional[str] = None, lambda_rest: Optional[float] = None, wavel_range: Optional[Tuple[float, float]] = None)[source]#

Bases: object

Class for the analysis of emission lines.

Parameters
• object_name (str) – Object name as stored in the database with add_object() or add_companion().

• spec_name (str) – Name of the spectrum that is stored at the object data of object_name.

• hydrogen_line (str, None) – Name of the hydrogen line that will be analyzed. The names available lines can be checked with the  list_hydrogen_lines() method. If the argument is set to None then provide the rest wavelength as argument of lambda_rest.

• lambda_rest (float, None) – Rest wavelength (um) of the emission line. The parameter if used for calculating the radial velocity and its uncertainty. The argument can be set to None and will be ignored if the argument of hydrogen_line is used.

• wavel_range (tuple(float, float), None) – Wavelength range (um) that is cropped from the spectrum. The full spectrum is used if the argument is set to None.

Returns

None

Return type

NoneType

accretion_luminosity(line_lum: Union[float, numpy.ndarray]) Union[float, numpy.ndarray][source]#

Method for calculating the accretion luminosity from the (hydrogen) line luminosity with the relation from Aoyama et al. (2021).

Parameters

line_lum (float, np.array) – Line luminosity ($L_odot$) or array with line luminosities.

Returns

Accretion luminosity ($L_odot$) or array with accretion luminosities.

Return type

float, np.ndarray

fit_gaussian(tag: str, min_num_live_points: float = 400, bounds: Dict[str, Tuple[float, float]] = None, output: str = 'ultranest/', plot_filename: Optional[str] = 'line_fit.pdf', show_status: bool = True, double_gaussian: bool = False) None[source]#

Method for fitting a Gaussian profile to an emission line and using UltraNest for sampling the posterior distributions and estimating the evidence.

Parameters
• tag (str) – Database tag where the posterior samples will be stored.

• min_num_live_points (int) – Minimum number of live points (see https://johannesbuchner.github.io/UltraNest/issues.html).

• bounds (dict(str, tuple(float, float)), None) – The boundaries that are used for the uniform priors of the 3 Gaussian parameters (gauss_amplitude, gauss_mean, and gauss_sigma). Conservative prior boundaries will be estimated from the spectrum if the argument is set to None or if any of the required parameters is missing in the bounds dictionary.

• output (str) – Path that is used for the output files from UltraNest.

• plot_filename (str) – Filename for the plot with the best-fit line profile. The plot is shown in an interface window if the argument is set to None.

• show_status (bool) – Print information about the convergence.

• double_gaussian (bool) – Set to True for fitting a double instead of a single Gaussian. In that case, the bounds dictionary may also contain 'gauss_amplitude_2', 'gauss_mean_2', and 'gauss_sigma_2' (otherwise conservative parameter boundaries are estimated from the data).

Returns

None

Return type

NoneType

integrate_flux(wavel_int: Tuple[float, float], interp_kind: str = 'linear', plot_filename: Optional[str] = 'int_line.pdf') numpy.float64[source]#

Method for calculating the integrated line flux and error. The spectrum is first interpolated to $$R = 100000$$ and then integrated across the specified wavelength range with the composite trapezoidal rule of np.trapz. The error is estimated with a Monte Carlo approach from 1000 samples. The accretion luminosity is also calculated with the relation from Aoyama et al. (2021) if the argument of hydrogen_line was set.

Parameters
• wavel_int (tuple(float, float)) – Wavelength range (um) across which the flux will be integrated.

• interp_kind (str) – Kind of interpolation kind for scipy.interpolate.interp1d (default: ‘linear’).

• plot_filename (str, None) – Filename for the plot with the interpolated line profile. The plot is shown in an interface window if the argument is set to None.

Returns

• float – Integrated line flux (W m-2).

• float – Flux error (W m-2).

list_hydrogen_lines() List[str][source]#

Function to list the available hydrogen lines for which the accretion luminosity relation is available. These names can be set as argument hydrogen_line. In that case, the measured line luminosity from fit_gaussian() will be automatically converted to an accretion luminosity with the relation from Aoyama et al. (2021).

Returns

List with the names of the hydrogen lines for which there are coefficients available for the accretion relation.

Return type

list(str)

subtract_continuum(poly_degree: int = 3, plot_filename: Optional[str] = 'continuum.pdf', spec_filename: Optional[str] = None) None[source]#

Method for fitting the continuum with a polynomial function of the following form: $$P = \sum_{i=0}^{i=n}C_{i} * x^{i}$$. The spectrum is first smoothed with a median filter and then fitted with a linear least squares algorithm.

Parameters
• poly_degree (int) – Degree of the polynomial series.

• plot_filename (str, None) – Filename for the plots with the continuum fit and the continuum-subtracted spectrum. The plot is shown in an interface window if the argument is set to None.

• spec_filename (str, None) – Output text file for writing the continuum-subtracted spectrum. The file will not be created if the argument is set to None.

Returns

None

Return type

NoneType

## species.analysis.fit_model module#

Module with functionalities for fitting atmospheric model spectra.

class species.analysis.fit_model.FitModel(object_name: str, model: str, bounds: Optional[Dict[str, Union[Tuple[float, float], Tuple[Optional[Tuple[float, float]], Optional[Tuple[float, float]]], List[Tuple[float, float]]]]] = None, inc_phot: Union[bool, List[str]] = True, inc_spec: Union[bool, List[str]] = True, fit_corr: Optional[List[str]] = None, weights: Optional[Dict[str, float]] = None)[source]#

Bases: object

Class for fitting atmospheric model spectra to spectra and/or photometric fluxes, and using Bayesian inference (with MultiNest or UltraNest) to estimate the posterior distribution and marginalized likelihood (i.e. “evidence”). A grid of model spectra is linearly interpolated for each spectrum and photometric flux, while taking into account the filter profile, spectral resolution, and wavelength sampling. The computation time depends mostly on the number of free parameters and the resolution / number of data points of the spectra.

Parameters
• object_name (str) – Object name of the companion as stored in the database with add_object() or add_companion().

• model (str) – Name of the atmospheric model (e.g. ‘bt-settl’, ‘exo-rem’, ‘planck’, or ‘powerlaw’).

• bounds (dict(str, tuple(float, float)), None) –

The boundaries that are used for the uniform or log-uniform priors. Fixing a parameter is possible by providing the same value as lower and upper boundary of the parameter (e.g. bounds={'logg': (4., 4.). An explanation of the various parameters can be found below.

Atmospheric model parameters (e.g. with model='bt-settl-cifist'; see docstring of add_model() for the available model grids):

• Boundaries are provided as tuple of two floats. For example, bounds={'teff': (1000, 1500.), 'logg': (3.5, 5.)}.

• The grid boundaries (i.e. the maximum range) are adopted as priors if a parameter range is set to None or if a mandatory parameter is not included in the dictionary of bounds. For example, bounds={'teff': (1000., 1500.), 'logg': None}. The default range for the radius is $$0.5-5.0~R_\mathrm{J}$$. With bounds=None, automatic priors will be set for all mandatory parameters.

• It is possible to fit a weighted combination of two atmospheric parameters from the same model. This can be useful to fit data of a spectroscopic binary or to account for atmospheric asymmetries of a single object. For each atmospheric parameter, a tuple of two tuples can be provided, for example bounds={'teff': ((1000., 1500.), (1300., 1800.))}. Mandatory parameters that are not included are assumed to be the same for both components. The grid boundaries are used as parameter range if a component is set to None. For example, bounds={'teff': (None, None), 'logg': (4.0, 4.0), (4.5, 4.5)} will use the full range for $$T_\mathrm{eff}$$ of both components and fixes $$\log{g}$$ to 4.0 and 4.5, respectively. The spec_weight parameter is automatically included in the fit, as it sets the weight of the two components.

Blackbody parameters (with model='planck'):

• Parameter boundaries have to be provided for ‘teff’ and ‘radius’.

• For a single blackbody component, the values are provided as a tuple with two floats. For example, bounds={'teff': (1000., 2000.), 'radius': (0.8, 1.2)}.

• For multiple blackbody components, the values are provided as a list with tuples. For example, bounds={'teff': [(1000., 1400.), (1200., 1600.)], 'radius': [(0.8, 1.5), (1.2, 2.)]}.

• When fitting multiple blackbody components, an additional prior is used for restricting the temperatures and radii to decreasing and increasing values, respectively, in the order as provided in bounds.

Power-law spectrum (model='powerlaw'):

• Parameter boundaries have to be provided for ‘log_powerlaw_a’, ‘log_powerlaw_b’, and ‘log_powerlaw_c’. For example, bounds={'log_powerlaw_a': (-20., 0.), 'log_powerlaw_b': (-20., 5.), 'log_powerlaw_c': (-20., 5.)}.

• The spectrum is parametrized as $$\log10{f} = a + b*\log10{\lambda}^c$$, where $$a$$ is log_powerlaw_a, $$b$$ is log_powerlaw_b, and $$c$$ is log_powerlaw_c.

• Only implemented for fitting photometric fluxes, for example the IR fluxes of a star with disk. In that way, synthetic photometry can be calculated afterwards for a different filter. Note that this option assumes that the photometric fluxes are dominated by continuum emission while spectral lines are ignored.

• The plot_mag_posterior() function can be used for calculating synthetic photometry and error bars from the posterior distributions.

Calibration parameters:

• For each spectrum/instrument, two optional parameters can be fitted to account for biases in the calibration: a scaling of the flux and a relative inflation of the uncertainties.

• For example, bounds={'SPHERE': ((0.8, 1.2), (0., 1.))} if the scaling is fitted between 0.8 and 1.2, and the error is inflated (relative to the sampled model fluxes) with a value between 0 and 1.

• The dictionary key should be the same as the database tag of the spectrum. For example, {'SPHERE': ((0.8, 1.2), (0., 1.))} if the spectrum is stored as 'SPHERE' with add_object().

• Each of the two calibration parameters can be set to None in which case the parameter is not used. For example, bounds={'SPHERE': ((0.8, 1.2), None)}.

• The errors of the photometric fluxes can be inflated to account for underestimated error bars. The error inflation is relative to the actual flux and is either fitted separately for a filter, or a single error inflation is applied to all filters from an instrument. For the first case, the keyword in the bounds dictionary should be provided in the following format: 'Paranal/NACO.Mp_error': (0., 1.). Here, the error of the NACO $$M'$$ flux is inflated up to 100 percent of the actual flux. For the second case, only the telescope/instrument part of the the filter name should be provided in the bounds dictionary, so in the following format: 'Paranal/NACO_error': (0., 1.). This will increase the errors of all NACO filters by the same (relative) amount.

• No calibration parameters are fitted if the spectrum name is not included in bounds.

ISM extinction parameters:

• There are three approaches for fitting extinction. The first is with the empirical relation from Cardelli et al. (1989) for ISM extinction.

• The extinction is parametrized by the $V$ band extinction, $A_V$ (ism_ext), and optionally the reddening, R_V (ism_red). If ism_red is not provided, its value is fixed to 3.1 and not fitted.

• The prior boundaries of ism_ext and ism_red should be provided in the bounds dictionary, for example bounds={'ism_ext': (0., 10.), 'ism_red': (0., 20.)}.

Log-normal size distribution:

• The second approach is fitting the extinction of a log-normal size distribution of grains with a crystalline MgSiO3 composition, and a homogeneous, spherical structure.

• The size distribution is parameterized with a mean geometric radius (lognorm_radius in um) and a geometric standard deviation (lognorm_sigma, dimensionless). The grid of cross sections has been calculated for mean geometric radii between 0.001 and 10 um, and geometric standard deviations between 1.1 and 10.

• The extinction (lognorm_ext) is fitted in the $V$ band ($A_V$ in mag) and the wavelength-dependent extinction cross sections are interpolated from a pre-tabulated grid.

• The prior boundaries of lognorm_radius, lognorm_sigma, and lognorm_ext should be provided in the bounds dictionary, for example bounds={'lognorm_radius': (0.001, 10.), 'lognorm_sigma': (1.1, 10.), 'lognorm_ext': (0., 5.)}.

• A uniform prior is used for lognorm_sigma and lognorm_ext, and a log-uniform prior for lognorm_radius.

Power-law size distribution:

• The third approach is fitting the extinction of a power-law size distribution of grains, again with a crystalline MgSiO3 composition, and a homogeneous, spherical structure.

• The size distribution is parameterized with a maximum radius (powerlaw_max in um) and a power-law exponent (powerlaw_exp, dimensionless). The minimum radius is fixed to 1 nm. The grid of cross sections has been calculated for maximum radii between 0.01 and 100 um, and power-law exponents between -10 and 10.

• The extinction (powerlaw_ext) is fitted in the $V$ band ($A_V$ in mag) and the wavelength-dependent extinction cross sections are interpolated from a pre-tabulated grid.

• The prior boundaries of powerlaw_max, powerlaw_exp, and powerlaw_ext should be provided in the bounds dictionary, for example {'powerlaw_max': (0.01, 100.), 'powerlaw_exp': (-10., 10.), 'powerlaw_ext': (0., 5.)}.

• A uniform prior is used for powerlaw_exp and powerlaw_ext, and a log-uniform prior for powerlaw_max.

Blackbody disk emission:

• Additional blackbody emission can be added to the atmospheric spectrum to account for thermal emission from a disk.

• Parameter boundaries have to be provided for ‘disk_teff’ and ‘disk_radius’. For example, bounds={'teff': (2000., 3000.), 'radius': (1., 5.), 'logg': (3.5, 4.5), 'disk_teff': (100., 2000.), 'disk_radius': (1., 100.)}.

• inc_phot (bool, list(str)) – Include photometric data in the fit. If a boolean, either all (True) or none (False) of the data are selected. If a list, a subset of filter names (as stored in the database) can be provided.

• inc_spec (bool, list(str)) – Include spectroscopic data in the fit. If a boolean, either all (True) or none (False) of the data are selected. If a list, a subset of spectrum names (as stored in the database with add_object()) can be provided.

• fit_corr (list(str), None) – List with spectrum names for which the covariances are modeled with a Gaussian process (see Wang et al. 2020). This option can be used if the actual covariances as determined from the data are not available for the spectra of object_name. The parameters that will be fitted are the correlation length and the fractional amplitude.

• weights (dict(str, float), None) – Weights to be applied to the log-likelihood components of the different spectroscopic and photometric data that are provided with inc_spec and inc_phot. This parameter can for example be used to bias the weighting of the photometric data points. An equal weighting is applied if the argument is set to None.

Returns

None

Return type

NoneType

lnlike_func(params, prior: Optional[Dict[str, Tuple[float, float]]]) numpy.float64[source]#

Function for calculating the log-likelihood for the sampled parameter cube.

Parameters
• params (np.ndarray, pymultinest.run.LP_c_double) – Cube with physical parameters.

• prior (dict(str, tuple(float, float))) – Dictionary with Gaussian priors for one or multiple parameters. The prior can be set for any of the atmosphere or calibration parameters, e.g. prior={'teff': (1200., 100.)}. Additionally, a prior can be set for the mass, e.g. prior={'mass': (13., 3.)} for an expected mass of 13 Mjup with an uncertainty of 3 Mjup.

Returns

Log-likelihood.

Return type

float

run_multinest(tag: str, n_live_points: int = 1000, output: str = 'multinest/', prior: Optional[Dict[str, Tuple[float, float]]] = None) None[source]#

Function to run the PyMultiNest wrapper of the MultiNest sampler. While PyMultiNest can be installed with pip from the PyPI repository, MultiNest has to to be build manually. See the PyMultiNest documentation for details: http://johannesbuchner.github.io/PyMultiNest/install.html. Note that the library path of MultiNest should be set to the environmental variable LD_LIBRARY_PATH on a Linux machine and DYLD_LIBRARY_PATH on a Mac. Alternatively, the variable can be set before importing the species package, for example:

>>> import os
>>> os.environ['DYLD_LIBRARY_PATH'] = '/path/to/MultiNest/lib'
>>> import species

Parameters
• tag (str) – Database tag where the samples will be stored.

• n_live_points (int) – Number of live points.

• output (str) – Path that is used for the output files from MultiNest.

• prior (dict(str, tuple(float, float)), None) – Dictionary with Gaussian priors for one or multiple parameters. The prior can be set for any of the atmosphere or calibration parameters, for example prior={'teff': (1200., 100.)}. Additionally, a prior can be set for the mass, for example prior={'mass': (13., 3.)} for an expected mass of 13 Mjup with an uncertainty of 3 Mjup. The parameter is not used if set to None.

Returns

None

Return type

NoneType

run_ultranest(tag: str, min_num_live_points=400, output: str = 'ultranest/', prior: Optional[Dict[str, Tuple[float, float]]] = None) None[source]#

Function to run UltraNest for constructing the posterior probability distributions on model parameters and computing the marginal likelihood (i.e. “evidence”).

Parameters
• tag (str) – Database tag where the samples will be stored.

• min_num_live_points (int) – Minimum number of live points. The default of 400 is a reasonable number (see https://johannesbuchner.github.io/UltraNest/issues.html). In principle, choosing a very low number allows nested sampling to make very few iterations and go to the peak quickly. However, the space will be poorly sampled, giving a large region and thus low efficiency, and potentially not seeing interesting modes. Therefore, a value above 100 is typically useful.

• output (str) – Path that is used for the output files from UltraNest.

• prior (dict(str, tuple(float, float)), None) – Dictionary with Gaussian priors for one or multiple parameters. The prior can be set for any of the atmosphere or calibration parameters, e.g. prior={'teff': (1200., 100.)}. Additionally, a prior can be set for the mass, e.g. prior={'mass': (13., 3.)} for an expected mass of 13 Mjup with an uncertainty of 3 Mjup. The parameter is not used if set to None.

Returns

None

Return type

NoneType

## species.analysis.fit_spectrum module#

Module with functionalities for photometric and spectroscopic calibration. The fitting routine can be used to fit photometric data with a calibration spectrum (e.g. extracted with get_model()) by simply fitting a scaling parameter.

class species.analysis.fit_spectrum.FitSpectrum(object_name: str, filters: Optional[List[str]], spectrum: str, bounds: Dict[str, Tuple[float, float]])[source]#

Bases: object

Class for fitting a calibration spectrum to photometric data.

Parameters
• object_name (str) – Object name in the database.

• filters (list(str)) – Filter names for which the photometry is selected. All available photometry of the object is selected if set to None.

• spectrum (str) – Calibration spectrum as labelled in the database. The calibration spectrum can be stored in the database with add_calibration().

• bounds (dict) – Boundaries of the scaling parameter, as {'scaling':(min, max)}.

Returns

None

Return type

NoneType

run_mcmc(nwalkers: int, nsteps: int, guess: Union[Dict[str, float], Dict[str, None]], tag: str) None[source]#

Function to run the MCMC sampler.

Parameters
• nwalkers (int) – Number of walkers.

• nsteps (int) – Number of steps per walker.

• guess (dict(str, float), dict(str, None)) – Guess of the scaling parameter.

• tag (str) – Database tag where the MCMC samples will be stored.

Returns

None

Return type

NoneType

species.analysis.fit_spectrum.lnprob(param: numpy.ndarray, bounds: Dict[str, Tuple[float, float]], modelpar: List[str], objphot: List[numpy.ndarray], specphot: Union[List[float], List[Tuple[species.analysis.photometry.SyntheticPhotometry, Tuple[numpy.float64, numpy.float64]]]]) float[source]#

Internal function for calculating the posterior probability.

Parameters
• param (np.ndarray) – Value of the scaling parameter.

• bounds (dict) – Boundaries of the main scaling parameter.

• modelpar (list(str)) – Parameter names.

• objphot (list(tuple(float, float))) – Photometry of the object.

• specphot (list(float), photometry.SyntheticPhotometry) – Synthetic photometry of the calibration spectrum for the same filters as the photometry of the object.

Returns

Log posterior probability.

Return type

float

## species.analysis.photometry module#

Module with functionalities for calculating synthetic photometry.

class species.analysis.photometry.SyntheticPhotometry(filter_name: str)[source]#

Bases: object

Class for calculating synthetic photometry from a spectrum and also for conversion between magnitudes and fluxes. Note that depending on the detector type (energy- or photon-counting) the integral for the filter-weighted flux contains an additional wavelength factor.

Parameters

filter_name (str) – Filter name as listed in the database. Filters from the SVO Filter Profile Service are automatically downloaded and added to the database.

Returns

None

Return type

NoneType

flux_to_magnitude(flux: float, error: Optional[Union[float, numpy.ndarray]] = None, parallax: Optional[Union[Tuple[float, Optional[float]], Tuple[numpy.ndarray, Optional[numpy.ndarray]]]] = None, distance: Optional[Union[Tuple[float, Optional[float]], Tuple[numpy.ndarray, Optional[numpy.ndarray]]]] = None) Tuple[Union[Tuple[float, Optional[float]], Tuple[numpy.ndarray, Optional[numpy.ndarray]]], Union[Tuple[Optional[float], Optional[float]], Tuple[Optional[numpy.ndarray], Optional[numpy.ndarray]]]][source]#

Function for converting a flux into a magnitude.

Parameters
• flux (float, np.ndarray) – Flux (W m-2 um-1).

• error (float, np.ndarray, None) – Uncertainty (W m-2 um-1). Not used if set to None.

• parallax (tuple(float, float), , tuple(np.ndarray, np.ndarray), None) – Parallax and uncertainty (mas). The returned absolute magnitude is set to None in case parallax and distance are set to None. The error is not propagated into the error on the absolute magnitude in case the parallax uncertainty is set to None, for example parallax=(10., None).

• distance (tuple(float, float), tuple(np.ndarray, np.ndarray), None) – Distance and uncertainty (pc). The returned absolute magnitude is set to None in case distance and parallax are set to None. The error is not propagated into the error on the absolute magnitude in case the distance uncertainty is set to None, for example distance=(20., None). This parameter is ignored if the parallax parameter is used.

Returns

• tuple(float, float), tuple(np.ndarray, np.ndarray) – Apparent magnitude and uncertainty.

• tuple(float, float), tuple(np.ndarray, np.ndarray) – Absolute magnitude and uncertainty.

magnitude_to_flux(magnitude: float, error: Optional[float] = None, zp_flux: Optional[float] = None) Tuple[numpy.float64, numpy.float64][source]#

Function for converting a magnitude to a flux.

Parameters
• magnitude (float) – Magnitude.

• error (float, None) – Error on the magnitude. Not used if set to None.

• zp_flux (float, None) – Zero-point flux (W m-2 um-1). The value is calculated if the argument is set to None.

Returns

• float – Flux (W m-2 um-1).

• float – Error (W m-2 um-1).

spectrum_to_flux(wavelength: numpy.ndarray, flux: numpy.ndarray, error: Optional[numpy.ndarray] = None, threshold: Optional[float] = 0.05) Tuple[Union[numpy.float32, numpy.float64], Union[numpy.float32, None, numpy.float64]][source]#

Function for calculating the average flux from a spectrum and a filter profile. The uncertainty is propagated by sampling 200 random values from the error distributions.

Parameters
• wavelength (np.ndarray) – Wavelength points (um).

• flux (np.ndarray) – Flux (W m-2 um-1).

• error (np.ndarray, None) – Uncertainty (W m-2 um-1). Not used if set to None.

• threshold (float, None) – Transmission threshold (value between 0 and 1). If the minimum transmission value is larger than the threshold, a NaN is returned. This will happen if the input spectrum does not cover the full wavelength range of the filter profile. The parameter is not used if set to None.

Returns

• float – Average flux (W m-2 um-1).

• float, None – Uncertainty (W m-2 um-1).

spectrum_to_magnitude(wavelength: numpy.ndarray, flux: numpy.ndarray, error: Optional[Union[numpy.ndarray, List[numpy.ndarray]]] = None, parallax: Optional[Tuple[float, Optional[float]]] = None, distance: Optional[Tuple[float, Optional[float]]] = None, threshold: Optional[float] = 0.05) Tuple[Tuple[float, Optional[float]], Optional[Tuple[Optional[float], Optional[float]]]][source]#

Function for calculating the apparent and absolute magnitude from a spectrum and a filter profile. The uncertainty is propagated by sampling 200 random values from the error distributions.

Parameters
• wavelength (np.ndarray) – Wavelength points (um).

• flux (np.ndarray) – Flux (W m-2 um-1).

• error (np.ndarray, list(np.ndarray), None) – Uncertainty (W m-2 um-1).

• parallax (tuple(float, float), None) – Parallax and uncertainty (mas). No absolute magnitude is calculated if set to None. No error on the absolute magnitude is calculated if the error parameter is set to None.

• distance (tuple(float, float), None) – Distance and uncertainty (pc). No absolute magnitude is calculated if set to None. No error on the absolute magnitude is calculated if the error parameter is set to None. This parameter is ignored if the parallax parameter is used.

• threshold (float, None) – Transmission threshold (value between 0 and 1). If the minimum transmission value is larger than the threshold, a NaN is returned. This will happen if the input spectrum does not cover the full wavelength range of the filter profile. The parameter is not used if set to None.

Returns

• tuple(float, float) – Apparent magnitude and uncertainty.

• tuple(float, float) – Absolute magnitude and uncertainty.

zero_point() numpy.float64[source]#

Internal function for calculating the zero point of the provided filter_name.

Returns

Zero-point flux (W m-2 um-1).

Return type

float

## species.analysis.retrieval module#

Module with a frontend for atmospheric retrieval with the radiative transfer and retrieval code petitRADTRANS (see https://petitradtrans.readthedocs.io).

class species.analysis.retrieval.AtmosphericRetrieval(object_name: str, line_species: Optional[List[str]] = None, cloud_species: Optional[List[str]] = None, output_folder: str = 'multinest', wavel_range: Optional[Tuple[float, float]] = None, scattering: bool = True, inc_spec: Union[bool, List[str]] = True, inc_phot: Union[bool, List[str]] = False, pressure_grid: str = 'smaller', weights: Optional[Dict[str, float]] = None, lbl_species: Optional[List[str]] = None, max_pressure: float = 1000.0)[source]#

Bases: object

Class for atmospheric retrievals of self-luminous atmospheres of giant planets and brown dwarfs within a Bayesian framework. This class provides a frontend for petitRADTRANS, with a variety of P-T profiles, cloud models, priors, and more.

Parameters
• object_name (str) – Name of the object as stored in the database with add_object().

• line_species (list, None) – List with the line species. A minimum of one line species should be included.

• cloud_species (list, None) – List with the cloud species. No cloud species are used if the argument is to None.

• output_folder (str) – Folder name that is used for the output files from MultiNest. The folder is created if it does not exist.

• wavel_range (tuple(float, float), None) – The wavelength range (um) that is used for the forward model. Should be a bit broader than the minimum and maximum wavelength of the data. If photometric fluxes are included (see inc_phot), it is important that wavel_range encompasses the full filter profile, which can be inspected with the functionalities of ReadFilter. The wavelength range is set automatically if the argument is set to None.

• scattering (bool) – Turn on scattering in the radiative transfer. Only recommended at infrared wavelengths when clouds are included in the forward model. Using scattering will increase the computation time significantly.

• inc_spec (bool, list(str)) – Include spectroscopic data in the fit. If a boolean, either all (True) or none (False) of the available data are selected. If a list, a subset of spectrum names (as stored in the database with add_object()) can be provided.

• inc_phot (bool, list(str)) – Include photometric data in the fit. If a boolean, either all (True) or none (False) of the available data are selected. If a list, a subset of filter names (as stored in the database with add_object()) can be provided.

• pressure_grid (str) – The type of pressure grid that is used for the radiative transfer. Either ‘standard’, to use 180 layers both for the atmospheric structure (e.g. when interpolating the abundances) and 180 layers with the radiative transfer, or ‘smaller’ to use 60 (instead of 180) with the radiative transfer, or ‘clouds’ to start with 1440 layers but resample to ~100 layers (depending on the number of cloud species) with a refinement around the cloud decks. For cloudless atmospheres it is recommended to use ‘smaller’, which runs faster than ‘standard’ and provides sufficient accuracy. For cloudy atmosphere, it is recommended to test with ‘smaller’ but it might be required to use ‘clouds’ to improve the accuracy of the retrieved parameters, at the cost of a long runtime.

• weights (dict(str, float), None) – Weights to be applied to the log-likelihood components of the different spectroscopic and photometric data that are provided with inc_spec and inc_phot. This parameter can for example be used to increase the weighting of the photometric data points relative to the spectroscopic data. An equal weighting is applied if the argument is set to None.

• lbl_species (list, None) – List with the line species that will be used for calculating line-by-line spectra for the list of high-resolution spectra that are provided as argument of cross_corr when starting the retrieval with species.analysis.retrieval.AtmosphericRetrieval.run_multinest(). The argument can be set to None when cross_corr=None. The lbl_species and cross_corr parameters should only be used if the log-likelihood component should be determined with a cross-correlation instead of a direct comparison of data and model.

• max_pressure (float) – Maximum pressure (bar) that is used for the P-T profile. The default is set to 1000 bar.

Returns

None

Return type

NoneType

rebin_opacities(wavel_bin: float, out_folder: str = 'rebin_out') None[source]#

Function for downsampling the c-k opacities from $$\lambda/\Delta\lambda = 1000$$ to a smaller wavelength binning. The downsampled opacities should be stored in the opacities/lines/corr_k/ folder of pRT_input_data_path.

Parameters
• wavel_bin (float) – Wavelength binning, $$\lambda/\Delta\lambda$$, to which the opacities will be downsampled.

• out_folder (str) – Path of the output folder where the downsampled opacities will be stored.

Returns

None

Return type

NoneType

run_multinest(bounds: dict, chemistry: str = 'equilibrium', quenching: Optional[str] = 'pressure', pt_profile: str = 'molliere', fit_corr: Optional[List[str]] = None, cross_corr: Optional[List[str]] = None, n_live_points: int = 2000, resume: bool = False, plotting: bool = False, check_isothermal: bool = False, pt_smooth: Optional[float] = 0.3, check_flux: Optional[float] = None, temp_nodes: Optional[int] = None, prior: Optional[Dict[str, Tuple[float, float]]] = None, check_phot_press: Optional[float] = None) None[source]#

Function for running the atmospheric retrieval. The parameter estimation and computation of the marginalized likelihood (i.e. model evidence), is done with PyMultiNest wrapper of the MultiNest sampler. While PyMultiNest can be installed with pip from the PyPI repository, MultiNest has to to be compiled manually. See the PyMultiNest documentation: http://johannesbuchner.github.io/PyMultiNest/install.html. Note that the library path of MultiNest should be set to the environment variable LD_LIBRARY_PATH on a Linux machine and DYLD_LIBRARY_PATH on a Mac. Alternatively, the variable can be set before importing the species toolkit, for example:

>>> import os
>>> os.environ['DYLD_LIBRARY_PATH'] = '/path/to/MultiNest/lib'
>>> import species


When using MPI, it is also required to install mpi4py (e.g. pip install mpi4py), otherwise an error may occur when the output_folder is created by multiple processes.

Parameters
• bounds (dict) –

The boundaries that are used for the uniform or log-uniform priors. The dictionary contains the parameters as key and the boundaries as value. The boundaries are provided as a tuple with two values (lower and upper boundary). Fixing a parameter is possible by providing the same value as lower and upper boundary of the parameter (e.g. bounds={'logg': (4., 4.). An explanation of the mandatory and optional parameters can be found in the description of the model_param parameter of species.read.read_radtrans.ReadRadtrans.get_model(). Additional parameters that can specifically be used for a retrieval are listed below.

Scaling parameters (mandatory):

• The radius ($$R_\mathrm{J}$$), radius, is a mandatory parameter to include. It is used for scaling the flux from the planet surface to the observer.

• The parallax (mas), parallax, is also used for scaling the flux. However, this parameter is automatically included in the retrieval with a Gaussian prior (based on the object data of object_name). So this parameter does not need to be included in bounds).

Calibration parameters (optional):

• For each spectrum/instrument, three optional parameters can be fitted to account for biases in the calibration: a scaling of the flux, a constant inflation of the uncertainties, and a constant offset in the wavelength solution.

• For example, bounds={'SPHERE': ((0.8, 1.2), (-16., -14.), (-0.01, 0.01))} if the scaling is fitted between 0.8 and 1.2, each uncertainty is inflated with a constant value between $$10^{-16}$$ and $$10^{-14}$$ W $$\mathrm{m}^{-2}$$ $$\mu\mathrm{m}^{-1}$$, and a constant wavelength offset between -0.01 and 0.01 $$\mu\mathrm{m}$$

• The dictionary key should be the same as to the database tag of the spectrum. For example, {'SPHERE': ((0.8, 1.2), (-16., -14.), (-0.01, 0.01))} if the spectrum is stored as 'SPHERE' with add_object().

• Each of the three calibration parameters can be set to None in which case the parameter is not used. For example, bounds={'SPHERE': ((0.8, 1.2), None, None)}.

• No calibration parameters are fitted if the spectrum name is not included in bounds.

Prior parameters (optional):

• The log_sigma_alpha parameter can be used when pt_profile='molliere'. This prior penalizes samples if the parametrized, pressure-dependent opacity is not consistent with the atmosphere’s non-gray opacity structure (see GRAVITY Collaboration et al. 2020 for details).

• The log_gamma_r and log_beta_r parameters can be included when pt_profile='monotonic' or pt_profile='free'. A prior will be applied that penalizes wiggles in the P-T profile through the second derivative of the temperature structure (see Line et al. (2015) for details).

• chemistry (str) – The chemistry type: ‘equilibrium’ for equilibrium chemistry or ‘free’ for retrieval of free abundances (but constant with altitude).

• quenching (str, None) – Quenching type for CO/CH4/H2O abundances. Either the quenching pressure (bar) is a free parameter (quenching='pressure') or the quenching pressure is calculated from the mixing and chemical timescales (quenching='diffusion'). The quenching is not applied if the argument is set to None.

• pt_profile (str) – The parametrization for the pressure-temperature profile (‘molliere’, ‘free’, ‘monotonic’, ‘eddington’).

• fit_corr (list(str), None) – List with spectrum names for which the correlation lengths and fractional amplitudes are fitted (see Wang et al. 2020) to model the covariances in case these are not available.

• cross_corr (list(str), None) – List with spectrum names for which a cross-correlation to log-likelihood mapping is used (see Brogi & Line 2019) instead of a direct comparison of model an data with a least-squares approach. This parameter should only be used for high-resolution spectra. Currently, it only supports spectra that have been shifted to the planet’s rest frame.

• n_live_points (int) – Number of live points used for the nested sampling.

• resume (bool) – Resume from a previous run.

• plotting (bool) – Plot sample results for testing purpose. Not recommended to use when running the full retrieval.

• check_isothermal (bool) – Check if there is an isothermal region below 1 bar. If so, discard the sample. This parameter is experimental and has not been properly implemented.

• pt_smooth (float, None) – Standard deviation of the Gaussian kernel that is used for smoothing the P-T profile, after the temperature nodes have been interpolated to a higher pressure resolution. Only required with pt_profile=’free’ or pt_profile=’monotonic’. The argument should be given as $$\log10{P/\mathrm{bar}}$$, with the default value set to 0.3 dex. No smoothing is applied if the argument if set to 0 or None. The pt_smooth parameter can also be included in bounds, in which case the value is fitted and the pt_smooth argument is ignored.

• check_flux (float, None) – Relative tolerance for enforcing a constant bolometric flux at all pressures layers. By default, only the radiative flux is used for the bolometric flux. The convective flux component is also included if the mix_length parameter (relative to the pressure scale height) is included in the bounds dictionary. To use check_flux, the opacities should be recreated with rebin_opacities() at $R = 10$ (i.e. spec_res=10) and placed in the folder of pRT_input_data_path. This parameter is experimental and has not been fully tested.

• temp_nodes (int, None) – Number of free temperature nodes that are used with pt_profile='monotonic' or pt_profile='free'.

• prior (dict(str, tuple(float, float)), None) – Dictionary with Gaussian priors for one or multiple parameters. The prior can be set for any of the atmosphere or calibration parameters, for example prior={'logg': (4.2, 0.1)}. Additionally, a prior can be set for the mass, for example prior={'mass': (13., 3.)} for an expected mass of 13 Mjup with an uncertainty of 3 Mjup. The parameter is not used if set to None.

• check_phot_press (float, None) – Remove the sample if the photospheric pressure that is calculated for the P-T profile is more than a factor check_phot_press larger or smaller than the photospheric pressure that is calculated from the Rosseland mean opacity of the non-gray opacities of the atmospheric structure (see Eq. 7 in GRAVITY Collaboration et al. 2020, where a factor of 5 was used). This parameter can only in combination with pt_profile='molliere'. The parameter is not used used if set to None. Finally, since samples are removed when not full-filling this requirement, the runtime of the retrieval may increase significantly.

Returns

None

Return type

NoneType

set_parameters(bounds: dict, chemistry: str, quenching: Optional[str], pt_profile: str, fit_corr: List[str], rt_object) None[source]#

Function to add the model parameters to the list of the parameters attribute.

Parameters
• bounds (dict) – Dictionary with the boundaries that are used as uniform priors for the parameters.

• chemistry (str) – The chemistry type: ‘equilibrium’ for equilibrium chemistry or ‘free’ for retrieval of free abundances (but constant with altitude).

• quenching (str, None) – Quenching type for CO/CH4/H2O abundances. Either the quenching pressure (bar) is a free parameter (quenching='pressure') or the quenching pressure is calculated from the mixing and chemical timescales (quenching='diffusion'). The quenching is not applied if the argument is set to None.

• pt_profile (str) – The parametrization for the pressure-temperature profile (‘molliere’, ‘free’, ‘monotonic’, ‘eddington’).

• fit_corr (list(str), None) –

List with spectrum names for which the correlation lengths and fractional amplitudes are fitted (see Wang et al. 2020) to model the covariances in case these are not available.

• rt_object (petitRADTRANS.radtrans.Radtrans) – Instance of Radtrans from petitRADTRANS.

Returns

None

Return type

NoneType