Fitting the SED of an unresolved binary star with different parallaxe…

…s, several additional improvements for fitting a combination of two spectra

species/data/database.py

@@ -2339,7 +2339,7 @@ def get_probable_sample(
 
                 prob_sample[par_key] = par_value
 
-            if "parallax" not in prob_sample and "parallax" in dset.attrs:
+            if "parallax" not in prob_sample and "parallax_0" not in prob_sample and "parallax" in dset.attrs:
                 prob_sample["parallax"] = dset.attrs["parallax"]
 
             elif "distance" not in prob_sample and "distance" in dset.attrs:
@@ -2426,7 +2426,7 @@ def get_median_sample(
                 par_value = np.median(samples[:, i])
                 median_sample[par_key] = par_value
 
-            if "parallax" not in median_sample and "parallax" in dset.attrs:
+            if "parallax" not in median_sample and "parallax_0" not in median_sample and "parallax" in dset.attrs:
                 median_sample["parallax"] = dset.attrs["parallax"]
 
             elif "distance" not in median_sample and "distance" in dset.attrs:
@@ -2686,7 +2686,7 @@ def get_mcmc_spectra(
                 if param[j] not in ignore_param:
                     model_param[param[j]] = samples[i, j]
 
-            if "parallax" not in model_param and parallax is not None:
+            if "parallax" not in model_param and "parallax_0" not in model_param and parallax is not None:
                 model_param["parallax"] = parallax
 
             elif "distance" not in model_param and distance is not None:
@@ -2733,10 +2733,17 @@ def get_mcmc_spectra(
                             ext_filter=ext_filter,
                         )
 
-                        flux_comb = (
-                            model_param["spec_weight"] * specbox_0.flux
-                            + (1.0 - model_param["spec_weight"]) * specbox_1.flux
-                        )
+                        # Weighted flux of two spectra for atmospheric asymmetries
+                        # Or simply the same in case of an actual binary system
+
+                        if "spec_weight" in model_param:
+                            flux_comb = (
+                                model_param["spec_weight"] * specbox_0.flux
+                                + (1.0 - model_param["spec_weight"]) * specbox_1.flux
+                            )
+
+                        else:
+                            flux_comb = specbox_0.flux + specbox_1.flux
 
                         specbox = create_box(
                             boxtype="model",
@@ -2880,7 +2887,7 @@ def get_mcmc_photometry(
             for j in range(n_param):
                 model_param[param[j]] = samples[i, j]
 
-            if "parallax" not in model_param and parallax is not None:
+            if "parallax" not in model_param and "parallax_0" not in model_param and parallax is not None:
                 model_param["parallax"] = parallax
 
             elif "distance" not in model_param and distance is not None:
@@ -2914,10 +2921,17 @@ def get_mcmc_photometry(
                             param_1 = binary_to_single(model_param, 1)
                             mcmc_phot_1, _ = readmodel.get_magnitude(param_1)
 
-                            mcmc_phot[i] = (
-                                model_param["spec_weight"] * mcmc_phot_0
-                                + (1.0 - model_param["spec_weight"]) * mcmc_phot_1
-                            )
+                            # Weighted flux of two spectra for atmospheric asymmetries
+                            # Or simply the same in case of an actual binary system
+
+                            if "spec_weight" in model_param:
+                                mcmc_phot[i] = (
+                                    model_param["spec_weight"] * mcmc_phot_0
+                                    + (1.0 - model_param["spec_weight"]) * mcmc_phot_1
+                                )
+
+                            else:
+                                mcmc_phot[i] = mcmc_phot_0 + mcmc_phot_1
 
                         else:
                             mcmc_phot[i], _ = readmodel.get_magnitude(model_param)
@@ -2932,10 +2946,17 @@ def get_mcmc_photometry(
                             param_1 = binary_to_single(model_param, 1)
                             mcmc_phot_1, _ = readmodel.get_flux(param_1)
 
-                            mcmc_phot[i] = (
-                                model_param["spec_weight"] * mcmc_phot_0
-                                + (1.0 - model_param["spec_weight"]) * mcmc_phot_1
-                            )
+                            # Weighted flux of two spectra for atmospheric asymmetries
+                            # Or simply the same in case of an actual binary system
+
+                            if "spec_weight" in model_param:
+                                mcmc_phot[i] = (
+                                    model_param["spec_weight"] * mcmc_phot_0
+                                    + (1.0 - model_param["spec_weight"]) * mcmc_phot_1
+                                )
+
+                            else:
+                                mcmc_phot[i] = mcmc_phot_0 + mcmc_phot_1
 
                         else:
                             mcmc_phot[i], _ = readmodel.get_flux(model_param)

species/fit/fit_model.py

@@ -472,7 +472,7 @@ def __init__(
         # Set attributes
 
         self.object = ReadObject(object_name)
-        self.parallax = self.object.get_parallax()
+        self.obj_parallax = self.object.get_parallax()
         self.binary = False
         self.ext_filter = ext_filter
 
@@ -545,12 +545,6 @@ def __init__(
                         self.bounds[f"{key}_1"] = bounds_grid[key]
                         del self.bounds[key]
 
-                    elif isinstance(self.bounds[key][0], tuple):
-                        self.binary = True
-                        self.bounds[f"{key}_0"] = self.bounds[key][0]
-                        self.bounds[f"{key}_1"] = self.bounds[key][1]
-                        del self.bounds[key]
-
                     else:
                         if self.bounds[key][0] < bounds_grid[key][0]:
                             warnings.warn(
@@ -620,6 +614,10 @@ def __init__(
                 readmodel = ReadModel(self.model, None, None)
                 self.bounds = readmodel.get_bounds()
 
+            print(f"Object name: {object_name}")
+            print(f"Model tag: {model}")
+            print(f"Binary star: {self.binary}")
+
             self.modelpar = readmodel.get_parameters()
 
             if "flux_scaling" in self.bounds:
@@ -634,7 +632,22 @@ def __init__(
 
             else:
                 self.modelpar.append("radius")
-                self.modelpar.append("parallax")
+
+                if self.binary:
+                    if "parallax" in self.bounds:
+                        if isinstance(self.bounds["parallax"][0], tuple):
+                            self.modelpar.append("parallax_0")
+                            self.modelpar.append("parallax_1")
+                            self.bounds["parallax_0"] = self.bounds["parallax"][0]
+                            self.bounds["parallax_1"] = self.bounds["parallax"][1]
+                            del self.bounds["parallax"]
+
+                    if "parallax_0" in self.normal_prior:
+                            self.modelpar.append("parallax_0")
+                            self.modelpar.append("parallax_1")
+
+                if "parallax_0" not in self.modelpar:
+                    self.modelpar.append("parallax")
 
             if "flux_offset" in self.bounds:
                 self.modelpar.append("flux_offset")
@@ -698,10 +711,16 @@ def __init__(
                         self.modelpar[par_index] = key[:-2] + "_0"
                         self.modelpar.insert(par_index, key[:-2] + "_1")
 
-                self.modelpar.append("spec_weight")
+
+                if "radius" in self.modelpar:
+                    # Fit a weighting for the two spectra in case this
+                    # is a single object, so not an actual binary star.
+                    # In that case the combination of two spectra is
+                    # used to account for atmospheric assymetries
+                    self.modelpar.append("spec_weight")
 
-                if "spec_weight" not in self.bounds:
-                    self.bounds["spec_weight"] = (0.0, 1.0)
+                    if "spec_weight" not in self.bounds:
+                        self.bounds["spec_weight"] = (0.0, 1.0)
 
         # Select filters and spectra
 
@@ -737,6 +756,8 @@ def __init__(
         self.filter_name = []
         self.instr_name = []
 
+        print()
+
         for item in inc_phot:
             if self.model == "planck":
                 # Create SyntheticPhotometry objects when fitting a Planck function
@@ -1035,8 +1056,8 @@ def __init__(
 
         # Add parallax to dictionary with Gaussian priors
 
-        if "parallax" in self.modelpar and "parallax" not in self.fix_param:
-            self.normal_prior["parallax"] = (self.parallax[0], self.parallax[1])
+        if "parallax" in self.modelpar and "parallax" not in self.fix_param and "parallax" not in self.bounds:
+            self.normal_prior["parallax"] = (self.obj_parallax[0], self.obj_parallax[1])
 
         # Printing uniform and normal priors
 
@@ -1313,6 +1334,20 @@ def _lnlike_func(
         # not be provided in the bounds dictionary
 
         if self.model != "powerlaw":
+            if "parallax_0" in self.cube_index:
+                parallax_0 = params[self.cube_index["parallax_0"]]
+            elif "parallax_0" in self.fix_param:
+                parallax_0 = self.fix_param["parallax_0"]
+            else:
+                parallax_0 = params[self.cube_index["parallax"]]
+
+            if "parallax_1" in self.cube_index:
+                parallax_1 = params[self.cube_index["parallax_1"]]
+            elif "parallax_0" in self.fix_param:
+                parallax_1 = self.fix_param["parallax_1"]
+            else:
+                parallax_1 = params[self.cube_index["parallax"]]
+
             if "parallax" in self.cube_index:
                 parallax = params[self.cube_index["parallax"]]
             elif "parallax" in self.fix_param:
@@ -1398,17 +1433,19 @@ def _lnlike_func(
         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 / (
-                    1e3 * constants.PARSEC / parallax
+                    1e3 * constants.PARSEC / parallax_0
                 ) ** 2
 
                 flux_scaling_1 = (param_dict["radius_1"] * constants.R_JUP) ** 2 / (
-                    1e3 * constants.PARSEC / parallax
+                    1e3 * constants.PARSEC / parallax_1
                 ) ** 2
 
                 # The scaling is applied manually because of the interpolation
                 del param_dict["radius_0"]
                 del param_dict["radius_1"]
 
+                flux_offset = 0.0
+
             else:
                 if parallax is None:
                     if "flux_scaling" in self.cube_index:
@@ -1571,12 +1608,17 @@ def _lnlike_func(
                         phot_flux_1 *= flux_scaling_1
                         phot_flux_1 += flux_offset
 
-                    # Weighted flux of two stars
+                    # Weighted flux of two spectra for atmospheric asymmetries
+                    # Or simply the same in case of an actual binary system
 
-                    phot_flux = (
-                        params[self.cube_index["spec_weight"]] * phot_flux_0
-                        + (1.0 - params[self.cube_index["spec_weight"]]) * phot_flux_1
-                    )
+                    if "spec_weight" in self.cube_index:
+                        phot_flux = (
+                            params[self.cube_index["spec_weight"]] * phot_flux_0
+                            + (1.0 - params[self.cube_index["spec_weight"]]) * phot_flux_1
+                        )
+
+                    else:
+                        phot_flux = phot_flux_0 + phot_flux_1
 
                 else:
                     phot_flux = self.modelphot[i].spectrum_interp(
@@ -1763,12 +1805,17 @@ def _lnlike_func(
                         model_flux_1 *= flux_scaling_1
                         model_flux_1 += flux_offset
 
-                    # Weighted flux of two stars
+                    # Weighted flux of two spectra for atmospheric asymmetries
+                    # Or simply the same in case of an actual binary system
+
+                    if "spec_weight" in self.cube_index:
+                        model_flux = (
+                            params[self.cube_index["spec_weight"]] * model_flux_0
+                            + (1.0 - params[self.cube_index["spec_weight"]]) * model_flux_1
+                        )
+                    else:
+                        model_flux = model_flux_0 + model_flux_1
 
-                    model_flux = (
-                        params[self.cube_index["spec_weight"]] * model_flux_0
-                        + (1.0 - params[self.cube_index["spec_weight"]]) * model_flux_1
-                    )
 
                 else:
                     model_flux = self.modelspec[i].spectrum_interp(
@@ -2236,7 +2283,7 @@ def _lnlike_multinest(
             "spec_type": "model",
             "spec_name": self.model,
             "ln_evidence": (ln_z, ln_z_error),
-            "parallax": self.parallax[0],
+            "parallax": self.obj_parallax[0],
         }
 
         if self.ext_filter is not None:
@@ -2502,7 +2549,7 @@ def _lnlike_ultranest(params: np.ndarray) -> Union[float, np.float64]:
             "spec_type": "model",
             "spec_name": self.model,
             "ln_evidence": (ln_z, ln_z_error),
-            "parallax": self.parallax[0],
+            "parallax": self.obj_parallax[0],
         }
 
         if self.ext_filter is not None:
@@ -2865,7 +2912,7 @@ def run_dynesty(
             "spec_type": "model",
             "spec_name": self.model,
             "ln_evidence": (ln_z, ln_z_error),
-            "parallax": self.parallax[0],
+            "parallax": self.obj_parallax[0],
         }
 
         if self.ext_filter is not None:

species/plot/plot_mcmc.py

@@ -725,6 +725,14 @@ def plot_posterior(
         if object_type == "star":
             samples[:, radius_index] *= constants.R_JUP / constants.R_SUN
 
+    for radius_idx in range(100):
+        if f"radius_{radius_idx}" in box.parameters:
+            radius_index = np.argwhere(np.array(box.parameters) == f"radius_{radius_idx}")[0]
+            if object_type == "star":
+                samples[:, radius_index] *= constants.R_JUP / constants.R_SUN
+        else:
+            break
+
     if "mass" in box.parameters:
         mass_index = np.argwhere(np.array(box.parameters) == "mass")[0]
         if object_type == "star":

species/read/read_model.py

@@ -1948,10 +1948,17 @@ def binary_spectrum(
             wavel_resample=wavel_resample,
         )
 
-        flux_comb = (
-            model_param["spec_weight"] * model_box_0.flux
-            + (1.0 - model_param["spec_weight"]) * model_box_1.flux
-        )
+        # Weighted flux of two spectra for atmospheric asymmetries
+        # Or simply the same in case of an actual binary system
+
+        if "spec_weight" in model_param:
+            flux_comb = (
+                model_param["spec_weight"] * model_box_0.flux
+                + (1.0 - model_param["spec_weight"]) * model_box_1.flux
+            )
+
+        else:
+            flux_comb = model_box_0.flux + model_box_1.flux
 
         model_box = create_box(
             boxtype="model",

species/util/fit_util.py

@@ -127,10 +127,17 @@ def multi_photometry(
                         param_1 = binary_to_single(parameters, 1)
                         model_flux_1 = readmodel.get_flux(param_1)[0]
 
-                        flux[item] = (
-                            parameters["spec_weight"] * model_flux_0
-                            + (1.0 - parameters["spec_weight"]) * model_flux_1
-                        )
+                        # Weighted flux of two spectra for atmospheric asymmetries
+                        # Or simply the same in case of an actual binary system
+
+                        if "spec_weight" in parameters:
+                            flux[item] = (
+                                parameters["spec_weight"] * model_flux_0
+                                + (1.0 - parameters["spec_weight"]) * model_flux_1
+                            )
+
+                        else:
+                            flux[item] = model_flux_0 + model_flux_1
 
                     else:
                         # Single object
@@ -368,10 +375,17 @@ def get_residuals(
                             wavel_resample=wl_new,
                         )
 
-                        flux_comb = (
-                            parameters["spec_weight"] * model_spec_0.flux
-                            + (1.0 - parameters["spec_weight"]) * model_spec_1.flux
-                        )
+                        # Weighted flux of two spectra for atmospheric asymmetries
+                        # Or simply the same in case of an actual binary system
+
+                        if "spec_weight" in parameters:
+                            flux_comb = (
+                                parameters["spec_weight"] * model_spec_0.flux
+                                + (1.0 - parameters["spec_weight"]) * model_spec_1.flux
+                            )
+
+                        else:
+                            flux_comb = model_spec_0.flux + model_spec_1.flux
 
                         model_spec = create_box(
                             boxtype="model",