Added wavel_resample parameter in Database.get_mcmc_spectra and ReadP…

…lanck.get_spectrum, added phot_type parameter in Database.get_mcmc_photometry, minor improvements in docstrings

species/data/database.py

@@ -1376,27 +1376,34 @@ def get_mcmc_spectra(self,
                          random: int,
                          burnin: Optional[int] = None,
                          wavel_range: Optional[Union[Tuple[float, float], str]] = None,
-                         spec_res: Optional[float] = None) -> Union[List[box.ModelBox],
-                                                                    List[box.SpectrumBox]]:
+                         spec_res: Optional[float] = None,
+                         wavel_resample: Optional[np.ndarray] = None) -> \
+                             Union[List[box.ModelBox], List[box.SpectrumBox]]:
         """
+        Function for drawing random spectra from the sampled posterior distributions.
+
         Parameters
         ----------
         tag : str
             Database tag with the posterior samples.
         random : int
             Number of random samples.
         burnin : int, None
-            Number of burnin steps. No burnin is removed if set to ``None``.
+            Number of burnin steps. No burnin is removed if set to ``None``. Not required when
+            using nested sampling.
         wavel_range : tuple(float, float), str, None
             Wavelength range (um) or filter name. Full spectrum is used if set to ``None``.
         spec_res : float
             Spectral resolution that is used for the smoothing with a Gaussian kernel. No smoothing
-            is applied if set to ``None``.
+            is applied if the argument set to ``None``.
+        wavel_resample : np.ndarray, None
+            Wavelength points (um) to which the model spectrum will be resampled. The resampling is
+            applied after the optional smoothing to the resolution of ``spec_res``.
 
         Returns
         -------
         list(species.core.box.ModelBox)
-            Boxes with the randomly sampled spectra.
+            List with ``ModelBox`` objects.
         """
 
         if burnin is None:
@@ -1498,13 +1505,23 @@ def get_mcmc_spectra(self,
 
             if spectrum_type == 'model':
                 if spectrum_name == 'planck':
-                    specbox = readmodel.get_spectrum(model_param, spec_res)
+                    specbox = readmodel.get_spectrum(model_param,
+                                                     spec_res,
+                                                     smooth=True,
+                                                     wavel_resample=wavel_resample)
 
                 elif spectrum_name == 'powerlaw':
+                    if wavel_resample is not None:
+                        warnings.warn('The \'wavel_resample\' parameter is not support by the '
+                                      '\'powerlaw\' model so the argument will be ignored.')
+
                     specbox = read_util.powerlaw_spectrum(wavel_range, model_param)
 
                 else:
-                    specbox = readmodel.get_model(model_param, spec_res=spec_res, smooth=True)
+                    specbox = readmodel.get_model(model_param,
+                                                  spec_res=spec_res,
+                                                  wavel_resample=wavel_resample,
+                                                  smooth=True)
 
             elif spectrum_type == 'calibration':
                 specbox = readcalib.get_spectrum(model_param)
@@ -1519,23 +1536,33 @@ def get_mcmc_spectra(self,
     def get_mcmc_photometry(self,
                             tag: str,
                             filter_name: str,
-                            burnin: Optional[int] = None) -> np.ndarray:
+                            burnin: Optional[int] = None,
+                            phot_type: str = 'magnitude') -> np.ndarray:
         """
+        Function for calculating synthetic magnitudes or fluxes from the posterior samples.
+
         Parameters
         ----------
         tag : str
             Database tag with the posterior samples.
         filter_name : str
-            Filter name for which the photometry will be computed.
+            Filter name for which the synthetic photometry will be computed.
         burnin : int, None
-            Number of burnin steps. No burnin is removed if set to ``None``.
+            Number of burnin steps. No burnin is removed if set to ``None``. Not required when
+            using nested sampling.
+        phot_type : str
+            Photometry type ('magnitude' or 'flux').
 
         Returns
         -------
         np.ndarray
-            Synthetic magnitudes.
+            Synthetic magnitudes or fluxes (W m-2 um-1).
         """
 
+        if phot_type not in ['magnitude', 'flux']:
+            raise ValueError('The argument of \'phot_type\' is not recognized and should be '
+                             'set to \'magnitude\' or \'flux\'.')
+
         if burnin is None:
             burnin = 0
 
@@ -1586,6 +1613,7 @@ def get_mcmc_photometry(self,
 
         for i in tqdm.tqdm(range(samples.shape[0]), desc='Getting MCMC photometry'):
             model_param = {}
+
             for j in range(n_param):
                 model_param[param[j]] = samples[i, j]
 
@@ -1596,16 +1624,27 @@ def get_mcmc_photometry(self,
                 if spectrum_name == 'powerlaw':
                     pl_box = read_util.powerlaw_spectrum(synphot.wavel_range, model_param)
 
-                    app_mag, _ = synphot.spectrum_to_magnitude(pl_box.wavelength, pl_box.flux)
+                    if phot_type == 'magnitude':
+                        app_mag, _ = synphot.spectrum_to_magnitude(pl_box.wavelength, pl_box.flux)
+                        mcmc_phot[i] = app_mag[0]
 
-                    mcmc_phot[i] = app_mag[0]
+                    elif phot_type == 'flux':
+                        mcmc_phot[i], _ = synphot.spectrum_to_flux(pl_box.wavelength, pl_box.flux)
 
                 else:
-                    mcmc_phot[i], _ = readmodel.get_magnitude(model_param)
+                    if phot_type == 'magnitude':
+                        mcmc_phot[i], _ = readmodel.get_magnitude(model_param)
+
+                    elif phot_type == 'flux':
+                        mcmc_phot[i], _ = readmodel.get_flux(model_param)
 
             elif spectrum_type == 'calibration':
-                app_mag, _ = readcalib.get_magnitude(model_param=model_param, distance=None)
-                mcmc_phot[i] = app_mag[0]
+                if phot_type == 'magnitude':
+                    app_mag, _ = readcalib.get_magnitude(model_param=model_param, distance=None)
+                    mcmc_phot[i] = app_mag[0]
+
+                elif phot_type == 'flux':
+                    mcmc_phot[i], _ = readcalib.get_flux(model_param=model_param)
 
         return mcmc_phot
 

species/read/read_calibration.py

@@ -263,8 +263,12 @@ def get_flux(self,
 
         Returns
         -------
-        tuple(float, float)
-            Average flux and uncertainty (W m-2 um-1).
+        Returns
+        -------
+        float
+            Average flux (W m-2 um-1).
+        float
+            Uncertainty (W m-2 um-1).
         """
 
         specbox = self.get_spectrum(model_param=model_param)

species/read/read_planck.py

@@ -2,13 +2,14 @@
 Module with reading functionalities for Planck spectra.
 """
 
-import os
-import math
 import configparser
+import math
+import os
 
-from typing import Optional, Union, Dict, Tuple, List
+from typing import Dict, List, Optional, Tuple, Union
 
 import numpy as np
+import spectres
 
 from typeguard import typechecked
 
@@ -129,9 +130,12 @@ def update_parameters(model_param: Dict[str, Union[float, List[float]]]) -> Dict
     def get_spectrum(self,
                      model_param: Dict[str, Union[float, List[float]]],
                      spec_res: float,
-                     smooth: bool = False) -> box.ModelBox:
+                     smooth: bool = False,
+                     wavel_resample: Optional[np.ndarray] = None) -> box.ModelBox:
         """
-        Function for calculating a Planck spectrum or a combination of multiple Planck spectra.
+        Function for calculating a Planck spectrum or a combination of multiple Planck spectra. The
+        spectrum is calculated at :math:`R = 500`. Afterwards, an optional smoothing and wavelength
+        resampling can be applied.
 
         Parameters
         ----------
@@ -141,8 +145,13 @@ def get_spectrum(self,
             of multiple Planck functions, e.g.
             {'teff': [1500., 1000.], 'radius': [1., 2.], 'distance': 10.}.
         spec_res : float
-            Spectral resolution.
+            Spectral resolution that is used for smoothing the spectrum with a Gaussian kernel when
+            ``smooth=True``.
         smooth : bool
+            If ``True``, the spectrum is smoothed to the spectral resolution of ``spec_res``.
+        wavel_resample : np.ndarray, None
+            Wavelength points (um) to which the spectrum will be resampled. The resampling is
+            applied after the optional smoothing to ``spec_res`` when ``smooth=True``.
 
         Returns
         -------
@@ -153,7 +162,7 @@ def get_spectrum(self,
         if 'teff' in model_param and isinstance(model_param['teff'], list):
             model_param = self.update_parameters(model_param)
 
-        wavel_points = read_util.create_wavelengths(self.wavel_range, spec_res)
+        wavel_points = read_util.create_wavelengths(self.wavel_range, 500.)
 
         n_planck = 0
 
@@ -198,6 +207,17 @@ def get_spectrum(self,
                                    parameters=model_param,
                                    quantity='flux')
 
+        if wavel_resample is not None:
+            flux = spectres.spectres(wavel_resample,
+                                     wavel_points,
+                                     flux,
+                                     spec_errs=None,
+                                     fill=np.nan,
+                                     verbose=True)
+
+            model_box.wavelength = wavel_resample
+            model_box.flux = flux
+
         if n_planck == 1 and 'radius' in model_param:
             model_box.parameters['luminosity'] = 4. * np.pi * (
                 model_box.parameters['radius'] * constants.R_JUP)**2 * constants.SIGMA_SB * \