Fixed issue with flux scaling when using flux ratio prior, including …

…the luminosity of both stars in plot_posterior when fitting the combined SED
tomasstolker · May 3, 2024 · a1cd551 · a1cd551
1 parent 01f917c
commit a1cd551
Show file tree

Hide file tree

Showing 2 changed files with 52 additions and 10 deletions.
diff --git a/species/fit/fit_model.py b/species/fit/fit_model.py
@@ -1143,7 +1143,8 @@ def __init__(
 
                     else:
                         print(
-                            f"   - {spec_item} = {np.amin(self.weights[spec_item]):.2e} - {np.amax(self.weights[spec_item]):.2e}"
+                            f"   - {spec_item} = {np.amin(self.weights[spec_item]):.2e} "
+                            f"- {np.amax(self.weights[spec_item]):.2e}"
                         )
 
                 for phot_item in inc_phot:
@@ -1186,7 +1187,8 @@ def __init__(
 
                 else:
                     print(
-                        f"   - {spec_item} = {np.amin(self.weights[spec_item]):.2e} - {np.amax(self.weights[spec_item]):.2e}"
+                        f"   - {spec_item} = {np.amin(self.weights[spec_item]):.2e} "
+                        f"- {np.amax(self.weights[spec_item]):.2e}"
                     )
 
             for phot_item in inc_phot:
@@ -1563,22 +1565,28 @@ def _lnlike_func(
                         )
 
             elif key[:6] == "ratio_":
-                filt_name = key[6:]
+                filter_name = key[6:]
 
                 param_0 = binary_to_single(param_dict, 0)
-                phot_flux_0 = self.flux_ratio[filt_name].spectrum_interp(
+
+                phot_flux_0 = self.flux_ratio[filter_name].spectrum_interp(
                     list(param_0.values())
                 )[0][0]
 
+                phot_flux_0 *= flux_scaling_0
+
                 param_1 = binary_to_single(param_dict, 1)
-                phot_flux_1 = self.flux_ratio[filt_name].spectrum_interp(
+
+                phot_flux_1 = self.flux_ratio[filter_name].spectrum_interp(
                     list(param_1.values())
                 )[0][0]
 
+                phot_flux_1 *= flux_scaling_1
+
                 # Uniform prior for the flux ratio
 
-                if f"ratio_{filt_name}" in self.bounds:
-                    ratio_prior = self.bounds[f"ratio_{filt_name}"]
+                if f"ratio_{filter_name}" in self.bounds:
+                    ratio_prior = self.bounds[f"ratio_{filter_name}"]
 
                     if ratio_prior[0] > phot_flux_1 / phot_flux_0:
                         return -np.inf
@@ -1587,8 +1595,8 @@ def _lnlike_func(
 
                 # Normal prior for the flux ratio
 
-                if f"ratio_{filt_name}" in self.normal_prior:
-                    ratio_prior = self.normal_prior[f"ratio_{filt_name}"]
+                if f"ratio_{filter_name}" in self.normal_prior:
+                    ratio_prior = self.normal_prior[f"ratio_{filter_name}"]
 
                     ln_like += (
                         -0.5

diff --git a/species/plot/plot_mcmc.py b/species/plot/plot_mcmc.py
@@ -608,7 +608,41 @@ def plot_posterior(
                 box.parameters.append("luminosity")
                 ndim += 1
 
-        elif "teff_0" in box.parameters and "radius_0" in box.parameters:
+        if "teff_0" in box.parameters and "radius_0" in box.parameters:
+            teff_index = np.argwhere(np.array(box.parameters) == "teff_0")
+            radius_index = np.argwhere(np.array(box.parameters) == "radius_0")
+
+            luminosity = (
+                4.0
+                * np.pi
+                * (samples[..., radius_index[0]] * constants.R_JUP) ** 2
+                * constants.SIGMA_SB
+                * samples[..., teff_index[0]] ** 4.0
+                / constants.L_SUN
+            )
+
+            samples = np.append(samples, np.log10(luminosity), axis=-1)
+            box.parameters.append("luminosity_0")
+            ndim += 1
+
+        if "teff_1" in box.parameters and "radius_1" in box.parameters:
+            teff_index = np.argwhere(np.array(box.parameters) == "teff_1")
+            radius_index = np.argwhere(np.array(box.parameters) == "radius_1")
+
+            luminosity = (
+                4.0
+                * np.pi
+                * (samples[..., radius_index[0]] * constants.R_JUP) ** 2
+                * constants.SIGMA_SB
+                * samples[..., teff_index[0]] ** 4.0
+                / constants.L_SUN
+            )
+
+            samples = np.append(samples, np.log10(luminosity), axis=-1)
+            box.parameters.append("luminosity_1")
+            ndim += 1
+
+        if "teff_0" in box.parameters and "radius_0" in box.parameters:
             luminosity = 0.0
 
             for i in range(100):