Supernovae plus extended source simulation

In this notebook, we simulate population of lensed supernovae and simulate image of a random lensed supernovae. It follows following steps:

1. Simulate lensed supernovae population

2. Choose a lens at random

3. Set observation time and other image configuration

4. Simulate image of a selected lens

5. Visualize it

Before running this notebook, please download the “scotch_SNIa_host_galaxies.fits”

file from the following link: https://github.com/LSST-strong-lensing/data_public.git.

This file contains type Ia supernovae host galaxies.

Simulate lensed supernovae population

import os
from astropy.cosmology import FlatLambdaCDM
from astropy.units import Quantity
from slsim.lens_pop import LensPop
import numpy as np
from slsim.image_simulation import lens_image_series
from slsim.Plots.plot_functions import create_image_montage_from_image_list
from slsim.image_simulation import point_source_coordinate_properties
import matplotlib.pyplot as plt
import corner

data_path = os.path.abspath("/Users/sanchez/Devel/DESC/data_public")

# define a cosmology
cosmo = FlatLambdaCDM(H0=70, Om0=0.3)


# define a sky area
sky_area = Quantity(value=1, unit="deg2")


# define limits in the intrinsic deflector and source population (in addition to the
# skypy config
# file)
kwargs_deflector_cut = {"z_min": 0.01, "z_max": 2.5}
kwargs_source_cut = {}
## create a point plus extended source lens population.
supernova_lens_pop = LensPop(
    deflector_type="elliptical",  # type of the deflector. It could be elliptical or
    # all-galaxies.
    source_type="supernovae_plus_galaxies",  # keyword for source type. it can be
    # galaxies, quasar, quasar_plus_galaxies, and supernovae_plus_galaxies.
    kwargs_deflector_cut=kwargs_deflector_cut,  # cuts that one wants to apply for the
    # deflector.
    kwargs_source_cut=kwargs_source_cut,  # cuts that one wants to apply for the
    # source.
    variability_model="light_curve",  # keyword for the variability model.
    kwargs_variability={"supernovae_lightcurve", "i"},  # specify kewords for
    # lightcurve. "i" is a band for the lightcurve.
    sn_type="Ia",  # supernovae type.
    sn_absolute_mag_band="bessellb",  # Band used to normalize to absolute magnitude
    sn_absolute_zpsys="ab",  # magnitude system. It can be Optional, AB or Vega.
    kwargs_mass2light=None,  # mass-to-light relation for the deflector galaxy.
    skypy_config=None,  # Sky configuration for the simulation. If None, lsst-like
    # configuration will be used.
    sky_area=sky_area,  # Sky area for the simulation
    cosmo=cosmo,  # astropy cosmology
    source_light_profile="double_sersic",  # light profile for the source galaxy
    catalog_type="scotch",  # catalog type. It can be None or scotch
    lightcurve_time=np.linspace(
        -20, 100, 1000
    ),  # array of light curve observation time.
    catalog_path=os.path.join(
        data_path, "SupernovaeHostcatalog/scotch_SNIa_host_galaxies.fits"
    ),
    # path for catalog. If not provided, small size catalog from
    # /slsim/Source/SupernovaeCatalog will be used for
    # source_type="supernovae_plus_galaxies" case. For other cases, we do not need to
    # provide outside catalog. One can download scotch_SNIa_host_galaxies.fits from
    # https://github.com/LSST-strong-lensing/data_public.git
)

Choose a random lens

# specifying cuts of the population
kwargs_lens_cuts = {}
# drawing population
supernovae_lens_population = supernova_lens_pop.draw_population(
    kwargs_lens_cuts=kwargs_lens_cuts
)

print("Number of lenses:", len(supernovae_lens_population))

lens_samples = []
labels = [
    r"$\sigma_v$",
    r"$\log(M_{*})$",
    r"$\theta_E$",
    r"$z_{\rm l}$",
    r"$z_{\rm s}$",
    r"$m_{\rm host}$",
    r"$m_{\rm ps}$",
    r"$m_{\rm lens}$",
]

for supernovae_lens in supernovae_lens_population:
    vel_disp = supernovae_lens.deflector_velocity_dispersion()
    m_star = supernovae_lens.deflector_stellar_mass()
    theta_e = supernovae_lens.einstein_radius
    zl = supernovae_lens.deflector_redshift
    zs = supernovae_lens.source_redshift
    source_mag = supernovae_lens.extended_source_magnitude(band="i", lensed=True)
    ps_source_mag = supernovae_lens.point_source_magnitude(band="i")
    deflector_mag = supernovae_lens.deflector_magnitude(band="i")
    lens_samples.append(
        [
            vel_disp,
            np.log10(m_star),
            theta_e,
            zl,
            source_mag,
            ps_source_mag,
            deflector_mag,
        ]
    )

Number of lenses: 26

/Users/sanchez/Devel/DESC/slsim/slsim/Sources/source.py:362: RuntimeWarning: divide by zero encountered in log10
  mag_source0 = -2.5 * np.log10(w0 * flux)
/Users/sanchez/.virtualenvs/slsim/lib/python3.12/site-packages/sncosmo/models.py:189: RuntimeWarning: divide by zero encountered in log10
  result[i] = -2.5 * np.log10(f / zpf)

hist2dkwargs = {
    "plot_density": False,
    "plot_contours": False,
    "plot_datapoints": True,
    "color": "b",
    "data_kwargs": {"ms": 5},
}
corner.corner(
    np.array(lens_samples), labels=labels, label_kwargs={"fontsize": 20}, **hist2dkwargs
)
plt.show()

Choose a lens to simulate an image

kwargs_lens_cut = {"min_image_separation": 1, "max_image_separation": 10}
rgb_band_list = ["i", "r", "g"]
lens_class = supernovae_lens_population[-5]
# (
#     lens_class.source.source_dict["z"],
#     lens_class.einstein_radius,
#     lens_class.source.source_dict["mag_i"],
#     lens_class.source.source_dict["ps_mag_i"],
#     lens_class._deflector_dict["mag_i"],
#     lens_class._deflector_dict["z"],
# )

pix_coord = point_source_coordinate_properties(
    lens_class,
    band="i",
    mag_zero_point=27,
    delta_pix=0.2,
    num_pix=32,
    transform_pix2angle=np.array([[0.2, 0], [0, 0.2]]),
)["image_pix"]

/Users/sanchez/Devel/DESC/slsim/slsim/Sources/source.py:362: RuntimeWarning: divide by zero encountered in log10
  mag_source0 = -2.5 * np.log10(w0 * flux)

pix_coord

array([[14.00699877, 16.4677447 ],
       [17.07628824, 15.42222337]])

See the light curve of a selected supernovae

light_curve = lens_class.source.variability_class.kwargs_model

plt.plot(light_curve["MJD"], light_curve["ps_mag_i"])
# plt.ylim(12, 18)
plt.gca().invert_yaxis()
plt.ylabel("Magnitude")
plt.xlabel("Time" "[Days]")
plt.xlim(-22, 100)

(-22.0, 100.0)

Set observation time and image configuration

time = np.array([-19.5, -15, -11.35135135135135, 0, 10, 20, 25, 30, 40, 44.86])
# time = sorted(np.random.uniform(-20, 100, 10))
# time = np.array([0, 50, 70, 120])
repeats = 10
# load your psf kernel and transform matrix. If you have your own psf, please provide
# it here.
path = "../../tests/TestData/psf_kernels_for_deflector.npy"
psf_kernel = 1 * np.load(path)
psf_kernel[psf_kernel < 0] = 0
transform_matrix = np.array([[0.2, 0], [0, 0.2]])

# let's set up psf kernel for each exposure. Here we have taken the same psf that we
# extracted above. However, each exposure can have different psf kernel and user should
# provide corresponding psf kernel to each exposure.
psf_kernel_list = [psf_kernel]
transform_matrix_list = [transform_matrix]
psf_kernels_all = psf_kernel_list * repeats
# psf_kernels_all = np.array([dp0["psf_kernel"][:10]])[0]

# let's set pixel to angle transform matrix. Here we have taken the same matrix for
# each exposure but user should provide corresponding transform matrix to each exposure.
transform_matrix_all = transform_matrix_list * repeats

# provide magnitude zero point for each exposures. Here we have taken the same magnitude
#  zero point for each exposure but user should provide the corresponding magnitude
# zero point for each exposure.
mag_list = [31.0]
mag_zero_points_all = mag_list * repeats
# mag_zero_points_all = np.array([dp0["zero_point"][:10]])[0]

expo_list = [30]
exposure_time_all = expo_list * repeats

Simulate Image

# Simulate a lens image
image_lens_series = lens_image_series(
    lens_class=lens_class,
    band="i",
    mag_zero_point=mag_zero_points_all,
    num_pix=32,
    psf_kernel=psf_kernels_all,
    transform_pix2angle=transform_matrix_all,
    exposure_time=exposure_time_all,
    t_obs=time,
    with_deflector=True,
    with_source=True,
)

/Users/sanchez/Devel/DESC/slsim/slsim/Sources/source.py:362: RuntimeWarning: divide by zero encountered in log10
  mag_source0 = -2.5 * np.log10(w0 * flux)

## Images in log scale
log_images = []
for i in range(len(image_lens_series)):
    log_images.append(np.log10(image_lens_series[i]))

/var/folders/61/jz5zv4qn1ml4kvb_mzhpz5500000gn/T/ipykernel_81151/422990630.py:4: RuntimeWarning: divide by zero encountered in log10
  log_images.append(np.log10(image_lens_series[i]))

Visualize simulated images

plot_montage = create_image_montage_from_image_list(
    num_rows=2, num_cols=5, images=image_lens_series, time=time, image_center=pix_coord
)