Ice/NVR and DetectorOps updates

 - Update DetectorOps pix_noise function to have more flexibility to manually insert custom noise properties
 - Add scale factors of nominal ice and NVR levels for the OTE and NIRCam, based on simulations from Chuck Bowers.

pynrc/nrc_utils.py

@@ -191,7 +191,8 @@ def bp_igood(bp, min_trans=0.001, fext=0.05):
 
 
 def read_filter(filter, pupil=None, mask=None, module=None, ND_acq=False,
-    ice_scale=None, nvr_scale=None, grism_order=1, **kwargs):
+    ice_scale=None, nvr_scale=None, ote_scale=None, nc_scale=None,
+    grism_order=1, **kwargs):
     """Read filter bandpass.
 
     Read in filter throughput curve from file generated by STScI.
@@ -213,12 +214,22 @@ def read_filter(filter, pupil=None, mask=None, module=None, ND_acq=False,
         ND acquisition square in coronagraphic mask.
     ice_scale : float
         Add in additional OTE H2O absorption. This is a scale factor
-        relative to 0.0131 um thickness.
+        relative to 0.0131 um thickness. Also includes about 0.0150 um of
+        photolyzed Carbon.
     nvr_scale : float
-        Add in additional NIRCam non-volatile residue. This is a scale
-        factor relative to 0.280 um thickness. If set to None, then 
-        assumes a scale factor of 1.0 as is contained in the NIRCam
-        filter curves. Setting nvr_scale=0 will remove these contributions.
+        Modify NIRCam non-volatile residue. This is a scale factor relative 
+        to 0.280 um thickness already built into filter throughput curves. 
+        If set to None, then assumes a scale factor of 1.0. 
+        Setting nvr_scale=0 will remove these contributions.
+    ote_scale : float
+        Scale factor of OTE contaminants relative to End of Life model. 
+        This is the same as setting ice_scale. Will override ice_scale value.
+    nc_scale : float
+        Scale factor for NIRCam contaminants relative to End of Life model.
+        This model assumes 0.189 um of NVR and 0.050 um of water ice on
+        the NIRCam optical elements. Setting this keyword will remove all
+        NVR contributions built into the NIRCam filter curves.
+        Overrides nvr_scale value.
     grism_order : int
         Option to use 2nd order grism throughputs instead. Useful if
         someone wanted to overlay the 2nd order contributions onto a 
@@ -444,7 +455,11 @@ def read_filter(filter, pupil=None, mask=None, module=None, ND_acq=False,
         w2 = wgood.max()
         wrange = w2 - w1
 
-
+    # OTE scaling (use ice_scale keyword)
+    if ote_scale is not None:
+        ice_scale = ote_scale
+    if nc_scale is not None:
+        nvr_scale = 0
     # Water ice and NVR additions (for LW channel only)
     if ((ice_scale is not None) or (nvr_scale is not None)) and ('LW' in channel):
         fname = conf.PYNRC_PATH + 'throughputs/ote_nc_sim_1.00.txt'
@@ -461,21 +476,51 @@ def read_filter(filter, pupil=None, mask=None, module=None, ND_acq=False,
             ttemp = np.insert(ttemp, 0, [1.0]) # Estimates for w<2.5um
             ttemp = np.append(ttemp, [1.0])    # Estimates for w>5.0um
             # Interpolate transmission onto filter wavelength grid
-            ttemp = np.interp(bp.wave/1e4, wtemp, ttemp, left=0, right=0)
-            th_ice = np.exp(ice_scale * np.log(ttemp))
-            th_new = th_ice * th_new
+            ttemp = np.interp(bp.wave/1e4, wtemp, ttemp)#, left=0, right=0)
+
+            # Scale is fraction of absorption feature depth, not of layer thickness
+            th_new = th_new * (1 - ice_scale * (1 - ttemp))
+            # th_ice = np.exp(ice_scale * np.log(ttemp))
+            # th_new = th_ice * th_new
 
         if nvr_scale is not None:
             ttemp = data['t_nvr']
             ttemp = np.insert(ttemp, 0, [1.0]) # Estimates for w<2.5um
             ttemp = np.append(ttemp, [1.0])    # Estimates for w>5.0um
             # Interpolate transmission onto filter wavelength grid
-            ttemp = np.interp(bp.wave/1e4, wtemp, ttemp, left=0, right=0)
+            ttemp = np.interp(bp.wave/1e4, wtemp, ttemp)#, left=0, right=0)
+
+            # Scale is fraction of absorption feature depth, not of layer thickness
+            # First, remove NVR contributions already included in throughput curve
+            th_new = th_new / ttemp
+            th_new = th_new * (1 - nvr_scale * (1 - ttemp))
 
             # The "-1" removes NVR contributions already included in
             # NIRCam throughput curves
-            th_nvr = np.exp((nvr_scale-1) * np.log(ttemp))
-            th_new = th_nvr * th_new
+            # th_nvr = np.exp((nvr_scale-1) * np.log(ttemp))
+            # th_new = th_nvr * th_new
+
+        if nc_scale is not None:
+            names = ['Wave', 'coeff'] # coeff is per um path length
+            path = conf.PYNRC_PATH
+            data_ice  = ascii.read(path + 'throughputs/h2o_abs.txt', names=names)
+            data_nvr  = ascii.read(path + 'throughputs/nvr_abs.txt', names=names)
+
+            w_ice = data_ice['Wave']
+            a_ice = data_ice['coeff']
+            a_ice = np.interp(bp.wave/1e4, w_ice, a_ice)
+
+            w_nvr = data_nvr['Wave']
+            a_nvr = data_nvr['coeff']
+            a_nvr = np.interp(bp.wave/1e4, w_nvr, a_nvr)
+
+            ttemp = np.exp(-0.189 * a_nvr - 0.050 * a_ice)
+            th_new = th_new * (1 - nc_scale * (1 - ttemp))
+
+            # ttemp = np.exp(-nc_scale*(a_nvr*0.189 + a_ice*0.05))
+            # th_new = ttemp * th_new
+
+
 
         # Create new bandpass
         bp = S.ArrayBandpass(bp.wave, th_new)
@@ -1080,7 +1125,7 @@ def psf_coeff(filter_or_bp, pupil=None, mask=None, module='A',
     # Get filter throughput and create bandpass
     if isinstance(filter_or_bp, six.string_types):
         filter = filter_or_bp
-        bp = read_filter(filter, pupil=pupil, mask=mask, module=module)
+        bp = read_filter(filter, pupil=pupil, mask=mask, module=module, **kwargs)
     else:
         bp = filter_or_bp
         filter = bp.name

pynrc/pynrc_core.py

@@ -184,7 +184,8 @@ def channel(self):
         """Detector channel 'SW' or 'LW' (inferred from detector ID)"""
         return 'LW' if self.detid.endswith('5') else 'SW'
 
-    def pixel_noise(self, fsrc=0.0, fzodi=0.0, fbg=0.0, ng=None, verbose=False, **kwargs):
+    def pixel_noise(self, fsrc=0.0, fzodi=0.0, fbg=0.0, rn=None, ktc=None, idark=None,
+        p_excess=None, ng=None, verbose=False, **kwargs):
         """Noise values per pixel.
         
         Return theoretical noise calculation for the specified MULTIACCUM exposure 
@@ -202,11 +203,22 @@ def pixel_noise(self, fsrc=0.0, fzodi=0.0, fbg=0.0, ng=None, verbose=False, **kw
             Flux of the zodiacal background in e-/sec/pix
         fbg : float or image
             Flux of telescope background in e-/sec/pix
+        idark : float or image
+            Option to specify dark current in e-/sec/pix.
+        rn : float
+            Option to specify Read Noise per pixel (e-).
+        ktc : float
+            Option to specify kTC noise (in e-). Only valid for single frame (n=1)
+        p_excess : array-like
+            Optional. An array or list of two elements that holds the parameters
+            describing the excess variance observed in effective noise plots.
+            By default these are both 0. For NIRCam detectors, recommended
+            values are [1.0,5.0] for SW and [1.5,10.0] for LW.
         ng : None or int or image
             Option to explicitly states number of groups. This is specifically
             used to enable the ability of only calculating pixel noise for
             unsaturated groups for each pixel. If a numpy array, then it should
-            be the same shape as `fsrc` image. By default will us `self.ngroup`.
+            be the same shape as `fsrc` image. By default will use `self.ngroup`.
         verbose : bool
             Print out results at the end.
 
@@ -227,11 +239,19 @@ def pixel_noise(self, fsrc=0.0, fzodi=0.0, fbg=0.0, ng=None, verbose=False, **kw
         ma = self.multiaccum
         if ng is None:
             ng = ma.ngroup
+        if rn is None:
+            rn = self.read_noise
+        if ktc is None:
+            ktc = self.ktc
+        if p_excess is None:
+            p_excess = self.p_excess
+        if idark is None:
+            idark = self.dark_current
 
         # Pixel noise per ramp (e-/sec/pix)
-        pn = pix_noise(ng, ma.nf, ma.nd2, tf=self.time_frame, \
-                       rn=self.read_noise, ktc=self.ktc, p_excess=self.p_excess, \
-                       idark=self.dark_current, fsrc=fsrc, fzodi=fzodi, fbg=fbg, **kwargs)
+        pn = pix_noise(ng, ma.nf, ma.nd2, tf=self.time_frame, 
+                       rn=rn, ktc=ktc, p_excess=p_excess, 
+                       idark=idark, fsrc=fsrc, fzodi=fzodi, fbg=fbg, **kwargs)
 
         # Divide by sqrt(Total Integrations)
         final = pn / np.sqrt(ma.nint)
@@ -436,6 +456,8 @@ def __init__(self, filter='F210M', pupil=None, mask=None, module='A', ND_acq=Fal
         # Specify ice and nvr scalings
         self._ice_scale = kwargs['ice_scale'] if 'ice_scale' in kwargs.keys() else None
         self._nvr_scale = kwargs['nvr_scale'] if 'nvr_scale' in kwargs.keys() else None
+        self._ote_scale = kwargs['ote_scale'] if 'ote_scale' in kwargs.keys() else None
+        self._nc_scale = kwargs['nc_scale'] if 'nc_scale' in kwargs.keys() else None
 
         # Let's figure out what keywords the user has set and try to 
         # interpret what he/she actually wants. If certain values have
@@ -623,7 +645,8 @@ def _update_bp(self):
         """Update bandpass based on filter, pupil, and module, etc."""
         self._bandpass = read_filter(self._filter, self._pupil, self._mask, 
                                      self.module, self.ND_acq,
-                                     ice_scale=self._ice_scale, nvr_scale=self._nvr_scale)
+                                     ice_scale=self._ice_scale, nvr_scale=self._nvr_scale,
+                                     ote_scale=self._ote_scale, nc_scale=self._nc_scale)
 
     def plot_bandpass(self, ax=None, color=None, title=None, **kwargs):
         """