Retrieve the f_sed parameter for each cloud species individually by s…

…etting global_fsed=False

species/fit/retrieval.py

@@ -692,14 +692,19 @@ def _set_parameters(
                 self.parameters.append("albedo")
 
         elif len(self.cloud_species) > 0:
-            self.parameters.append("fsed")
+            if self.global_fsed:
+                self.parameters.append("fsed")
+
             self.parameters.append("log_kzz")
             self.parameters.append("sigma_lnorm")
 
             for item in self.cloud_species:
                 if "log_p_base" in bounds:
                     self.parameters.append(f"log_p_base_{item}")
 
+                if not self.global_fsed:
+                    self.parameters.append(f"fsed_{item}")
+
                 cloud_lower = item[:-3].lower()
 
                 if f"{cloud_lower}_tau" in bounds:
@@ -1475,16 +1480,33 @@ def _prior_transform(
             # mixing velocities of the cloud particles
             # (used in Eq. 3 of Mollière et al. 2020)
 
-            if "fsed" in bounds:
-                fsed = (
-                    bounds["fsed"][0]
-                    + (bounds["fsed"][1] - bounds["fsed"][0]) * cube[cube_index["fsed"]]
-                )
+            if self.global_fsed:
+                if "fsed" in bounds:
+                    fsed = (
+                        bounds["fsed"][0]
+                        + (bounds["fsed"][1] - bounds["fsed"][0])
+                        * cube[cube_index["fsed"]]
+                    )
+                else:
+                    # Default: 0 - 10
+                    fsed = 10.0 * cube[cube_index["fsed"]]
+
+                cube[cube_index["fsed"]] = fsed
+
             else:
-                # Default: 0 - 10
-                fsed = 10.0 * cube[cube_index["fsed"]]
+                for item in self.cloud_species:
+                    if "fsed" in bounds:
+                        fsed = (
+                            bounds["fsed"][0]
+                            + (bounds["fsed"][1] - bounds["fsed"][0])
+                            * cube[cube_index[f"fsed_{item}"]]
+                        )
 
-            cube[cube_index["fsed"]] = fsed
+                    else:
+                        # Default: 0 - 10
+                        fsed = 10.0 * cube[cube_index[f"fsed_{item}"]]
+
+                    cube[cube_index[f"fsed_{item}"]] = fsed
 
             # Log10 of the eddy diffusion coefficient (cm2 s-1)
 
@@ -2183,7 +2205,7 @@ def _lnlike(
 
             # Create dictionary with cloud parameters
 
-            if "fsed" in self.parameters:
+            if "log_kzz" in self.parameters:
                 cloud_param = [
                     "fsed",
                     "log_kzz",
@@ -2204,15 +2226,15 @@ def _lnlike(
                     if item in self.parameters:
                         cloud_dict[item] = cube[cube_index[item]]
 
-                    # elif item in ['log_kzz', 'sigma_lnorm']:
-                    #     cloud_dict[item] = None
-
                 for item in self.cloud_species:
                     if f"log_p_base_{item}" in self.parameters:
                         cloud_dict[f"log_p_base_{item}"] = cube[
                             cube_index[f"log_p_base_{item}"]
                         ]
 
+                    if f"fsed_{item}" in self.parameters:
+                        cloud_dict[f"fsed_{item}"] = cube[cube_index[f"fsed_{item}"]]
+
             elif "fsed_1" in self.parameters and "fsed_2" in self.parameters:
                 cloud_param_1 = [
                     "fsed_1",
@@ -3126,6 +3148,7 @@ def setup_retrieval(
         check_phot_press: Optional[float] = None,
         apply_rad_vel: Optional[List[str]] = None,
         apply_vsini: Optional[List[str]] = None,
+        global_fsed: bool = True,
     ) -> None:
         """
         Function for running the atmospheric retrieval. The parameter
@@ -3376,6 +3399,11 @@ def setup_retrieval(
             ``bounds``. The :math:`v \\sin(i)` is applied to all
             spectra by setting the argument of ``apply_vsini``
             to ``None``.
+        global_fsed : bool
+            Retrieve a global ``fsed`` parameter when set to ``True``
+            or retrieve the ``fsed`` parameter for each cloud species
+            individually in the ``cloud_species`` list when set to
+            ``False``.
 
         Returns
         -------
@@ -3398,6 +3426,7 @@ def setup_retrieval(
         self.cross_corr = cross_corr
         self.apply_rad_vel = apply_rad_vel
         self.apply_vsini = apply_vsini
+        self.global_fsed = global_fsed
 
         print(f"P-T profile: {self.pt_profile}")
         print(f"Chemistry: {self.chemistry}")

species/util/retrieval_util.py

@@ -1210,10 +1210,14 @@ def calc_spectrum_clouds(
 
             abund_in[f"{cloud_item}(c)"] = np.zeros_like(temperature)
 
+            if "fsed" in cloud_dict:
+                f_sed = cloud_dict["fsed"]
+            else:
+                f_sed = cloud_dict[f"fsed_{cloud_item}(c)"]
+
             abund_in[f"{cloud_item}(c)"][pressure < p_base_item] = (
                 10.0 ** log_x_base[cloud_item]
-                * (pressure[pressure <= p_base_item] / p_base_item)
-                ** cloud_dict["fsed"]
+                * (pressure[pressure <= p_base_item] / p_base_item) ** f_sed
             )
 
     # Adaptive pressure refinement around the cloud base
@@ -1237,12 +1241,14 @@ def calc_spectrum_clouds(
     # Use the same value for all cloud species
 
     fseds = {}
-    for item in rt_object.cloud_species:
-        # The item has the form of e.g. MgSiO3(c)
-        # For parametrized cloud opacities,
-        # then number of cloud_species is zero
-        # so the fseds dictionary remains empty
-        fseds[item] = cloud_dict["fsed"]
+    for cloud_item in rt_object.cloud_species:
+        # The cloud_item has the form of e.g. MgSiO3(c). For
+        # parametrized cloud opacities, then the number of
+        # cloud_species is zero so the fseds dictionary is empty
+        if "fsed" in cloud_dict:
+            fseds[cloud_item] = cloud_dict["fsed"]
+        else:
+            fseds[cloud_item] = cloud_dict[f"fsed_{cloud_item}"]
 
     # Create an array with a constant eddy diffusion coefficient (cm2 s-1)
 
@@ -1337,7 +1343,10 @@ def calc_spectrum_clouds(
             plt.ylim(1e3, 1e-6)
             plt.xlim(1e-10, 1.0)
             log_x_base_item = log_x_base[item]
-            fsed = cloud_dict["fsed"]
+            if "fsed" in cloud_dict:
+                fsed = cloud_dict["fsed"]
+            else:
+                fsed = cloud_dict[f"fsed_{item}(c)"]
             log_kzz = cloud_dict["log_kzz"]
             plt.title(
                 f"fsed = {fsed:.2f}, log(Kzz) = {log_kzz:.2f}, "