Fix readnoise pixfrac scaling (#350)

* Make astropy times use isot under-the-hood.

* Fix read noise scaling with pixscalefrac in L3 generation.

Author: schlafly

romanisim/l3.py

@@ -398,9 +398,10 @@ def simulate(shape, wcs, efftimes, filter_name, catalog, nexposures=1,
     # Set effective read noise
     if effreadnoise is None:
         readnoise = np.median(romanisim.parameters.reference_data['readnoise'])
+        # read noise in DN
         gain = np.median(romanisim.parameters.reference_data['gain'])
         effreadnoise = (
-            np.sqrt(2) * readnoise * gain)
+            np.sqrt(2) * readnoise * gain)  # electron, difference of two reads
         # sqrt(2) from subtracting one read from another
         effreadnoise /= (np.median(efftimes * pixscalefrac ** 2) / nexposures)
         # divided by the typical exposure length
@@ -412,10 +413,10 @@ def simulate(shape, wcs, efftimes, filter_name, catalog, nexposures=1,
         # note that we are ignoring all of the individual reads, which also
         # contribute to reducing the effective read noise.  Pass --effreadnoise
         # if you want to do better than this!
-        effreadnoise = effreadnoise * etomjysr  # electron -> MJy/sr
+        effreadnoise = effreadnoise * etomjysr * pixscalefrac ** 2  # electron -> MJy/sr
+        # factor of pixfrac ** 2 because etomjysr is per L3 pixel, but
+        # this effreadnoise calculation is per native pixel
         # converting to MJy/sr units
-    else:
-        effreadnoise = 0
 
     chromatic = False
     if (len(catalog) > 0 and not isinstance(catalog, astropy.table.Table)
@@ -481,7 +482,7 @@ def simulate_cps(image, filter_name, efftimes, objlist=None, psf=None,
                  xpos=None, ypos=None, coord=None, sky=0, bandpass=None,
                  effreadnoise=None, maggytoes=None, etomjysr=None,
                  rng=None, seed=None, ignore_distant_sources=10,):
-    """Simulate average MegaJanskies per steradian in a single SCA.
+    """Simulate average MJy/sr in a single SCA.
 
     Parameters
     ----------
@@ -507,10 +508,10 @@ def simulate_cps(image, filter_name, efftimes, objlist=None, psf=None,
     effreadnoise : float
         Effective read noise for mosaic (MJy / sr)
     maggytoes: float
-        Factor to convert electrons to MJy / sr; one maggy makes
+        Factor to convert e/s to MJy / sr; one maggy makes
         this many e/s.
     etomjysr : float
-        Factor to convert electron to MJy/sr;  one e/s/pix corresponds
+        Factor to convert e/s to MJy / sr; one e/s/coadd pix corresponds
         to this MJy/sr.
     rng : galsim.BaseDeviate
         random number generator

romanisim/tests/test_l3.py

@@ -532,7 +532,7 @@ def test_scaling():
     twcs1 = romanisim.wcs.create_tangent_plane_gwcs(
         (npix / 2, npix / 2), pscale, coord)
     twcs2 = romanisim.wcs.create_tangent_plane_gwcs(
-        (npix / 2, npix / 2), pscale / 2, coord)
+        (npix, npix), pscale / 2, coord)
 
     im1, extras1 = l3.simulate(
         (npix, npix), twcs1, exptime, imdict['filter_name'],
@@ -600,3 +600,24 @@ def test_scaling():
 
     # these all match to a few percent; worst case in initial test run
     # was err3fracdiff of 0.039.
+
+    # test read noise scaling
+    im1, extras1 = l3.simulate(
+        (npix, npix), twcs1, exptime, imdict['filter_name'],
+        imdict['tabcatalog'], seed=rng_seed, effreadnoise=None, sky=0)
+
+    # half pixel scale
+    im2, extras2 = l3.simulate(
+        (npix * 2, npix * 2), twcs2, exptime, imdict['filter_name'],
+        imdict['tabcatalog'], seed=rng_seed, effreadnoise=None, sky=0)
+
+    for im in [im1, im2]:
+        assert np.abs(mad_std(im1.data) / np.median(im1.err) - 1) < 0.1
+        assert np.abs(np.median(im1.err) /
+                      np.median(np.sqrt(im1.var_rnoise)) - 1) < 0.1
+    assert np.abs(np.median(im2.err) / np.median(im1.err) / 4 - 1) < 0.1
+    # 2x finer sampling -> 4x higher read noise in per-unit-time units
+    # because constant read noise is divided by 4x smaller exposure time
+    # the 2x finer sampling effectively means that you need to take the
+    # four exposures and interleave them, each using 1/4 of the exposure
+    # time