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

slsim/Deflectors/DeflectorPopulation/galaxy_deflectors.py

@@ -3,6 +3,7 @@
     vel_disp_abundance_matching,
 )
 from slsim.Deflectors.DeflectorPopulation.deflectors_base import DeflectorsBase
+from slsim.Util.color_gradient import attach_foreground_deflector_color_gradient
 from astropy.table import vstack
 
 
@@ -26,6 +27,8 @@ def __init__(
         catalog_type="skypy",
         mass_type="EPL",
         light_type="single_sersic",
+        foreground_color_gradient=None,
+        foreground_component_weights=(0.4, 0.6),
     ):
         """
         :param red_galaxy_list: list of dictionary with elliptical galaxy
@@ -50,13 +53,41 @@ def __init__(
         :type mass_type: string
         :param light_type: type of Source() model class for the light distribution
         :type light_type: string
+        :param foreground_color_gradient: Optional configuration for a
+         band-dependent two-component deflector light profile. When provided,
+         the light is handled by the standard ``Source``/``DoubleSersic`` path.
+        :type foreground_color_gradient: dict or None
+        :param foreground_component_weights: Reference-band flux weights of the
+         two Sersic components.
+        :type foreground_component_weights: tuple or list
         :param catalog_type: type of the catalog. If user is using deflector catalog
          other than generated from skypy pipeline, we require them to provide angular
          size of the galaxy in arcsec and specify catalog_type as None. Otherwise, by
          default, this class considers deflector catalog is generated using skypy
          pipeline.
         :type catalog_type: str. "skypy" or None.
         """
+        if foreground_color_gradient is not None:
+            if light_type not in ("single_sersic", "double_sersic"):
+                raise ValueError(
+                    "foreground_color_gradient requires a single_sersic or "
+                    "double_sersic light_type."
+                )
+            light_type = "double_sersic"
+            red_galaxy_list = red_galaxy_list.copy()
+            attach_foreground_deflector_color_gradient(
+                red_galaxy_list,
+                foreground_color_gradient,
+                component_weights=foreground_component_weights,
+            )
+            if blue_galaxy_list is not None:
+                blue_galaxy_list = blue_galaxy_list.copy()
+                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/deflector_util.py

@@ -121,7 +121,35 @@ def light2mass(
             )
         else:
             # scale light to mass ellipticity
-            e1_light, e2_light = kwargs_source["e1"], kwargs_source["e2"]
+            if "e1" in kwargs_source and "e2" in kwargs_source:
+                e1_light, e2_light = kwargs_source["e1"], kwargs_source["e2"]
+            elif all(
+                key in kwargs_source
+                for key in ("e1_0", "e2_0", "e1_1", "e2_1", "w0", "w1")
+            ):
+                # Use the total reference-band light shape rather than an
+                # arbitrarily selected DoubleSersic component.
+                weights = np.asarray(
+                    [kwargs_source["w0"], kwargs_source["w1"]], dtype=float
+                )
+                if np.any(weights < 0) or np.sum(weights) <= 0:
+                    raise ValueError(
+                        "DoubleSersic weights 'w0' and 'w1' must be non-negative "
+                        "and have a positive sum when deriving mass ellipticity."
+                    )
+                weights /= np.sum(weights)
+                e1_light = np.dot(
+                    weights, [kwargs_source["e1_0"], kwargs_source["e1_1"]]
+                )
+                e2_light = np.dot(
+                    weights, [kwargs_source["e2_0"], kwargs_source["e2_1"]]
+                )
+            else:
+                raise ValueError(
+                    "Cannot derive mass ellipticity from the light model. Provide "
+                    "either ('e1', 'e2') or all DoubleSersic fields "
+                    "('e1_0', 'e2_0', 'e1_1', 'e2_1', 'w0', 'w1')."
+                )
             e1_mass, e2_mass = (
                 light2mass_e_scaling * e1_light,
                 light2mass_e_scaling * e2_light,

slsim/ImageSimulation/image_quality_lenstronomy.py

@@ -1,7 +1,9 @@
+import numpy as np
+import speclite.filters
+
 from lenstronomy.SimulationAPI.ObservationConfig.LSST import LSST
 from lenstronomy.SimulationAPI.ObservationConfig.Roman import Roman
 from lenstronomy.SimulationAPI.ObservationConfig.Euclid import Euclid
-import speclite.filters
 
 _OBSERVATORY_REGISTRY = {}
 
@@ -10,6 +12,10 @@
 LSST_BAND_LIST = ["u", "g", "r", "i", "z", "y"]
 EUCLID_BAND_LIST = ["VIS", "Y", "J", "H"]
 
+# Ancillary bandpasses that are used by catalog sources but are not imaging
+# observatories registered below.
+_ADDITIONAL_BAND_EFFECTIVE_WAVELENGTH_MICRON = {"F814W": 0.805}
+
 
 def check_speclite_name(band):
     """Checks if the raw band name is a valid speclite filter.
@@ -33,6 +39,7 @@ def register_observatory(
     bands: list,
     speclite_fmt=check_speclite_name,
     sncosmo_fmt=None,
+    effective_wavelengths=None,
 ):
     """Register a new observatory to integrate it with image simulation tools.
 
@@ -56,6 +63,10 @@ def register_observatory(
     :param sncosmo_fmt: A callable function that takes a ``band`` string and returns the corresponding
         sncosmo bandpass name. Set to ``None`` to use the raw band name as the sncosmo bandpass name.
     :type sncosmo_fmt: callable, optional
+    :param effective_wavelengths: Optional mapping from registered band names to
+        throughput-weighted effective wavelengths in microns. Use this only
+        when the responses are not available through ``speclite_fmt``.
+    :type effective_wavelengths: dict or None
 
     Given below is a simple example of how to define a custom observatory and register it using this function.
     A sophisticated example demonstrating full image simulation capabilities can be found at https://github.com/timedilatesme/MidEx-sims/blob/main/v1/lagn_sims.ipynb
@@ -106,6 +117,7 @@ def kwargs_single_band(self):
         "bands": list(bands),
         "speclite_fmt": speclite_fmt,
         "sncosmo_fmt": sncosmo_fmt,
+        "effective_wavelengths": dict(effective_wavelengths or {}),
     }
 
 
@@ -123,6 +135,16 @@ def kwargs_single_band(self):
     bands=ROMAN_BAND_LIST,
     speclite_fmt=lambda band: f"Roman-{band}",
     sncosmo_fmt=lambda band: f"{band}",
+    effective_wavelengths={
+        "F062": 0.620,
+        "F087": 0.870,
+        "F106": 1.060,
+        "F129": 1.290,
+        "F146": 1.460,
+        "F158": 1.580,
+        "F184": 1.840,
+        "F213": 2.130,
+    },
 )
 register_observatory(
     name="Euclid",
@@ -256,3 +278,54 @@ def get_all_supported_bands():
     for info in _OBSERVATORY_REGISTRY.values():
         all_bands.extend(info["bands"])
     return all_bands
+
+
+def get_band_effective_wavelength(band):
+    """Return the throughput-weighted effective wavelength of a band.
+
+    The registered speclite response is used when available. Roman and
+    HST bands that are not shipped by speclite use explicitly configured
+    fallback values. Band-list position is deliberately not used as a
+    wavelength proxy.
+
+    :param band: Imaging band name.
+    :type band: str
+    :return: Effective wavelength in microns.
+    :rtype: float
+    :raises ValueError: if the band is not registered.
+    """
+    if band in _ADDITIONAL_BAND_EFFECTIVE_WAVELENGTH_MICRON:
+        return _ADDITIONAL_BAND_EFFECTIVE_WAVELENGTH_MICRON[band]
+
+    obs_name = get_observatory(band)
+    configured_wavelengths = _OBSERVATORY_REGISTRY[obs_name]["effective_wavelengths"]
+    if band in configured_wavelengths:
+        return float(configured_wavelengths[band])
+
+    filter_name = get_speclite_filtername(band)
+    if filter_name is None:
+        raise ValueError(
+            f"Band '{band}' has neither an explicitly configured effective "
+            f"wavelength nor a registered speclite filter response."
+        )
+    response = speclite.filters.load_filter(filter_name)
+    return float(response.effective_wavelength.to("micron").value)
+
+
+def get_band_central_wavelength(band):
+    """Backward-compatible alias for :func:`get_band_effective_wavelength`."""
+    return get_band_effective_wavelength(band)
+
+
+def get_band_log_wavelength_ratio(band, reference_band):
+    """Return ``log(lambda_band / lambda_reference)``.
+
+    Both wavelengths are throughput-weighted effective wavelengths. This
+    is the chromatic coordinate used by the local power-law SED
+    approximation.
+    """
+    wavelength = get_band_effective_wavelength(band)
+    reference_wavelength = get_band_effective_wavelength(reference_band)
+    if wavelength <= 0 or reference_wavelength <= 0:
+        raise ValueError("Band effective wavelengths must be positive.")
+    return float(np.log(wavelength / reference_wavelength))

slsim/Lenses/lens.py

@@ -1388,8 +1388,8 @@ def lenstronomy_kwargs(
             field_galaxies_lens_model_list, kwargs_field_galaxies = (
                 self.field_galaxy_light_model_lenstronomy(band=band)
             )
-            lens_light_model_list += field_galaxies_lens_model_list
-            kwargs_lens_light += kwargs_field_galaxies
+            lens_light_model_list.extend(field_galaxies_lens_model_list)
+            kwargs_lens_light.extend(kwargs_field_galaxies)
 
         kwargs_model = {
             "lens_light_model_list": lens_light_model_list,