branch ris-custom-ramps contains edits to include CRNL and custom ramps

romanisim/image.py

@@ -23,8 +23,6 @@
 from romanisim import models
 from romanisim.models import wcs, parameters
 
-
-
 # galsim fluxes are in photons / cm^2 / s
 # we need to specify the area and exposure time in drawImage if
 # specifying fluxes in physical units.
@@ -894,6 +892,7 @@ def simulate(metadata, objlist,
              usecrds=True, psftype='galsim', level=2, crparam=dict(),
              persistence=None, seed=None, rng=None,
              psf_keywords=dict(), extra_counts=None,
+             fdnl_params=None, custom_ramps=None,
              **kwargs
              ):
     """Simulate a sequence of observations on a field in different bandpasses.
@@ -937,7 +936,22 @@ def simulate(metadata, objlist,
         An additional array that just gets added into the counts image.
         Useful for wrapping idealized images into L1/L2 images + the
         Roman datamodel.
+    fdnl_params : a 2-tuple or 2-element ndarray/list to be passed as follows to
+        galsim.detectors.addReciprocityFailure(exptime, *fdnl_params).  This
+        introduces flux-dependent nonlinearity (FDNL), a.k.a. Count Rate Dependent
+        Nonlinearity (CRNL) or Reciprocity Failure, to the image.  See
+        documentation of above Galsim method for further description of the
+        parameters.
+    custom_ramps : ndarray, an up-the-ramp sampled image cube that
+        will contribute additively to each read of pixel ramps
+        generated in l1.apportion_counts_to_resultants().  Ideally has
+        dtype='i4'.  The array shape must be (nr, 4088, 4088) where
+        nr>= the total number of reads (not resultants) of the ma
+        table specified in metadata.  If nr>nreads, the additional
+        frames will be ignored.  Useful for simulating ramp
+        nonlinearity effects such as the brighter fatter effect.
 
+
     Returns
     -------
     image : roman_datamodels model
@@ -946,8 +960,19 @@ def simulate(metadata, objlist,
         Dictionary of additionally valuable quantities.  Includes at least
         simcatobj, the image positions and fluxes of simulated objects.  It may
         also include information on persistence and cosmic ray hits.
+
     """
 
+    # populate keyword arguments if they are passed in through **kwargs, for example
+    # from romanisim.ris_make_utils.simulate_image_file()
+    if fdnl_params is None:
+        fdnl_params = kwargs.pop('fdnl_params', None)
+    if custom_ramps is None:
+        custom_ramps = kwargs.pop('custom_ramps', None)
+    if extra_counts is None:
+        extra_counts = kwargs.pop('extra_counts', None)
+
+
     if not usecrds:
         log.warning('--usecrds is not set.  romanisim will not use reference '
                     'files from CRDS.  The WCS may be incorrect and up-to-date '
@@ -998,12 +1023,20 @@ def simulate(metadata, objlist,
     log.info('Simulating filter {0}...'.format(filter_name))
     counts, simcatobj = simulate_counts(
         image_mod.meta, objlist, rng=rng, usecrds=usecrds, darkrate=darkrate,
-        psftype=psftype, flat=flat, gain=gain, psf_keywords=psf_keywords)
+        psftype=psftype, flat=flat, gain=gain, psf_keywords=psf_keywords
+    )
 
     # If extra_counts is passed in, add directly to counts
     if extra_counts is not None:
         counts += extra_counts
 
+    # add flux dependent nonlinearity
+    if fdnl_params is not None:
+        exptime = parameters.read_time * read_pattern[-1][-1]
+        counts.addReciprocityFailure(exptime, *fdnl_params)
+        counts.quantize()
+        log.info('Included Flux Dependent Nonlinearity')
+
     util.update_pointing_and_wcsinfo_metadata(image_mod.meta, counts.wcs)
     if level == 0:
         im = dict(data=counts.array, meta=dict(image_mod.meta.items()))
@@ -1019,6 +1052,7 @@ def simulate(metadata, objlist,
             persistence=persistence,
             saturation=saturation,
             darkdecaysignal=darkdecaysignal,
+            custom_ramps=custom_ramps,
             ipc_model=ipc_model,
             **kwargs)
     if level == 0:

romanisim/l1.py

@@ -268,7 +268,9 @@ def apply_dark_decay(resultants, darkdecaysignal, read_pattern, sign=1):
 def apportion_counts_to_resultants(
         counts, tij, pedestal=None, pedestal_extra_noise=None,
         inv_linearity=None, crparam=None, persistence=None,
-        tstart=None, rng=None, seed=None):
+        tstart=None, rng=None, seed=None,
+        custom_ramps=None
+        ):
     """Apportion counts to resultants given read times.
 
     This finds a statistically appropriate assignment of electrons to each
@@ -321,6 +323,14 @@ def apportion_counts_to_resultants(
         random number generator
     seed : int
         seed to use for random number generator
+    custom_ramps : ndarray, an up-the-ramp sampled image cube that
+        will contribute additively to each read of the pixel ramps
+        generated here.  Ideally has dtype='i4'.  The array shape must
+        be (nr, 4088, 4088) where nr>=nreads, the total number of
+        reads (not resultants) of the ma table specified in the
+        metadata of the observation.  If nr>nreads, the additional
+        frames will be ignored.  Useful for simulating ramp
+        nonlinearity effects such as the brighter-fatter effect.
 
     Returns
     -------
@@ -329,6 +339,7 @@ def apportion_counts_to_resultants(
         array of n_resultant images giving each resultant
     dq : np.ndarray[n_resultant, nx, ny]
         dq array marking CR hits in resultants
+
     """
     if not np.all(counts == np.round(counts)):
         raise ValueError('apportion_counts_to_resultants expects the counts '
@@ -382,6 +393,10 @@ def apportion_counts_to_resultants(
     if persistence is not None:
         tstart = tstart.mjd
 
+    # initialize custom_ramps counter
+    if custom_ramps is not None:
+        cr_count = 0
+
     # Loop over read probabilities
     for i, pi in enumerate(pij):
         # Reset resultant counts
@@ -393,6 +408,13 @@ def apportion_counts_to_resultants(
             read = rng_numpy.binomial(counts - counts_so_far, p)
             counts_so_far += read
 
+            # Add custom_reads if supplied and increment its frame counter
+            if custom_ramps is not None:
+                custom_so_far = custom_ramps[cr_count]
+                cr_count += 1
+            else:
+                custom_so_far = 0
+
             # Apply cosmic rays
             if crparam is not None:
                 old_instrumental_so_far = instrumental_so_far.copy()
@@ -411,9 +433,10 @@ def apportion_counts_to_resultants(
             if inv_linearity is not None:
                 # Apply inverse linearity
                 resultant_counts += inv_linearity.apply(
-                    counts_so_far + instrumental_so_far, electrons=True)
+                    counts_so_far + instrumental_so_far
+                    + custom_so_far, electrons=True)
             else:
-                resultant_counts += counts_so_far + instrumental_so_far
+                resultant_counts += counts_so_far + instrumental_so_far + custom_so_far
 
         # set the read count to the average of the resultant count
         resultants[i, ...] = resultant_counts / len(pi)
@@ -583,7 +606,8 @@ def make_l1(counts, read_pattern,
             rng=None, seed=None,
             gain=None, inv_linearity=None, crparam=None,
             persistence=None, tstart=None, saturation=None,
-            darkdecaysignal=None, ipc_model=None):
+            darkdecaysignal=None, custom_ramps=None, ipc_model=None):
+
     """Make an L1 image from a total electrons image.
 
     This apportions the total electrons among the different resultants and adds
@@ -623,13 +647,21 @@ def make_l1(counts, read_pattern,
         Dictionary with keys 'amplitude', 'time_constant', and 'sca'
         describing the dark decay signal.  If None, no dark decay is
         added.
+    custom_ramps : ndarray, an up-the-ramp sampled image cube that
+        will contribute additively to each read of pixel ramps
+        generated in apportion_counts_to_resultants().  Ideally has
+        dtype='i4'.  The array shape must be (nr, 4088, 4088) where
+        nr>=nreads, the total number of reads (not resultants) of the
+        ma table specified the metadata of the observation.  If
+        nr>nreads, the additional frames will be ignored.
 
     Returns
     -------
     l1 : np.ndarray[n_resultant, ny, nx]
         Resultants image array in DN including systematic effects
     dq : np.ndarray[n_resultant, ny, nx]
         DQ array marking saturated pixels and cosmic rays
+
     """
 
     tij = read_pattern_to_tij(read_pattern)
@@ -646,7 +678,7 @@ def make_l1(counts, read_pattern,
         pedestal=pedestal, pedestal_extra_noise=pedestal_extra_noise,
         inv_linearity=inv_linearity, crparam=crparam,
         persistence=persistence, tstart=tstart,
-        rng=rng, seed=seed)
+        rng=rng, seed=seed, custom_ramps=custom_ramps)
 
     # roman.addReciprocityFailure(resultants_object)