Removed interp_method from configuration so linear interpolation is u…

…sed, added rad_vel and vsini as optional retrieval parameters in AtmosphericRetrieval
tomasstolker · Feb 16, 2024 · 3f07949 · 3f07949
1 parent 57d7eca
commit 3f07949
Show file tree

Hide file tree

Showing 10 changed files with 432 additions and 196 deletions.
diff --git a/docs/tutorials/atmospheric_retrieval.ipynb b/docs/tutorials/atmospheric_retrieval.ipynb
diff --git a/species/core/species_init.py b/species/core/species_init.py
@@ -71,11 +71,11 @@ def __init__(self) -> None:
                 file_obj.write("; Folder where data will be downloaded\n")
                 file_obj.write("data_folder = ./data/\n\n")
 
-                file_obj.write("; Method for the grid interpolation\n")
-                file_obj.write(
-                    "; Options: linear, nearest, slinear, " "cubic, quintic, pchip\n"
-                )
-                file_obj.write("interp_method = linear\n")
+                # file_obj.write("; Method for the grid interpolation\n")
+                # file_obj.write(
+                #     "; Options: linear, nearest, slinear, " "cubic, quintic, pchip\n"
+                # )
+                # file_obj.write("interp_method = linear\n")
 
             print(" [DONE]")
 
@@ -102,18 +102,18 @@ def __init__(self) -> None:
                 file_obj.write("\n; Folder where data will be downloaded\n")
                 file_obj.write("data_folder = ./data/\n")
 
-        if "interp_method" in config["species"]:
-            interp_method = config["species"]["interp_method"]
-
-        else:
-            interp_method = "linear"
-
-            with open(config_file, "a", encoding="utf-8") as file_obj:
-                file_obj.write("\n; Method for the grid interpolation\n")
-                file_obj.write(
-                    "; Options: linear, nearest, slinear, " "cubic, quintic, pchip\n"
-                )
-                file_obj.write("interp_method = linear\n")
+        # if "interp_method" in config["species"]:
+        #     interp_method = config["species"]["interp_method"]
+        #
+        # else:
+        #     interp_method = "linear"
+        #
+        #     with open(config_file, "a", encoding="utf-8") as file_obj:
+        #         file_obj.write("\n; Method for the grid interpolation\n")
+        #         file_obj.write(
+        #             "; Options: linear, nearest, slinear, " "cubic, quintic, pchip\n"
+        #         )
+        #         file_obj.write("interp_method = linear\n")
 
         if "vega_mag" in config["species"]:
             vega_mag = config["species"]["vega_mag"]
@@ -128,7 +128,7 @@ def __init__(self) -> None:
         print("\nConfiguration settings:")
         print(f"   - Database: {database_file}")
         print(f"   - Data folder: {data_folder}")
-        print(f"   - Interpolation method: {interp_method}")
+        # print(f"   - Interpolation method: {interp_method}")
         print(f"   - Magnitude of Vega: {vega_mag}")
 
         if not os.path.isfile(database_file):

diff --git a/species/data/database.py b/species/data/database.py
@@ -1504,7 +1504,7 @@ def add_object(
                     print(f"      - Data shape: {read_cov[spec_item].shape}")
 
             if verbose:
-                print("   - Spectral resolution:")
+                print("   - Resolution:")
 
             for spec_item, spec_value in spectrum.items():
                 hdf5_file.create_dataset(
@@ -2211,7 +2211,7 @@ def get_mcmc_spectra(
         print(f"Database tag: {tag}")
         print(f"Number of samples: {random}")
         print(f"Wavelength range (um): {wavel_range}")
-        print(f"Spectral resolution: {spec_res}")
+        print(f"Resolution: {spec_res}")
 
         if burnin is None:
             burnin = 0
@@ -2603,6 +2603,7 @@ def get_object(
         object_name: str,
         inc_phot: Union[bool, List[str]] = True,
         inc_spec: Union[bool, List[str]] = True,
+        verbose: bool = True,
     ) -> ObjectBox:
         """
         Function for extracting the photometric and/or spectroscopic
@@ -2625,18 +2626,21 @@ def get_object(
             database with
             :func:`~species.data.database.Database.add_object`) can
             be provided.
+        verbose : bool
+            Print output.
 
         Returns
         -------
         species.core.box.ObjectBox
             Box with the object's data.
         """
 
-        print_section(f"Get object")
+        if verbose:
+            print_section(f"Get object")
 
-        print(f"Object name: {object_name}")
-        print(f"Include photometry: {inc_phot}")
-        print(f"Include spectra: {inc_spec}")
+            print(f"Object name: {object_name}")
+            print(f"Include photometry: {inc_phot}")
+            print(f"Include spectra: {inc_spec}")
 
         with h5py.File(self.database, "r") as hdf5_file:
             if f"objects/{object_name}" not in hdf5_file:

diff --git a/species/fit/fit_evolution.py b/species/fit/fit_evolution.py
@@ -135,7 +135,6 @@ def __init__(
 
         self.database = config["species"]["database"]
         self.database_path = config["species"]["database"]
-        self.interp_method = config["species"]["interp_method"]
 
         # Add grid with evolution data
 

diff --git a/species/fit/fit_model.py b/species/fit/fit_model.py
@@ -134,6 +134,11 @@ def __init__(
                  automatic priors will be set for all mandatory
                  parameters.
 
+               - Radial velocity can be included with the ``rad_vel``
+                 parameter (km/s). This parameter will only be relevant
+                 if the radial velocity shift can be spectrally
+                 resolved given the instrument resolution. 
+
                - Rotational broadening can be fitted by including the
                  ``vsini`` parameter (km/s). This parameter will only
                  be relevant if the rotational broadening is stronger
@@ -1700,11 +1705,13 @@ def _lnlike_func(
                 wavel_shift = (
                     rad_vel[item] * 1e3 * self.spectrum[item][0][:, 0] / constants.LIGHT
                 )
+
                 spec_interp = interpolate.interp1d(
                     self.spectrum[item][0][:, 0] + wavel_shift,
                     model_flux,
                     fill_value="extrapolate",
                 )
+
                 model_flux = spec_interp(self.spectrum[item][0][:, 0])
 
             # Apply rotational broadening
@@ -2313,13 +2320,13 @@ def _lnlike_ultranest(params: np.ndarray) -> Union[float, np.float64]:
 
         # Best-fit parameters
 
-        print("\nBest-fit parameters (mean +/- std):")
+        print("\nBest-fit parameters (mean +/- sigma):")
 
         for i, item in enumerate(self.modelpar):
             mean = np.mean(result["samples"][:, i])
-            std = np.std(result["samples"][:, i])
+            sigma = np.std(result["samples"][:, i])
 
-            print(f"   - {item} = {mean:.2e} +/- {std:.2e}")
+            print(f"   - {item} = {mean:.2e} +/- {sigma:.2e}")
 
         # Maximum likelihood sample