Support for only including photometric fluxes when comparing a grid o…

…f model spectra with CompareSpectra

species/analysis/compare_spectra.py

@@ -444,7 +444,11 @@ def compare_model(
             else:
                 w_i[spec_item] = np.ones(obj_wavel.shape[0])
 
+        # Open the HDF5 databse
+        species_db = database.Database()
+
         if isinstance(inc_phot, bool):
+            # The argument of inc_phot is a boolean
             if inc_phot:
                 # Select all filters if inc_phot=True
                 species_db = database.Database()
@@ -454,6 +458,10 @@ def compare_model(
             else:
                 inc_phot = []
 
+        else:
+            # The argument of inc_phot is a list with filters
+            object_box = species_db.get_object(self.object_name, inc_phot=inc_phot)
+
         if scale_spec is None:
             scale_spec = []
 
@@ -592,21 +600,26 @@ def compare_model(
 
                                     model_spec[spec_item] = model_flux
 
+                                g_fit = 0.0
+
                                 model_list = []
                                 data_list = []
                                 weights_list = []
+
                                 for spec_item in self.spec_name:
                                     if spec_item not in scale_spec:
                                         model_list.append(model_spec[spec_item])
                                         data_list.append(spec_data[spec_item][0])
                                         weights_list.append(w_i[spec_item])
 
+                                model_phot = {}
+
                                 for phot_item in inc_phot:
                                     syn_phot = photometry.SyntheticPhotometry(phot_item)
-                                    model_phot = syn_phot.spectrum_to_flux(
+                                    model_phot[phot_item] = syn_phot.spectrum_to_flux(
                                         model_box_full.wavelength, model_box_full.flux
-                                    )
-                                    model_list.append(np.array([model_phot[0]]))
+                                    )[0]
+                                    model_list.append(np.array([model_phot[phot_item]]))
 
                                     phot_flux = object_box.flux[phot_item]
                                     phot_data = np.array(
@@ -632,6 +645,7 @@ def compare_model(
                                     * model_list
                                     / data_list[:, 2] ** 2
                                 )
+
                                 c_denom = (
                                     weights_list
                                     * model_list**2
@@ -689,7 +703,18 @@ def compare_model(
                                         spec_idx,
                                     ] = spec_scaling
 
-                                g_fit = 0.0
+                                for phot_item in inc_phot:
+                                    phot_flux = object_box.flux[phot_item]
+
+                                    g_fit += np.nansum(
+                                        w_i[phot_item]
+                                        * (
+                                            phot_flux[0]
+                                            - scaling * model_phot[phot_item]
+                                        )
+                                        ** 2
+                                        / phot_flux[1] ** 2
+                                    )
 
                                 for spec_item in self.spec_name:
                                     if spec_item in scale_spec:
@@ -766,4 +791,5 @@ def compare_model(
             model=model,
             scale_spec=scale_spec,
             extra_scaling=extra_scaling,
+            inc_phot=inc_phot,
         )

species/data/database.py

@@ -2567,21 +2567,26 @@ def get_object(
         inc_spec: Union[bool, List[str]] = True,
     ) -> box.ObjectBox:
         """
-        Function for extracting the photometric and/or spectroscopic data of an object from the
-        database. The spectroscopic data contains optionally the covariance matrix and its inverse.
+        Function for extracting the photometric and/or spectroscopic
+        data of an object from the database. The spectroscopic data
+        contains optionally the covariance matrix and its inverse.
 
         Parameters
         ----------
         object_name : str
             Object name in the database.
         inc_phot : bool, list(str)
-            Include photometric data. If a boolean, either all (``True``) or none (``False``) of
-            the data are selected. If a list, a subset of filter names (as stored in the database)
-            can be provided.
+            Include photometric data. If a boolean, either all
+            (``True``) or none (``False``) of the data are selected.
+            If a list, a subset of filter names (as stored in the
+            database) can be provided.
         inc_spec : bool, list(str)
-            Include spectroscopic data. If a boolean, either all (``True``) or none (``False``) of
-            the data are selected. If a list, a subset of spectrum names (as stored in the database
-            with :func:`~species.data.database.Database.add_object`) can be provided.
+            Include spectroscopic data. If a boolean, either all
+            (``True``) or none (``False``) of the data are selected.
+            If a list, a subset of spectrum names (as stored in the
+            database with
+            :func:`~species.data.database.Database.add_object`) can
+            be provided.
 
         Returns
         -------
@@ -3064,6 +3069,7 @@ def add_comparison(
         model: str,
         scale_spec: List[str],
         extra_scaling: Optional[np.ndarray],
+        inc_phot: List[str],
     ) -> None:
         """
         Function for adding results obtained with
@@ -3100,6 +3106,9 @@ def add_comparison(
             Array with extra scalings that have been applied to the
             spectra of ``scale_spec``. The argument can be set to
             ``None`` if no extra scalings have been applied.
+        inc_phot : list(str)
+            List with filter names of which photometric data
+            was included with the comparison.
 
         Returns
         -------
@@ -3130,6 +3139,7 @@ def add_comparison(
             dset.attrs["n_spec_name"] = len(spec_name)
             dset.attrs["n_scale_spec"] = len(scale_spec)
             dset.attrs["parallax"] = parallax
+            dset.attrs["n_inc_phot"] = len(inc_phot)
 
             for i, item in enumerate(model_param):
                 dset.attrs[f"parameter{i}"] = item
@@ -3140,6 +3150,9 @@ def add_comparison(
             for i, item in enumerate(scale_spec):
                 dset.attrs[f"scale_spec{i}"] = item
 
+            for i, item in enumerate(inc_phot):
+                dset.attrs[f"inc_phot{i}"] = item
+
             h5_file.create_dataset(
                 f"results/comparison/{tag}/flux_scaling", data=flux_scaling
             )

species/data/model_data.json

@@ -138,7 +138,7 @@
         "wavelength range": [0.35, 250.0],
         "resolution": 500,
         "teff range": [400, 2000],
-        "reference": "Charney et al. (2018)",
+        "reference": "Charnay et al. (2018)",
         "url": "https://ui.adsabs.harvard.edu/abs/2018ApJ...854..172C/abstract"
     },
     "exo-rem-highres": {
@@ -148,7 +148,7 @@
         "wavelength range": [0.67, 250.0],
         "resolution": 20000,
         "teff range": [400, 2000],
-        "reference": "Charney et al. (2018)",
+        "reference": "Charnay et al. (2018)",
         "url": "https://ui.adsabs.harvard.edu/abs/2018ApJ...854..172C/abstract"
     },
     "koester-wd": {

species/plot/plot_comparison.py

@@ -963,10 +963,10 @@ def plot_model_spectra(
     flux_offset : float, None
         Offset to be applied such that the spectra do not overlap. No
         offset is applied if the argument is set to ``None``.
-    xlim : tuple(float, float)
-        Limits of the spectral type axis.
-    ylim : tuple(float, float)
-        Limits of the goodness-of-fit axis.
+    xlim : tuple(float, float), None
+        Limits of the wavelength axis.
+    ylim : tuple(float, float), None
+        Limits of the flux axis.
     title : str
         Plot title.
     offset : tuple(float, float)
@@ -1016,11 +1016,16 @@ def plot_model_spectra(
     n_param = dset.attrs["n_param"]
     n_scale_spec = dset.attrs["n_scale_spec"]
     parallax = dset.attrs["parallax"]
+    n_inc_phot = dset.attrs["n_inc_phot"]
 
     spec_name = []
     for i in range(n_spec_name):
         spec_name.append(dset.attrs[f"spec_name{i}"])
 
+    inc_phot = []
+    for i in range(n_inc_phot):
+        inc_phot.append(dset.attrs[f"inc_phot{i}"])
+
     model_param = []
     coord_points = []
     for i in range(n_param):
@@ -1174,6 +1179,19 @@ def plot_model_spectra(
                     color="black",
                 )
 
+        if n_inc_phot > 0 and n_spec_name == 0:
+            model_box = model_reader.get_data(
+                model_param=param_select,
+                spec_res=100.,
+            )
+
+            ax.plot(
+                model_box.wavelength,
+                (n_spectra - i - 1) * flux_offset + model_box.flux,
+                color="tomato",
+                lw=0.5,
+            )
+
         for spec_key, spec_item in obj_spec.items():
             if spec_key not in spec_name:
                 continue

species/plot/plot_spectrum.py

@@ -1057,7 +1057,11 @@ def plot_spectrum(
                                 )
 
                 finite = np.isfinite(residuals.photometry[item][1])
-                res_max = np.max(np.abs(residuals.photometry[item][1][finite]))
+
+                max_tmp = np.max(np.abs(residuals.photometry[item][1][finite]))
+
+                if max_tmp > res_max:
+                    res_max = max_tmp
 
         if residuals.spectrum is not None:
             for key, value in residuals.spectrum.items():

species/util/plot_util.py

@@ -1392,6 +1392,7 @@ def create_model_label(
     if model_name in model_list:
         read_mod = read_model.ReadModel(model_name)
         check_param = read_mod.get_parameters()
+        check_param += ['radius', 'mass', 'luminosity']
 
     else:
         check_param = None