Include twilight effects

doc/changes.rst

@@ -5,6 +5,7 @@ desisurvey change log
 0.10.1 (unreleased)
 -------------------
 
+* Add Ephemerides method to calculate sun altitude during twilight.
 * Set the ``EXTNAME`` keyword on the Table returned by ``Progress.get_exposures()``.
 
 0.10.0 (2017-11-09)

py/desisurvey/ephemerides.py

@@ -169,6 +169,12 @@ def __init__(self, start_date=None, stop_date=None, num_obj_steps=25,
         # tabulating the (ra,dec) of objects to avoid.
         t_obj = np.linspace(0., 1., num_obj_steps)
 
+        # Lookup sun altitude thresholds.
+        bright_max_alt = (
+            config.programs.BRIGHT.max_sun_altitude().to(u.rad).value)
+        dark_max_alt = (
+            config.programs.DARK.max_sun_altitude().to(u.rad).value)
+
         # Calculate ephmerides for each night.
         for day_offset in range(num_nights):
             day = self.start + day_offset * u.day
@@ -177,15 +183,13 @@ def __init__(self, start_date=None, stop_date=None, num_obj_steps=25,
             # Store local noon for this day.
             row['noon'] = day.mjd
             # Calculate bright twilight.
-            mayall.horizon = (
-                config.programs.BRIGHT.max_sun_altitude().to(u.rad).value)
+            mayall.horizon = bright_max_alt
             row['brightdusk'] = mayall.next_setting(
                 ephem.Sun(), use_center=True) + mjd0
             row['brightdawn'] = mayall.next_rising(
                 ephem.Sun(), use_center=True) + mjd0
             # Calculate dark / gray twilight.
-            mayall.horizon = (
-                config.programs.DARK.max_sun_altitude().to(u.rad).value)
+            mayall.horizon = dark_max_alt
             row['dusk'] = mayall.next_setting(
                 ephem.Sun(), use_center=True) + mjd0
             row['dawn'] = mayall.next_rising(
@@ -291,8 +295,7 @@ def get_moon_illuminated_fraction(self, mjd):
         Parameters
         ----------
         mjd : float or array
-            MJD values during a single night where the program should be
-            tabulated.
+            MJD values where the illuminated fraction should be tabulated.
 
         Returns
         -------
@@ -423,6 +426,60 @@ def is_full_moon(self, night, num_nights=None):
         return index - lo in sort_order[-num_nights:]
 
 
+def get_sun_altitude(row, mjd):
+    """Return the sun altitude at the specified time for twilight modeling.
+
+    This function only returns values in the range [-10, -20] deg
+    since values outside this range are either never used for observing
+    or else considered completely dark.
+
+    The calculation uses linear interpolation of the tabulated times when the
+    sun passes through the bright and dark maximum altitudes, nominally -13
+    and -15 deg, which provides an altitude accurate to better than 30"
+    between these limits.
+
+    Parameters
+    ----------
+    row : astropy.table.Row
+        A single row from the ephemerides astropy Table corresponding to the
+        night in question.
+    mjd : float
+        MJD value(s) during a single night where the sun altitude should be
+        calculated.
+
+    Returns
+    -------
+    astropy.units.Quantity
+        Sun altitude(s) in degrees at each input mjd value, clipped to
+        the range [-10, -20].
+    """
+    # Lookup the bright, dark threshold altitudes.
+    config = desisurvey.config.Configuration()
+    bright_max_alt = config.programs.BRIGHT.max_sun_altitude()
+    dark_max_alt = config.programs.DARK.max_sun_altitude()
+
+    # Validate input timestamp(s)
+    mjd = np.asarray(mjd)
+    scalar = mjd.ndim == 0
+    mjd = np.atleast_1d(mjd)
+    midnight = 0.5 * (row['dusk'] + row['dawn'])
+    if np.any(np.abs(mjd - midnight) > 0.5):
+        raise ValueError('mjd is not associated with the specified night')
+
+    # Use linear interpolation on the dusk data for times before midnight.
+    sun_alt = np.empty_like(mjd) * u.degree
+    at_dusk = mjd < midnight
+    sun_alt[at_dusk] = bright_max_alt + (dark_max_alt - bright_max_alt) * (
+        mjd[at_dusk] - row['brightdusk']) / (row['dusk'] - row['brightdusk'])
+
+    # Use linear interpolation on the dawn data for times after midnight.
+    at_dawn = ~at_dusk
+    sun_alt[at_dawn] = bright_max_alt + (bright_max_alt - dark_max_alt) * (
+        mjd[at_dawn] - row['brightdawn']) / (row['brightdawn'] - row['dawn'])
+
+    return np.clip(sun_alt, -20 * u.deg, -10 * u.deg)
+
+
 def get_object_interpolator(row, object_name, altaz=False):
     """Build an interpolator for object location during one night.
 

py/desisurvey/etc.py

@@ -121,11 +121,67 @@ def airmass_exposure_factor(airmass):
     return np.power((X / X0), 1.25)
 
 
+# Linear regression coefficients for converting scattered twilight r-band
+# magnitude into an exposure-time correction factor.
+_twilightCoefficients = np.array([
+    -8.70618444805, 64323409.928312987, 739.2405272968557,
+    -18763.675456604633, 158813.46020567254])
+
+def twilight_exposure_factor(sun_alt, sun_daz, airmass):
+    """Calculate exposure time factor due to scattered twilight.
+
+    The returned factor is relative to dark conditions when the moon is
+    below the local horizon.
+
+    This factor is based on a study of SNR for BGS threshold targets.
+    See ``surveysim:doc/nb/BGS_ETC_Study.ipynb`` for details.
+
+    Parameters
+    ----------
+    sun_alt : float
+        Altitude angle of the sun above the horizon in degrees. Must be in the
+        range [-90, -12] deg. Values below -18 deg are considered dark, with
+        no scattered sun.
+    sun_daz : float
+        Relative azimuth angle in degrees between the pointing and the sun.
+    airmass : float
+        Airmass used for observing this tile, must be >= 1.
+
+    Returns
+    -------
+    float
+        Dimensionless factor that exposure time should be increased to
+        account for increased sky brightness due to scattered twilight.
+        Will be 1 when the sun is below -18 degrees.
+    """
+    if (sun_alt < -90) or (sun_alt > -12):
+        raise ValueError('Got invalid sun_alt outside [-90,-12].')
+    if airmass < 1:
+        raise ValueError('Got invalid airmass < 1.')
+
+    # No exposure penalty when the sun is far enough below the horizon.
+    if sun_alt < -18:
+        return 1.
+
+    # Estimate the altitude angle corresponding to this observing airmass.
+    # We invert eqn.3 of KS1991 for this (instead of eqn.14).
+    obs_zenith = np.arcsin(np.sqrt((1 - airmass ** -2) / 0.96)) * u.rad
+    obj_altitude = 90 * u.deg - obs_zenith.to(u.deg)
+
+    # Calculate scattered twilight r-band brightness.
+    r = specsim.atmosphere.twilight_surface_brightness(
+        obj_altitude, sun_alt * u.deg, sun_daz * u.deg).value
+
+    # Evaluate the linear regression model.
+    X = np.array((1, np.exp(-r), 1/r, 1/r**2, 1/r**3))
+    return _twilightCoefficients.dot(X)
+
+
 # Linear regression coefficients for converting scattered moon V-band
 # magnitude into an exposure-time correction factor.
 _moonCoefficients = np.array([
-    -8.83964463188, -7372368.5041596508, 775.17763895781638,
-    -20185.959363990656, 174143.69095766739])
+    -9.6827757062, -5454732.2391116023, 805.44541757440129,
+    -20274.598621708101, 170436.98267100245])
 
 # V-band extinction coefficient to use in the scattered moonlight model.
 # See specsim.atmosphere.krisciunas_schaefer for details.
@@ -140,14 +196,8 @@ def moon_exposure_factor(moon_frac, moon_sep, moon_alt, airmass):
     This factor is based on a study of SNR for ELG targets and designed to
     achieve a median SNR of 7 for a typical ELG [OII] doublet at the lower
     flux limit of 8e-17 erg/(cm2 s A), averaged over the expected ELG target
-    redshift distribution 0.6 < z < 1.7.
-
-    TODO:
-    - Check the assumption that exposure time scales with SNR ** -0.5.
-    - Check if this ELG-based analysis is also valid for BGS targets.
-
-    For details, see the jupyter notebook doc/nb/ScatteredMoon.ipynb in
-    this package.
+    redshift distribution 0.6 < z < 1.7. See
+    ``surveysim:doc/nb/ScatteredMoon.ipynb`` for details.
 
     Parameters
     ----------
@@ -191,14 +241,16 @@ def moon_exposure_factor(moon_frac, moon_sep, moon_alt, airmass):
     # We invert eqn.3 of KS1991 for this (instead of eqn.14).
     obs_zenith = np.arcsin(np.sqrt((1 - airmass ** -2) / 0.96)) * u.rad
 
-    # Calculate scattered moon V-band brightness at each pixel.
+    # Calculate scattered moon V-band brightness.
     V = specsim.atmosphere.krisciunas_schaefer(
         obs_zenith, moon_zenith, separation_angle,
         moon_phase, _vband_extinction).value
 
     # Evaluate the linear regression model.
-    X = np.array((1, np.exp(-V), 1/V, 1/V**2, 1/V**3))
-    return _moonCoefficients.dot(X)
+    X = np.array((1., np.exp(-V), 1/V, 1/V**2, 1/V**3))
+
+    # Clip result to 1 since fits goes just below 1 at the faint end (V>26).
+    return np.maximum(1., _moonCoefficients.dot(X))
 
 
 def exposure_time(program, seeing, transparency, airmass, EBV,

py/desisurvey/progress.py

@@ -115,6 +115,20 @@ def __init__(self, restore=None, max_exposures=32):
             table['moonsep'] = astropy.table.Column(
                 length=num_tiles, shape=(max_exposures,), format='%.1f',
                 description='Moon-tile separation angle in degrees', unit='deg')
+            table['moonv'] = astropy.table.Column(
+                length=num_tiles, shape=(max_exposures,), format='%.1f',
+                description='V-band magnitude of scattered moonlight')
+            table['sunalt'] = astropy.table.Column(
+                length=num_tiles, shape=(max_exposures,), format='%.1f',
+                description='Sun altitude angle in degrees, clipped below -20',
+                unit='deg')
+            table['sundaz'] = astropy.table.Column(
+                length=num_tiles, shape=(max_exposures,), format='%.1f',
+                description='Sun-tile relative azimuth angle in degrees',
+                unit='deg')
+            table['sunr'] = astropy.table.Column(
+                length=num_tiles, shape=(max_exposures,), format='%.1f',
+                description='r-band magnitude of twilight scattered sun')
             # Copy tile data.
             table['tileid'] = tiles['TILEID']
             table['pass'] = tiles['PASS']
@@ -415,7 +429,8 @@ def copy_range(self, mjd_min=None, mjd_max=None):
         return Progress(restore=table)
 
     def add_exposure(self, tile_id, start, exptime, snr2frac, airmass, seeing,
-                     transparency, moonfrac, moonalt, moonsep):
+                     transparency, moonfrac, moonalt, moonsep, moonv,
+                     sunalt, sundaz, sunr):
         """Add a single exposure to the progress.
 
         Parameters
@@ -440,6 +455,16 @@ def add_exposure(self, tile_id, start, exptime, snr2frac, airmass, seeing,
             Moon altitude angle in degrees.
         moonsep : float
             Moon-tile separation angle in degrees.
+        moonv : float
+            The V-band magnitude of scattered moonlight, or -np.inf if
+            the moon is below the horizon.
+        sunalt : float
+            Sun altitude angle in degrees, clipped below -20.
+        sundaz : float
+            Sun-tile relative azimuth angle in degrees.
+        sunr : float
+            The r-band magnitude of twilight scattered sun, or -np.inf
+            if the sun is below -18 deg altitude.
         """
         mjd = start.mjd
         self.log.info(
@@ -478,6 +503,10 @@ def add_exposure(self, tile_id, start, exptime, snr2frac, airmass, seeing,
         row['moonfrac'][num_exp] = moonfrac
         row['moonalt'][num_exp] = moonalt
         row['moonsep'][num_exp] = moonsep
+        row['moonv'][num_exp] = moonv
+        row['sunalt'][num_exp] = sunalt
+        row['sundaz'][num_exp] = sundaz
+        row['sunr'][num_exp] = sunr
 
         # Update this tile's status.
         row['status'] = 1 if row['snr2frac'].sum() < self.min_snr2 else 2
@@ -486,7 +515,7 @@ def get_exposures(self, start=None, stop=None,
                       tile_fields='tileid,pass,ra,dec,ebmv',
                       exp_fields=('night,mjd,exptime,seeing,transparency,' +
                                   'airmass,moonfrac,moonalt,moonsep,' +
-                                  'program,flavor')):
+                                  'moonv,sunalt,sundaz,sunr,program,flavor')):
         """Create a table listing exposures in time order.
 
         Parameters
@@ -515,7 +544,8 @@ def get_exposures(self, start=None, stop=None,
         Returns
         -------
         astropy.table.Table
-            Table with the specified columns as uppercase and one row per exposure.
+            Table with the specified columns as uppercase and one row
+            per exposure.
         """
         # Get MJD range to show.
         if start is None:
@@ -552,6 +582,8 @@ def get_exposures(self, start=None, stop=None,
         output = astropy.table.Table()
         output.meta['EXTNAME'] = 'EXPOSURES'
         for name in tile_fields.split(','):
+            if name == '':
+                continue # handle tile_fields=''
             name = name.lower()
             if name == 'index':
                 output[name.upper()] = tile_index
@@ -569,6 +601,8 @@ def get_exposures(self, start=None, stop=None,
                         'Invalid tile field name: {0}.'.format(name))
                 output[name.upper()] = table[name][tile_index]
         for name in exp_fields.split(','):
+            if name == '':
+                continue # handle exp_fields=''
             name = name.lower()
             if name == 'snr2cum':
                 snr2cum = np.cumsum(
@@ -580,7 +614,8 @@ def get_exposures(self, start=None, stop=None,
                 mjd = table['mjd'].flatten()[order]
                 night = np.empty(len(mjd), dtype=(str, 8))
                 for i in range(len(mjd)):
-                    night[i] = str(desisurvey.utils.get_date(mjd[i])).replace('-', '')
+                    night[i] = str(desisurvey.utils.get_date(mjd[i])
+                                   ).replace('-', '')
                 output[name.upper()] = astropy.table.Column(
                     night,
                     description='Date at start of night when exposure taken')
@@ -605,8 +640,8 @@ def get_exposures(self, start=None, stop=None,
                     program[exppass < 4] = 'DARK'
                     program[exppass == 4] = 'GRAY'
 
-                output[name.upper()] = astropy.table.Column(program,
-                                                            description='Program name')
+                output[name.upper()] = astropy.table.Column(
+                    program, description='Program name')
             elif name == 'flavor':
                 flavor = np.empty(len(exppass), dtype=(str, 7))
                 flavor[:] = 'science'