diff --git a/romanisim/image.py b/romanisim/image.py index 466555dc..cf59f892 100644 --- a/romanisim/image.py +++ b/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: diff --git a/romanisim/l1.py b/romanisim/l1.py index b084f858..e92a1f6f 100644 --- a/romanisim/l1.py +++ b/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,6 +647,13 @@ 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 ------- @@ -630,6 +661,7 @@ def make_l1(counts, read_pattern, 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)