Added functionality to convert MJD time to days for transient events.

slsim/LsstSciencePipeline/opsim_pipeline.py

@@ -10,8 +10,8 @@ def opsim_time_series_images_data(
     dec_list,
     obs_strategy,
     MJD_min=60000,
-    MJD_max=64000,
-    size=101,
+    MJD_max=64500,
+    num_pix=101,
     moffat_beta=3.1,
     readout_noise=10,
     delta_pix=0.2,
@@ -28,7 +28,7 @@ def opsim_time_series_images_data(
         for example "baseline_v3.0_10yrs" (string)
     :param MJD_min: minimum MJD for the observations
     :param MJD_max: maximum MJD for the observations
-    :param size: cutout size of images (in pixels)
+    :param num_pix: cutout size of images (in pixels)
     :param moffat_beta: power index of the moffat psf kernel
     :param readout_noise: noise added per readout
     :param delta_pix: size of pixel in units arcseonds
@@ -111,7 +111,7 @@ def opsim_time_series_images_data(
 
         for psf in psf_fwhm:
             psf_kernel = kernel_util.kernel_moffat(
-                num_pix=size, delta_pix=delta_pix, fwhm=psf, moffat_beta=moffat_beta
+                num_pix=num_pix, delta_pix=delta_pix, fwhm=psf, moffat_beta=moffat_beta
             )
             psf_kernel = util.array2image(psf_kernel)
 

slsim/LsstSciencePipeline/util_lsst.py

@@ -1,7 +1,8 @@
 import numpy as np
 from astropy.table import Column
 from slsim.image_simulation import lens_image_series, lens_image
-from slsim.Util.param_util import fits_append_table
+from slsim.Util.param_util import (fits_append_table, convert_mjd_to_days,
+                                    transient_event_time_mjd)
 import os
 
 
@@ -28,17 +29,13 @@ def variable_lens_injection(
         in days for each single exposure images in time series images)
     :return: Astropy table of injected lenses and exposure information of dp0 data
     """
-    # the range of observation time of single exposure images might be outside of the
-    # lightcurve time. So, we use random observation time from the lens class lightcurve
-    # time to ensure simulation of reasonable images.
-    observation_time = np.random.uniform(
-        min(lens_class.max_redshift_source_class.lightcurve_time),
-        max(lens_class.max_redshift_source_class.lightcurve_time),
-        size=len(exposure_data["obs_time"]),
-    )
-    observation_time.sort()
-    new_obs_time = Column(name="obs_time", data=observation_time)
-    exposure_data.replace_column("obs_time", new_obs_time)
+    ## chose transient starting point randomly. 
+    start_point_mjd_time = transient_event_time_mjd(min(exposure_data["obs_time"]),
+                                        max(exposure_data["obs_time"]))
+    ## Convert mjd observation time to days. We should do this because lightcurves are
+    #  in the unit of days.
+    observation_time = convert_mjd_to_days(exposure_data["obs_time"],
+                                            start_point_mjd_time)
     lens_images = lens_image_series(
         lens_class,
         band=band,
@@ -47,7 +44,7 @@ def variable_lens_injection(
         psf_kernel=exposure_data["psf_kernel"],
         transform_pix2angle=transform_pix2angle,
         exposure_time=exposure_data["expo_time"],
-        t_obs=exposure_data["obs_time"],
+        t_obs=observation_time,
     )
 
     final_image = []
@@ -147,12 +144,16 @@ def opsim_variable_lens_injection(
     :return: Astropy table of injected lenses and exposure information of dp0 data
     """
 
+    ## chose transient starting point randomly. 
+    start_point_mjd_time = transient_event_time_mjd(min(exposure_data["obs_time"]),
+                                        max(exposure_data["obs_time"]))
     final_image = []
 
     for obs in range(len(exposure_data["obs_time"])):
 
         exposure_data_obs = exposure_data[obs]
-
+        observation_time = convert_mjd_to_days(exposure_data_obs["obs_time"],
+                                                start_point_mjd_time)
         if exposure_data_obs["band"] not in bands:
             continue
 
@@ -169,7 +170,7 @@ def opsim_variable_lens_injection(
             psf_kernel=exposure_data_obs["psf_kernel"],
             transform_pix2angle=transform_pix2angle,
             exposure_time=exposure_data_obs["expo_time"],
-            t_obs=exposure_data_obs["obs_time"],
+            t_obs=observation_time,
             std_gaussian_noise=std_gaussian_noise,
         )
 

slsim/Util/param_util.py

@@ -350,3 +350,26 @@ def catalog_with_angular_size_in_arcsec(galaxy_catalog, input_catalog_type="skyp
         )
         warnings.warn(warning_msg, category=UserWarning, stacklevel=2)
     return copied_galaxy_catalog
+
+def convert_mjd_to_days(reference_mjd, start_point_mjd):
+    """
+    Convert reference MJD(s) to days relative to a chosen zero-point MJD.
+
+    :param reference_mjd: The reference MJD(s) to convert. 
+    :type reference_mjd: float, list, or numpy.ndarray
+    :param start_point_mjd: The zero-point MJD to use as the reference.
+    :return: The time(s) in days relative to the zero-point MJD.
+    """
+    # Ensure input is a NumPy array for consistent handling
+    reference_mjd = np.array(reference_mjd)
+    return reference_mjd - start_point_mjd
+
+def transient_event_time_mjd(min_mjd, max_mjd):
+    """ Produces a random MJD time with in the given range
+
+    :param min_mjd: Minimum bound for the MJD time
+    :param max_mjd: Maximum bound for the MJD time
+    :return: A random MJD time between given min and max bounds.
+    """
+    start_mjd=np.random.randint(min_mjd, max_mjd)
+    return start_mjd

tests/test_param_util.py

@@ -16,6 +16,8 @@
     ellipticity_slsim_to_lenstronomy,
     fits_append_table,
     catalog_with_angular_size_in_arcsec,
+    convert_mjd_to_days,
+    transient_event_time_mjd,
 )
 from slsim.Sources.SourceVariability.variability import Variability
 from astropy.io import fits
@@ -250,6 +252,14 @@ def test_catalog_with_angular_size_in_arcsec():
     assert galaxy_cat2["angular_size"].unit == u.rad
     assert galaxy_cat["angular_size"].unit == u.arcsec
 
+def test_convert_mjd_to_days():
+    result = convert_mjd_to_days(60100, 60000)
+    assert result==100
+
+def test_start_point_mjd_time():
+    result = transient_event_time_mjd(60000, 60400)
+    assert 60000 <= result <=60400
+
 
 if __name__ == "__main__":
     pytest.main()