Add kilonova event source type

slsim/Sources/Events/BNSMerger/kilonova.py

@@ -0,0 +1,169 @@
+import numpy as np
+
+if not hasattr(np, "trapezoid"):
+    np.trapezoid = np.trapz
+
+from astropy import cosmology
+from redback.transient_models import kilonova_models
+
+
+class Kilonova:
+    """Class for initializing a kilonova light curve model.
+
+    If modeldir is provided, external kilonova model files are expected.
+    This option is currently not supported. If modeldir is not provided,
+    the model is retrieved from Redback's built-in kilonova models. By
+    default, the MOSFiT-based kilonova model is used. Information about
+    Redback can be found at
+    https://redback.readthedocs.io/en/latest/.
+
+    Following the GW170817 three-component convention, components 1, 2,
+    and 3 can represent the blue, purple, and red ejecta components when
+    assigned low, intermediate, and high opacities, respectively. This
+    correspondence is not enforced, and users may specify the parameters
+    of each component independently.
+    """
+
+    def __init__(
+        self,
+        redshift,
+        mej_1,
+        mej_2,
+        mej_3,
+        vej_1,
+        vej_2,
+        vej_3,
+        kappa_1,
+        kappa_2,
+        kappa_3,
+        temperature_floor_1,
+        temperature_floor_2,
+        temperature_floor_3,
+        model_name="mosfit_kilonova",
+        kappa_gamma=10,
+        mag_zpsys="AB",
+        cosmo=cosmology.FlatLambdaCDM(H0=70, Om0=0.3),
+        modeldir=None,
+        **kwargs,
+    ):
+        """
+        :param redshift: The redshift of the kilonova source.
+        :type redshift: float
+
+        :param mej_1: Ejecta mass of model component 1 in [M_sun], which
+            sets the amount of radiating material and affecting the photon diffusion
+            timescale.
+        :type mej_1: float
+        :param mej_2: Ejecta mass of model component 2, with the same
+            definition and units as ``mej_1``.
+        :type mej_2: float
+        :param mej_3: Ejecta mass of model component 3, with the same
+            definition and units as ``mej_1``.
+        :type mej_3: float
+        :param vej_1: Expansion velocity of model component 1 in units of the
+            speed of light [c], affecting its expansion and photon diffusion
+            timescale.
+        :type vej_1: float
+        :param vej_2: Expansion velocity of model component 2, with the same
+            definition and units as ``vej_1``.
+        :type vej_2: float
+        :param vej_3: Expansion velocity of model component 3, with the same
+            definition and units as ``vej_1``.
+        :type vej_3: float
+        :param kappa_1: Effective gray opacity of model component 1 in
+            [cm^2 g^-1], controlling how readily radiation escapes from the
+            ejecta.
+        :type kappa_1: float
+        :param kappa_2: Effective gray opacity of model component 2, with the
+            same definition and units as ``kappa_1``.
+        :type kappa_2: float
+        :param kappa_3: Effective gray opacity of model component 3, with the
+            same definition and units as ``kappa_1``.
+        :type kappa_3: float
+        :param temperature_floor_1: Minimum effective photospheric temperature of model
+            component 1 in [K], affecting its late-time spectral evolution.
+        :type temperature_floor_1: float
+        :param temperature_floor_2: Minimum effective photospheric temperature of model
+            component 2, with the same definition and units as ``temperature_floor_1``.
+        :type temperature_floor_2: float
+        :param temperature_floor_3: Minimum effective photospheric temperature of model
+            component 3, with the same definition and units as ``temperature_floor_1``.
+        :type temperature_floor_3: float
+
+        :param model_name: The kilonova light curve model to be used. If not provided,
+            the default model is the MOSFiT-based kilonova model.
+        :type model_name: str
+        :param kappa_gamma: Gamma-ray opacity shared by all three model components
+            in [cm^2 g^-1], controlling the trapping and escape of high-energy radiation
+            from radioactive decay.
+        :type kappa_gamma: float
+        :param mag_zpsys: Optional, AB or Vega (AB default).
+        :type mag_zpsys: str
+        :param cosmo: Cosmology for luminosity distance calculation.
+        :type cosmo: `~astropy.cosmology`
+        :param modeldir: Directory including files for external kilonova models.
+        :type modeldir: str or None
+        :param kwargs: Additional keyword arguments passed to the Redback kilonova model.
+        :type kwargs: dict
+        """
+
+        if modeldir is not None:
+            raise NotImplementedError(
+                "Only built-in Redback kilonova models are currently supported. "
+                "External kilonova model files are not supported."
+            )
+
+        if not hasattr(kilonova_models, model_name):
+            raise ValueError(
+                f"Unsupported kilonova model '{model_name}'. "
+                "The model must be available in "
+                "redback.transient_models.kilonova_models."
+            )
+
+        self._model = getattr(kilonova_models, model_name)
+
+        self._model_name = model_name
+        self._redshift = redshift
+        self._mag_zpsys = mag_zpsys
+        self._cosmo = cosmo
+        self._kwargs = kwargs
+
+        self._model_parameters = {
+            "mej_1": mej_1,
+            "mej_2": mej_2,
+            "mej_3": mej_3,
+            "vej_1": vej_1,
+            "vej_2": vej_2,
+            "vej_3": vej_3,
+            "kappa_1": kappa_1,
+            "kappa_2": kappa_2,
+            "kappa_3": kappa_3,
+            "temperature_floor_1": temperature_floor_1,
+            "temperature_floor_2": temperature_floor_2,
+            "temperature_floor_3": temperature_floor_3,
+            "kappa_gamma": kappa_gamma,
+        }
+
+    def get_apparent_magnitude(self, time, band, zpsys="AB"):
+        """Function to return apparent magnitude of a kilonova for a given band
+        and time.
+
+        :param time: The observer-frame time array to evaluate the model
+            (in days)
+        :type time: array-like
+        :param band: The band to evaluate the model over.
+        :type band: str or list
+        :param zpsys: Optional, AB or Vega (AB default)
+        :type zpsys: str
+        :return: magnitude of source
+        """
+
+        return self._model(
+            time=time,
+            redshift=self._redshift,
+            bands=band,
+            output_format="magnitude",
+            cosmology=self._cosmo,
+            **self._model_parameters,
+            **self._kwargs,
+        )

slsim/Sources/SourceTypes/kilonova_event.py

@@ -0,0 +1,158 @@
+import warnings
+import numpy as np
+from slsim.Sources.SourceTypes.source_base import SourceBase
+from slsim.Sources.Events.BNSMerger.kilonova import Kilonova
+from slsim.ImageSimulation.image_quality_lenstronomy import (
+    get_all_supported_bands,
+    get_sncosmo_filtername,
+)
+
+
+class KilonovaEvent(SourceBase):
+    """A class to manage a BNS merger."""
+
+    def __init__(
+        self,
+        lightcurve_time,
+        variability_model,
+        model_name="mosfit_kilonova",
+        mag_zpsys="AB",
+        modeldir=None,
+        kwargs_variability=None,
+        kwargs_kilonova=None,
+        cosmo=None,
+        **kwargs,
+    ):
+        """
+        :param lightcurve_time: Observation time array for the light curve in [days].
+        :type lightcurve_time: array-like
+        :param variability_model: Keyword for the variability model to be used. This is an
+            input for the Variability class.
+        :type variability_model: str
+        :param model_name: Kilonova light curve model to be used. If not provided, the
+            default model is the MOSFiT-based kilonova model.
+        :type model_name: str
+        :param mag_zpsys: Optional, AB or Vega (AB default).
+        :type mag_zpsys: str
+        :param modeldir: Directory including files for external kilonova models. This
+            option is currently not supported.
+        :type modeldir: str or None
+        :param kwargs_variability: List containing the variability keyword and the
+            bands for which the light curve should be generated.
+        :type kwargs_variability: list or None
+        :param kwargs_kilonova: Keyword arguments passed to the Kilonova class. This may
+            include ejecta_mass in [solar masses], ejecta_velocity in [c], opacity in
+            [cm^2 g^-1], temperature_floor in [K], and kappa_gamma in [cm^2 g^-1].
+        :type kwargs_kilonova: dict or None
+        :param cosmo: Astropy cosmology instance.
+        :type cosmo: `~astropy.cosmology`
+        :param kwargs: Keyword arguments passed to the SourceBase class. This may contain
+            source properties such as redshift and offsets from the host galaxy in [arcsec].
+        :type kwargs: dict
+        """
+        super().__init__(
+            extended_source=False,
+            point_source=True,
+            cosmo=cosmo,
+            variability_model=variability_model,
+            **kwargs,
+        )
+        self.name = "BNS"
+        self._variability_computed = False
+        self._kwargs_variability = kwargs_variability
+        self._lightcurve_time = lightcurve_time
+
+        self._model_name = model_name
+        self._mag_zpsys = mag_zpsys
+        self._modeldir = modeldir
+
+        if kwargs_kilonova is None:
+            self._kwargs_kilonova = {}
+        else:
+            self._kwargs_kilonova = kwargs_kilonova
+
+    @property
+    def light_curve(self):
+        """Provides lightcurves of a bns merger in each band."""
+        if self._kwargs_variability is not None:
+            kwargs_variab_extracted = {}
+            if self._cosmo is None:
+                raise ValueError(
+                    "Cosmology cannot be None for BNSMerger class. Please"
+                    "provide a suitable astropy cosmology."
+                )
+            else:
+                # Initialize BNS/Kilonova light curve model
+                lightcurve_class = Kilonova(
+                    redshift=self._z,
+                    model_name=self._model_name,
+                    mag_zpsys=self._mag_zpsys,
+                    cosmo=self._cosmo,
+                    modeldir=self._modeldir,
+                    **self._kwargs_kilonova,
+                )
+                self._lightcurve_class = lightcurve_class
+
+            supported_bands = get_all_supported_bands()
+            provided_bands = set(supported_bands) & set(self._kwargs_variability)
+
+            for element in provided_bands:
+                name = "ps_mag_" + element
+                times = self._lightcurve_time
+
+                # Convert SLSim short band labels, such as "i" and "r", to registered
+                # filter names required by Redback, such as "lssti" and "lsstr".
+                # We reuse the existing SLSim/SNcosmo filter-name helper for this mapping.
+                provided_band = get_sncosmo_filtername(element)
+
+                try:
+                    magnitudes = lightcurve_class.get_apparent_magnitude(
+                        time=times,
+                        band=provided_band,
+                        zpsys=self._mag_zpsys,
+                    )
+                    # make sure before and after the event, the flux is zero
+                    magnitudes = np.append(np.inf, magnitudes)
+                    magnitudes = np.append(magnitudes, np.inf)
+                    padded_times = np.append(times[0] - (times[1] - times[0]), times)
+                    padded_times = np.append(padded_times, 2 * times[-1] - times[-2])
+                except Exception as e:
+                    warnings.warn(
+                        f"Skipping band '{provided_band}': Failed to generate lightcurve. "
+                        f"(Error: {e})",
+                        UserWarning,
+                    )
+                    continue
+
+                if name not in self.source_dict:
+                    self.source_dict[name] = float(min(magnitudes))
+
+                kwargs_variab_extracted[element] = {
+                    "MJD": padded_times,
+                    name: magnitudes,
+                }
+        else:
+            kwargs_variab_extracted = {}
+
+        self._variability_computed = True
+        return kwargs_variab_extracted
+
+    def point_source_magnitude(self, band, image_observation_times=None):
+        """Get the magnitude of the BNS/kilonova point source in a specific
+        band.
+
+        :param band: Imaging band.
+        :type band: str
+        :param image_observation_times: Image observation times in
+            [days]. If None, takes the peak magnitude.
+        :type image_observation_times: array-like or None
+        :return: Magnitude of the point source in the specified band.
+        :rtype: float or array-like
+        """
+
+        if not self._variability_computed:
+            self._kwargs_variability_model = self.light_curve
+
+        return super().point_source_magnitude(
+            band=band, image_observation_times=image_observation_times
+        )

slsim/Sources/SourceTypes/point_source.py

@@ -1,4 +1,4 @@
-_SUPPORTED_POINT_SOURCES = ["supernova", "quasar", "general_lightcurve"]
+_SUPPORTED_POINT_SOURCES = ["supernova", "quasar", "general_lightcurve", "kilonova"]
 
 
 class PointSource(object):
@@ -25,6 +25,10 @@ def __init__(self, source_type, **source_dict):
             from slsim.Sources.SourceTypes.general_lightcurve import GeneralLightCurve
 
             self._point_source = GeneralLightCurve(**source_dict)
+        elif source_type in ["kilonova"]:
+            from slsim.Sources.SourceTypes.kilonova_event import KilonovaEvent
+
+            self._point_source = KilonovaEvent(**source_dict)
         else:
             raise ValueError(
                 "Point source type %s not supported. Chose among %s."

slsim/Sources/source.py

@@ -2,7 +2,7 @@
 from slsim.Sources.SourceTypes.source_base import SourceBase
 from copy import deepcopy
 
-_SUPPORTED_POINT_SOURCES = ["supernova", "quasar", "general_lightcurve"]
+_SUPPORTED_POINT_SOURCES = ["supernova", "quasar", "general_lightcurve", "kilonova"]
 _SUPPORTED_EXTENDED_SOURCES = [
     "single_sersic",
     "double_sersic",
@@ -27,7 +27,7 @@ def __init__(
          extended source types are 'single_sersic', 'double_sersic', 'catalog_source', and 'interpolated'.
         :type extended_source_type: str or None
         :param point_source_type: Keyword to specify type of point source. Supported point
-         source types are 'supernova', 'quasar', and 'general_lightcurve'.
+         source types are 'supernova', 'quasar', 'general_lightcurve', and 'kilonova'.
         :type point_source_type: str or None
         :param source_dict: Source properties. Can be a dictionary or an Astropy table.
          For a detailed description of this dictionary, please see the documentation for
@@ -69,6 +69,10 @@ def __init__(
             from slsim.Sources.SourceTypes.general_lightcurve import GeneralLightCurve
 
             self._source = GeneralLightCurve(**source_dict)
+        elif source_type in ["kilonova"]:
+            from slsim.Sources.SourceTypes.kilonova_event import KilonovaEvent
+
+            self._source = KilonovaEvent(**source_dict)
 
         # extended sources
         elif source_type in ["single_sersic"]:

test_requirements.txt

@@ -7,3 +7,4 @@ pyhalo<=1.4.6
 mcfit
 sncosmo
 stpsf
+redback>=1.17.0