Support for fitting multiple disk_teff and disk_radius parameters

species/data/spec_data/spec_vega.py

@@ -49,7 +49,7 @@ def add_vega(input_path, database):
                 progressbar=True,
             )
 
-        except requests.exceptions.HTTPError:
+        except (requests.exceptions.HTTPError, requests.exceptions.ConnectionError):
             url = (
                 "https://home.strw.leidenuniv.nl/~stolker/"
                 "species/alpha_lyr_stis_011.fits"

species/fit/fit_model.py

@@ -188,9 +188,34 @@ def __init__(
                  which adds a constant flux (in W m-2 um-1) to the
                  model spectrum.
 
-            Blackbody parameters (with ``model='planck'``):
+            Blackbody disk emission:
+
+               - Blackbody parameters can be fitted to account for
+                 thermal emission from one or multiple disk
+                 components, in addition to the atmospheric
+                 emission. These parameters should therefore be
+                 combined with an atmospheric model.
 
-               - Parameter boundaries have to be provided for 'teff'
+               - 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.)}`` for fitting a single
+                 blackbody component, in addition to the atmospheric
+                 parameters. Or, ``bounds={'teff': (2000., 3000.),
+                 'radius': (1., 5.), 'logg': (3.5, 4.5),
+                 'disk_teff': [(2000., 500.), (1000., 20.)],
+                 'disk_radius': [(1., 100.), (50., 1000.)]}`` for
+                 fitting two blackbody components. Any number of
+                 blackbody components can be fitted by including
+                 additional priors in the lists of ``'disk_teff'``
+                 and ``'disk_radius'``.
+
+            Blackbody parameters (only with ``model='planck'``):
+
+               - This implementation fits both the atmospheric emission
+                 and possible disk emission with blackbody components.
+                 Parameter boundaries have to be provided for 'teff'
                  and 'radius'.
 
                - For a single blackbody component, the values are
@@ -363,18 +388,6 @@ def __init__(
                    ``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
@@ -649,8 +662,32 @@ def __init__(
             self.n_planck = 0
 
             if "disk_teff" in self.bounds and "disk_radius" in self.bounds:
-                self.modelpar.append("disk_teff")
-                self.modelpar.append("disk_radius")
+                if isinstance(bounds["disk_teff"], list) and isinstance(
+                    bounds["disk_radius"], list
+                ):
+                    # Update temperature and radius parameters
+                    # in case of multiple disk components
+                    self.n_disk = len(bounds["disk_teff"])
+
+                    for i, item in enumerate(bounds["disk_teff"]):
+                        self.modelpar.append(f"disk_teff_{i}")
+                        self.modelpar.append(f"disk_radius_{i}")
+
+                        self.bounds[f"disk_teff_{i}"] = bounds["disk_teff"][i]
+                        self.bounds[f"disk_radius_{i}"] = bounds["disk_radius"][i]
+
+                    del bounds["disk_teff"]
+                    del bounds["disk_radius"]
+
+                else:
+                    # Fitting a single disk component
+                    self.n_disk = 1
+
+                    self.modelpar.append("disk_teff")
+                    self.modelpar.append("disk_radius")
+
+            else:
+                self.n_disk = 0
 
             if self.binary:
                 # Update list of model parameters
@@ -847,7 +884,7 @@ def __init__(
         self.diskphot = []
         self.diskspec = []
 
-        if "disk_teff" in self.bounds and "disk_radius" in self.bounds:
+        if self.n_disk > 0:
             for item in inc_phot:
                 print(f"Interpolating {item}...", end="", flush=True)
                 readmodel = ReadModel("blackbody", filter_name=item)
@@ -1182,7 +1219,6 @@ def _lnlike_func(
         corr_len = {}
         corr_amp = {}
         dust_param = {}
-        disk_param = {}
         veil_param = {}
         param_dict = {}
         rad_vel = {}
@@ -1223,12 +1259,6 @@ def _lnlike_func(
             elif self.ext_filter is not None and item == f"phot_ext_{self.ext_filter}":
                 dust_param[item] = params[self.cube_index[item]]
 
-            elif item == "disk_teff":
-                disk_param["teff"] = params[self.cube_index[item]]
-
-            elif item == "disk_radius":
-                disk_param["radius"] = params[self.cube_index[item]]
-
             elif item == "veil_a":
                 veil_param["veil_a"] = params[self.cube_index[item]]
 
@@ -1244,6 +1274,41 @@ def _lnlike_func(
             else:
                 param_dict[item] = params[self.cube_index[item]]
 
+        # Disk parameters
+
+        disk_param = {}
+
+        if self.n_disk == 1:
+            if "disk_teff" in self.fix_param:
+                disk_param["teff"] = self.fix_param["disk_teff"]
+            else:
+                disk_param["teff"] = params[self.cube_index["disk_teff"]]
+
+            if "disk_radius" in self.fix_param:
+                disk_param["radius"] = self.fix_param["disk_radius"]
+            else:
+                disk_param["radius"] = params[self.cube_index["disk_radius"]]
+
+        elif self.n_disk > 1:
+            for disk_idx in range(self.n_disk):
+                if f"disk_teff_{disk_idx}" in self.fix_param:
+                    disk_param[f"teff_{disk_idx}"] = self.fix_param[
+                        f"disk_teff_{disk_idx}"
+                    ]
+                else:
+                    disk_param[f"teff_{disk_idx}"] = params[
+                        self.cube_index[f"disk_teff_{disk_idx}"]
+                    ]
+
+                if f"disk_radius_{disk_idx}" in self.fix_param:
+                    disk_param[f"radius_{disk_idx}"] = self.fix_param[
+                        f"disk_radius_{disk_idx}"
+                    ]
+                else:
+                    disk_param[f"radius_{disk_idx}"] = params[
+                        self.cube_index[f"disk_radius_{disk_idx}"]
+                    ]
+
         # Add the parallax manually because it should
         # not be provided in the bounds dictionary
 
@@ -1285,18 +1350,14 @@ def _lnlike_func(
             elif item[:9] == "phot_ext_":
                 dust_param[item] = self.fix_param[item]
 
-            elif item == "disk_teff":
-                disk_param["teff"] = self.fix_param[item]
-
-            elif item == "disk_radius":
-                disk_param["radius"] = self.fix_param[item]
-
             elif item == "spec_weight":
                 pass
 
             else:
                 param_dict[item] = self.fix_param[item]
 
+        # Check if the blackbody temperatures/radii are decreasing/increasing
+
         if self.model == "planck" and self.n_planck > 1:
             for i in range(self.n_planck - 1):
                 if param_dict[f"teff_{i+1}"] > param_dict[f"teff_{i}"]:
@@ -1305,13 +1366,35 @@ def _lnlike_func(
                 if param_dict[f"radius_{i}"] > param_dict[f"radius_{i+1}"]:
                     return -np.inf
 
-        if disk_param:
+        if self.n_disk == 1:
             if disk_param["teff"] > param_dict["teff"]:
                 return -np.inf
 
             if disk_param["radius"] < param_dict["radius"]:
                 return -np.inf
 
+        elif self.n_disk > 1:
+            for disk_idx in range(self.n_disk):
+                if disk_idx == 0:
+                    if disk_param["teff_0"] > param_dict["teff"]:
+                        return -np.inf
+
+                    if disk_param["radius_0"] < param_dict["radius"]:
+                        return -np.inf
+
+                else:
+                    if (
+                        disk_param[f"teff_{disk_idx}"]
+                        > disk_param[f"teff_{disk_idx-1}"]
+                    ):
+                        return -np.inf
+
+                    if (
+                        disk_param[f"radius_{disk_idx}"]
+                        < disk_param[f"radius_{disk_idx-1}"]
+                    ):
+                        return -np.inf
+
         if self.model != "powerlaw":
             if "radius_0" in param_dict and "radius_1" in param_dict:
                 flux_scaling_0 = (param_dict["radius_0"] * constants.R_JUP) ** 2 / (
@@ -1503,7 +1586,9 @@ def _lnlike_func(
                     phot_flux *= flux_scaling
                     phot_flux += flux_offset
 
-            if disk_param:
+            # Add blackbody flux from disk components
+
+            if self.n_disk == 1:
                 phot_tmp = self.diskphot[i].spectrum_interp([disk_param["teff"]])[0][0]
 
                 phot_flux += (
@@ -1512,6 +1597,20 @@ def _lnlike_func(
                     / (1e3 * constants.PARSEC / parallax) ** 2
                 )
 
+            elif self.n_disk > 1:
+                for disk_idx in range(self.n_disk):
+                    phot_tmp = self.diskphot[i].spectrum_interp(
+                        [disk_param[f"teff_{disk_idx}"]]
+                    )[0][0]
+
+                    phot_flux += (
+                        phot_tmp
+                        * (disk_param[f"radius_{disk_idx}"] * constants.R_JUP) ** 2
+                        / (1e3 * constants.PARSEC / parallax) ** 2
+                    )
+
+            # Apply extinction
+
             if "lognorm_ext" in dust_param:
                 cross_tmp = self.cross_sections[phot_filter](
                     (10.0 ** dust_param["lognorm_radius"], dust_param["lognorm_sigma"])
@@ -1766,7 +1865,9 @@ def _lnlike_func(
                         self.spectrum[item][1] * sigma_i * sigma_j
                     )
 
-            if disk_param:
+            # Add blackbody flux from disk components
+
+            if self.n_disk == 1:
                 model_tmp = self.diskspec[i].spectrum_interp([disk_param["teff"]])[0, :]
 
                 model_tmp *= (disk_param["radius"] * constants.R_JUP) ** 2 / (
@@ -1775,6 +1876,20 @@ def _lnlike_func(
 
                 model_flux += model_tmp
 
+            elif self.n_disk > 1:
+                for disk_idx in range(self.n_disk):
+                    model_tmp = self.diskspec[i].spectrum_interp(
+                        [disk_param[f"teff_{disk_idx}"]]
+                    )[0, :]
+
+                    model_tmp *= (
+                        disk_param[f"radius_{disk_idx}"] * constants.R_JUP
+                    ) ** 2 / (1e3 * constants.PARSEC / parallax) ** 2
+
+                    model_flux += model_tmp
+
+            # Apply extinction
+
             if "lognorm_ext" in dust_param:
                 cross_tmp = self.cross_sections["spectrum"](
                     (
@@ -1821,6 +1936,8 @@ def _lnlike_func(
 
                 model_flux *= 10.0 ** (-0.4 * ext_spec)
 
+            # Calculate the likelihood
+
             if self.spectrum[item][2] is not None:
                 # Use the inverted covariance matrix
                 ln_like += -0.5 * np.dot(