Add band-dependent Double-sersic/Catalog-Source color gradients

slsim/Deflectors/DeflectorPopulation/all_lens_galaxies.py

@@ -11,6 +11,7 @@
 from astropy.table import vstack
 from slsim.Util.param_util import catalog_with_angular_size_in_arcsec
 from slsim.Deflectors.deflector import Deflector
+from slsim.Util.color_gradient import attach_foreground_deflector_color_gradient
 
 
 class AllLensGalaxies(DeflectorsBase):
@@ -26,6 +27,8 @@ def __init__(
         sky_area,
         gamma_pl=None,
         catalog_type="skypy",
+        foreground_color_gradient=None,
+        foreground_component_weights=(0.4, 0.6),
     ):
         """
         :param red_galaxy_list: list of dictionary with elliptical galaxy
@@ -52,13 +55,34 @@ def __init__(
          default, this class considers deflector catalog is generated using skypy
          pipeline.
         :type catalog_type: str. "skypy" or None.
+        :param foreground_color_gradient: Optional dictionary defining a
+         band-dependent foreground deflector light colour gradient. When provided,
+         the deflector light model receives the corresponding ``color_gradient``
+         plus two reference-band component weights.
+        :type foreground_color_gradient: dict or None
+        :param foreground_component_weights: Two reference-band flux weights for
+         the foreground light components.
+        :type foreground_component_weights: tuple or list
         """
         red_galaxy_list = catalog_with_angular_size_in_arcsec(
             galaxy_catalog=red_galaxy_list, input_catalog_type=catalog_type
         )
         blue_galaxy_list = catalog_with_angular_size_in_arcsec(
             galaxy_catalog=blue_galaxy_list, input_catalog_type=catalog_type
         )
+        if foreground_color_gradient is not None:
+            red_galaxy_list = red_galaxy_list.copy()
+            blue_galaxy_list = blue_galaxy_list.copy()
+            attach_foreground_deflector_color_gradient(
+                red_galaxy_list,
+                foreground_color_gradient,
+                component_weights=foreground_component_weights,
+            )
+            attach_foreground_deflector_color_gradient(
+                blue_galaxy_list,
+                foreground_color_gradient,
+                component_weights=foreground_component_weights,
+            )
         red_column_names = red_galaxy_list.colnames
         if "galaxy_type" not in red_column_names:
             red_galaxy_list["galaxy_type"] = "red"

slsim/Deflectors/DeflectorPopulation/elliptical_lens_galaxies.py

@@ -7,6 +7,7 @@
     vel_disp_abundance_matching,
 )
 from slsim.Deflectors.deflector import Deflector
+from slsim.Util.color_gradient import attach_foreground_deflector_color_gradient
 
 
 class EllipticalLensGalaxies(DeflectorsBase):
@@ -21,6 +22,8 @@ def __init__(
         sky_area,
         gamma_pl=None,
         catalog_type="skypy",
+        foreground_color_gradient=None,
+        foreground_component_weights=(0.4, 0.6),
     ):
         """
 
@@ -44,10 +47,25 @@ def __init__(
          gamma for uniform distribution. eg: gamma_pl=2.1, gamma_pl={"mean": a, "std_dev": b},
          gamma_pl={"gamma_min": c, "gamma_max": d}
         :type catalog_type: str. "skypy" or None.
+        :param foreground_color_gradient: Optional dictionary defining a
+         band-dependent foreground deflector light colour gradient. When provided,
+         the deflector light model receives the corresponding ``color_gradient``
+         plus two reference-band component weights.
+        :type foreground_color_gradient: dict or None
+        :param foreground_component_weights: Two reference-band flux weights for
+         the foreground light components.
+        :type foreground_component_weights: tuple or list
         """
         galaxy_list = param_util.catalog_with_angular_size_in_arcsec(
             galaxy_catalog=galaxy_list, input_catalog_type=catalog_type
         )
+        if foreground_color_gradient is not None:
+            galaxy_list = galaxy_list.copy()
+            attach_foreground_deflector_color_gradient(
+                galaxy_list,
+                foreground_color_gradient,
+                component_weights=foreground_component_weights,
+            )
         super().__init__(
             deflector_table=galaxy_list,
             kwargs_cut=kwargs_cut,

slsim/Deflectors/DeflectorTypes/epl_sersic.py

@@ -1,5 +1,7 @@
 from slsim.Deflectors.DeflectorTypes.epl import EPL
 from slsim.Util.param_util import ellipticity_slsim_to_lenstronomy
+from slsim.Util.color_gradient import component_weights_for_band
+import numpy as np
 
 
 class EPLSersic(EPL):
@@ -58,6 +60,15 @@ def light_model_lenstronomy(self, band=None):
             )
         )
         size_lens_arcsec = self.angular_size_light
+        if self._has_band_dependent_color_gradient():
+            return self._double_sersic_light_model_lenstronomy(
+                band=band,
+                mag_lens=mag_lens,
+                size_lens_arcsec=size_lens_arcsec,
+                e1_light_lens_lenstronomy=e1_light_lens_lenstronomy,
+                e2_light_lens_lenstronomy=e2_light_lens_lenstronomy,
+                center_lens=center_lens,
+            )
         lens_light_model_list = ["SERSIC_ELLIPSE"]
         kwargs_lens_light = [
             {
@@ -71,3 +82,104 @@ def light_model_lenstronomy(self, band=None):
             }
         ]
         return lens_light_model_list, kwargs_lens_light
+
+    def _has_band_dependent_color_gradient(self):
+        """Return whether the deflector light should use chromatic
+        components."""
+        color_gradient = self._deflector_dict.get("color_gradient")
+        return (
+            isinstance(color_gradient, dict)
+            and color_gradient.get("component_spectral_slopes") is not None
+        )
+
+    def _double_sersic_light_model_lenstronomy(
+        self,
+        band,
+        mag_lens,
+        size_lens_arcsec,
+        e1_light_lens_lenstronomy,
+        e2_light_lens_lenstronomy,
+        center_lens,
+    ):
+        """Return two Sersic components with band-dependent component
+        weights."""
+        flux = 10 ** (-mag_lens / 2.5)
+        w0, w1 = self._weights_for_band(band)
+        mag_lens0 = -2.5 * np.log10(w0 * flux)
+        mag_lens1 = -2.5 * np.log10(w1 * flux)
+        angular_size_0, angular_size_1 = self._component_angular_sizes(size_lens_arcsec)
+        n_sersic_0, n_sersic_1 = self._component_sersic_indices()
+
+        lens_light_model_list = ["SERSIC_ELLIPSE", "SERSIC_ELLIPSE"]
+        kwargs_lens_light = [
+            {
+                "magnitude": mag_lens0,
+                "R_sersic": angular_size_0,
+                "n_sersic": n_sersic_0,
+                "e1": e1_light_lens_lenstronomy,
+                "e2": e2_light_lens_lenstronomy,
+                "center_x": center_lens[0],
+                "center_y": center_lens[1],
+            },
+            {
+                "magnitude": mag_lens1,
+                "R_sersic": angular_size_1,
+                "n_sersic": n_sersic_1,
+                "e1": e1_light_lens_lenstronomy,
+                "e2": e2_light_lens_lenstronomy,
+                "center_x": center_lens[0],
+                "center_y": center_lens[1],
+            },
+        ]
+        return lens_light_model_list, kwargs_lens_light
+
+    def _weights_for_band(self, band):
+        """Return Sersic component weights for an imaging band."""
+        w0 = float(self._deflector_dict.get("w0", 0.5))
+        w1 = self._deflector_dict.get("w1")
+        if w1 is None:
+            w1 = 1 - w0
+        return component_weights_for_band(
+            base_weights=(w0, float(w1)),
+            band=band,
+            color_gradient=self._deflector_dict.get("color_gradient"),
+            source_dict=self._deflector_dict,
+            default_reference="i",
+        )
+
+    def _component_angular_sizes(self, size_lens_arcsec):
+        """Return component half-light radii, using conservative opt-in
+        defaults."""
+        angular_size_0 = self._deflector_dict.get("angular_size_0")
+        angular_size_1 = self._deflector_dict.get("angular_size_1")
+        if angular_size_0 is not None and angular_size_1 is not None:
+            return float(angular_size_0), float(angular_size_1)
+
+        radius_factors = self._deflector_dict.get("color_gradient", {}).get(
+            "component_radius_factors", (0.5, 1.5)
+        )
+        if len(radius_factors) != 2:
+            raise ValueError(
+                "color_gradient['component_radius_factors'] must contain two values."
+            )
+        return (
+            float(size_lens_arcsec) * float(radius_factors[0]),
+            float(size_lens_arcsec) * float(radius_factors[1]),
+        )
+
+    def _component_sersic_indices(self):
+        """Return component Sersic indices for the chromatic light model."""
+        n_sersic_0 = self._deflector_dict.get("n_sersic_0")
+        n_sersic_1 = self._deflector_dict.get("n_sersic_1")
+        if n_sersic_0 is not None and n_sersic_1 is not None:
+            return float(n_sersic_0), float(n_sersic_1)
+
+        n_sersic = float(self._deflector_dict["n_sersic"])
+        indices = self._deflector_dict.get("color_gradient", {}).get(
+            "component_sersic_indices", (n_sersic, n_sersic)
+        )
+        if len(indices) != 2:
+            raise ValueError(
+                "color_gradient['component_sersic_indices'] must contain two values."
+            )
+        return float(indices[0]), float(indices[1])

slsim/ImageSimulation/image_quality_lenstronomy.py

@@ -10,6 +10,28 @@
 LSST_BAND_LIST = ["u", "g", "r", "i", "z", "y"]
 EUCLID_BAND_LIST = ["VIS", "Y", "J", "H"]
 
+_DEFAULT_BAND_CENTRAL_WAVELENGTH_MICRON = {
+    "u": 0.367,
+    "g": 0.482,
+    "r": 0.622,
+    "i": 0.755,
+    "z": 0.869,
+    "y": 0.971,
+    "VIS": 0.715,
+    "Y": 1.063,
+    "J": 1.285,
+    "H": 1.577,
+    "F062": 0.620,
+    "F087": 0.870,
+    "F106": 1.060,
+    "F129": 1.290,
+    "F146": 1.460,
+    "F158": 1.580,
+    "F184": 1.840,
+    "F213": 2.130,
+    "F814W": 0.805,
+}
+
 
 def check_speclite_name(band):
     """Checks if the raw band name is a valid speclite filter.
@@ -256,3 +278,59 @@ def get_all_supported_bands():
     for info in _OBSERVATORY_REGISTRY.values():
         all_bands.extend(info["bands"])
     return all_bands
+
+
+def get_band_central_wavelength(band):
+    """Return an approximate central wavelength for a registered band.
+
+    The built-in LSST, Roman, and Euclid bands use fixed wavelength
+    values in microns. For custom registered observatories, this
+    function falls back to the band's order within the observatory
+    registry so callers can still build monotonic band-dependent
+    behavior.
+
+    :param band: Imaging band name.
+    :type band: str
+    :return: Approximate central wavelength in microns.
+    :rtype: float
+    :raises ValueError: if the band is not registered.
+    """
+    if band in _DEFAULT_BAND_CENTRAL_WAVELENGTH_MICRON:
+        return _DEFAULT_BAND_CENTRAL_WAVELENGTH_MICRON[band]
+
+    obs_name = get_observatory(band)
+    bands = _OBSERVATORY_REGISTRY[obs_name]["bands"]
+    if len(bands) == 1:
+        return 0.0
+    return float(bands.index(band)) / float(len(bands) - 1)
+
+
+def get_band_normalized_position(band, reference_band=None):
+    """Return a normalized wavelength position for a registered band.
+
+    Built-in LSST, Roman, and Euclid bands are normalized over the full set of
+    built-in bands, so mixed-observatory simulations are ordered consistently.
+    Custom bands fall back to their observatory-local registry order.
+
+    :param band: Imaging band name.
+    :type band: str
+    :param reference_band: Optional reference band. If provided, return
+        ``position(band) - position(reference_band)``.
+    :type reference_band: str or None
+    :return: Normalized position or position offset.
+    :rtype: float
+    :raises ValueError: if any requested band is not registered.
+    """
+    wavelength = get_band_central_wavelength(band)
+
+    if band in _DEFAULT_BAND_CENTRAL_WAVELENGTH_MICRON:
+        wavelength_values = list(_DEFAULT_BAND_CENTRAL_WAVELENGTH_MICRON.values())
+        min_wavelength = min(wavelength_values)
+        max_wavelength = max(wavelength_values)
+        position = (wavelength - min_wavelength) / (max_wavelength - min_wavelength)
+    else:
+        position = wavelength
+
+    if reference_band is None:
+        return position
+    return position - get_band_normalized_position(reference_band)