Support for photometric upper limits, fix in adding the luminosity sa…

…mples of multiple Planck functions
tomasstolker · Jul 2, 2020 · ac445dd · ac445dd
1 parent 617768a
commit ac445dd
Show file tree

Hide file tree

Showing 3 changed files with 16 additions and 5 deletions.
diff --git a/species/data/database.py b/species/data/database.py
@@ -1157,6 +1157,8 @@ def get_median_sample(self,
 
             samples = np.asarray(dset)
 
+            # samples = samples[samples[:, 2] > 100., ]
+
             if samples.ndim == 3:
                 if burnin > samples.shape[1]:
                     raise ValueError(f'The \'burnin\' value is larger than the number of steps '
@@ -1251,6 +1253,8 @@ def get_mcmc_spectra(self,
 
         samples = np.asarray(dset)
 
+        # samples = samples[samples[:, 2] > 100., ]
+
         if samples.ndim == 2:
             ran_index = np.random.randint(samples.shape[0], size=random)
             samples = samples[ran_index, ]

diff --git a/species/plot/plot_mcmc.py b/species/plot/plot_mcmc.py
@@ -206,9 +206,9 @@ def plot_posterior(tag: str,
                 else:
                     break
 
-            # samples = np.append(samples, np.log10(luminosity), axis=-1)
-            # box.parameters.append('luminosity')
-            # ndim += 1
+            samples = np.append(samples, np.log10(luminosity), axis=-1)
+            box.parameters.append('luminosity')
+            ndim += 1
 
             # teff_index = np.argwhere(np.array(box.parameters) == 'teff_0')
             # radius_index = np.argwhere(np.array(box.parameters) == 'radius_0')

diff --git a/species/plot/plot_spectrum.py b/species/plot/plot_spectrum.py
@@ -545,8 +545,15 @@ def plot_spectrum(boxes: list,
                                              yerr=flux_scaling*boxitem.flux[item][1, i]/scaling, zorder=3, **plot_kwargs[j][item][i])
 
                         else:
-                            ax1.errorbar(wavelength, flux_scaling*boxitem.flux[item][0]/scaling, xerr=fwhm/2.,
-                                         yerr=flux_scaling*boxitem.flux[item][1]/scaling, zorder=3, **plot_kwargs[j][item])
+                            if boxitem.flux[item][1] == 0.:
+                                ax1.errorbar(wavelength, flux_scaling*boxitem.flux[item][0]/scaling,
+                                             xerr=fwhm/2., yerr=0.5*flux_scaling*boxitem.flux[item][0]/scaling,
+                                             uplims=True, capsize=2., capthick=0., zorder=3, **plot_kwargs[j][item])
+
+                            else:
+                                ax1.errorbar(wavelength, flux_scaling*boxitem.flux[item][0]/scaling,
+                                             xerr=fwhm/2., yerr=flux_scaling*boxitem.flux[item][1]/scaling,
+                                             zorder=3, **plot_kwargs[j][item])
 
         elif isinstance(boxitem, box.SynphotBox):
             for i, find_item in enumerate(boxes):