revised exposure factor implementation for bright time

py/desisurvey/ephem.py

@@ -195,7 +195,8 @@ def __init__(self, start_date, stop_date, num_obj_steps=25, restore=None):
         # Add (ra,dec) arrays for each object that we need to avoid and
         # check that ephem has a model for it.
         models = {}
-        for name in config.avoid_bodies.keys:
+        print(config.avoid_bodies.keys)
+        for name in list(config.avoid_bodies.keys)+['sun']:
             models[name] = getattr(ephem, name.capitalize())()
             self._table[name + '_ra'] = astropy.table.Column(
                 length=num_nights, shape=(num_obj_steps,), format='%.2f',

py/desisurvey/etc.py

@@ -11,6 +11,7 @@
 from __future__ import print_function, division
 
 import numpy as np
+from itertools import chain, combinations_with_replacement
 
 import astropy.units as u
 
@@ -20,7 +21,7 @@
 
 import desisurvey.config
 import desisurvey.tiles
-
+        
 
 def seeing_exposure_factor(seeing):
     """Scaling of exposure time with seeing, relative to nominal seeing.
@@ -138,6 +139,7 @@ def airmass_exposure_factor(airmass):
 # See specsim.atmosphere.krisciunas_schaefer for details.
 _vband_extinction = 0.15154
 
+
 def moon_exposure_factor(moon_frac, moon_sep, moon_alt, airmass):
     """Calculate exposure time factor due to scattered moonlight.
 
@@ -208,6 +210,146 @@ def moon_exposure_factor(moon_frac, moon_sep, moon_alt, airmass):
     return _moonCoefficients.dot(X)
 
 
+def bright_exposure_factor(airmass, moon_frac, moon_sep, moon_alt, sun_sep, sun_alt):
+    """ Calculate exposure time factor based on airmass, moon, and sun parameters. 
+
+    :param moon_frac: 
+        Illuminated fraction of the moon within range [0,1].
+
+    :param moon_alt:
+        Altitude angle of the moon above the horizon in degrees within range [-90,90].
+
+    :param moon_sep:
+        Array of separation angles between field center and moon in degrees within the
+        range [0,180].
+
+    :param sun_alt:
+        Altitude angle of the sunin degrees
+
+    :param sun_sep: 
+        Arry of separations angles between field center and sun in degrees
+
+    :param airmass:
+        Array of airmass used for observing this tile, must be >= 1.
+
+    :returns expfactors: 
+        Dimensionless factors that exposure time should be increased to
+        account for increased sky brightness due to scattered moonlight.
+        Will be 1 when the moon is below the horizon.
+    """
+    # exposure_factor = 1 when moon is below the horizon and sun is below -20.
+    if moon_alt < 0 and sun_alt < -18.:
+        return np.ones(len(airmass))  
+
+    # check inputs 
+    moon_sep    = moon_sep.flatten() 
+    sun_sep     = sun_sep.flatten()
+    airmass     = airmass.flatten() 
+    if (moon_frac < 0) or (moon_frac > 1):
+        raise ValueError('Got invalid moon_frac outside [0,1].')
+    if (moon_alt < -90) or (moon_alt > 90):
+        raise ValueError('Got invalid moon_alt outside [-90,+90].')
+    if (moon_sep.min() < 0) or (moon_sep.max() > 180):
+        raise ValueError('Got invalid moon_sep outside [0,180].')
+    if airmass.min() < 1:
+        raise ValueError('Got invalid airmass < 1.')
+
+    # check size of inputs  
+    nexp = len(airmass) 
+    assert len(moon_sep) == nexp
+    assert len(sun_sep) == nexp
+
+    # calculate exposure factor 
+    config = desisurvey.config.Configuration()
+
+    # sky brightness at 5000A for reference BGS exposure  
+    Isky5000_ref = 3.6956611461286966 # 1e-17 erg/s/cm^2/A/arcsec^2
+
+    # exposure time for reference BGS exposure 
+    tref = getattr(config.nominal_exposure_time, 'BRIGHT')().to(u.s).value
+
+    # sky brightness at 5000A for observing conditions 
+    Isky5000_exp = bright_Isky5000_notwilight_regression(
+            airmass,
+            np.repeat(moon_frac, nexp), 
+            moon_sep, 
+            np.repeat(moon_alt, nexp))
+
+     # add twilight contribution 
+    if sun_alt >= -18.:
+        Isky5000_exp += np.clip(bright_Isky5000_twilight_regression(
+                airmass, 
+                sun_sep, 
+                np.repeat(sun_alt, nexp)), 0, None) 
+
+    # calculate exposure factor 
+    fibflux5000_ref = Isky5000_ref * 1.862089 # 1e-17 erg/s/cm^2/A
+    fibflux5000_exp = Isky5000_exp * 1.862089 # 1e-17 erg/s/cm^2/A
+
+    # 0.0629735016982807 = 1e-17 x (photons per bin) x throughput) at 5000A 
+    sky_photon_per_sec_ref = fibflux5000_ref * 0.0629735016982807
+    sky_photon_per_sec_exp = fibflux5000_exp * 0.0629735016982807
+
+    # (read noise)^2 at 5000A 
+    RNsq5000 = 56.329457658891435 # photon^2 
+
+    # solve the following: 
+    # S x tref / sqrt(sky_ref * tref + RN^2) = S x texp / sqrt(sky_exp * texp + RN^2)
+    texp = 0.5 * (
+            (tref * np.sqrt(4 * sky_photon_per_sec_ref * RNsq5000 * tref + sky_photon_per_sec_exp**2 * tref**2 + 4 * RNsq5000**2))/(sky_photon_per_sec_ref * tref + RNsq5000) +
+            (sky_photon_per_sec_exp * tref**2)/(sky_photon_per_sec_ref * tref + RNsq5000))
+    return texp/tref 
+
+
+def bright_Isky5000_notwilight_regression(airmass, moon_frac, moon_sep, moon_alt): 
+    ''' polynomial regression model for bright sky surface brightness at 5000A
+    *without twilight*. The regression model was fit using observed sky surface
+    brightnesses from 
+    * DESI SV1 bright exposures
+    * DESI CMX exposures with transparency > 0.9 
+    * BOSS exposures with sun alt < -18 
+
+    see
+    https://github.com/desi-bgs/bgs-cmxsv/blob/4c5f124164b649c595cd2dca87d14ba9f3b2c64d/doc/nb/sv1_sky_model_fit.ipynb
+    for detials.
+    '''
+    # polynomial regression cofficients for estimating exposure time factor during
+    # non-twilight from airmass, moon_frac, moon_sep, moon_alt  
+    coeffs = np.array([
+        1.22875526e+00,  1.91591548e-01,  3.17313759e+00,  5.22047416e-02,
+       -3.87652985e-02, -5.33224507e-01,  4.63261325e+00,  1.13410640e-02,
+       -1.01921126e-02,  1.06314395e+00,  7.26049602e-02, -1.08328690e-01,
+       -8.95312945e-04, -5.59394346e-04,  7.99072084e-04])
+
+    theta = np.atleast_2d(np.array([airmass, moon_frac, moon_alt, moon_sep]).T)
+
+    combs = chain.from_iterable(combinations_with_replacement(range(4), i) for i in range(0, 3))
+
+    theta_transform = np.empty((theta.shape[0], len(coeffs)))
+    for i, comb in enumerate(combs):
+        theta_transform[:, i] = theta[:, comb].prod(1)
+
+    return np.dot(theta_transform, coeffs.T)
+
+
+def bright_Isky5000_twilight_regression(airmass, sun_sep, sun_alt): 
+    ''' polynomial regression model for twilight contribution to the sky
+    surface brightness at 5000A. The regression model was fit using 
+    * BOSS exposures with sun alt > -18.
+    '''
+    norder = 1
+    skymodel_coeff = np.array([2.00474700e+00, 3.19827604e+00, 3.13212960e-01, 2.69673079e-03])
+
+    theta = np.atleast_2d(np.array([airmass, sun_alt, sun_sep]).T)
+
+    combs = chain.from_iterable(combinations_with_replacement(range(3), i) for i in range(0, norder+1))
+    theta_transform = np.empty((theta.shape[0], len(skymodel_coeff)))
+    for i, comb in enumerate(combs):
+        theta_transform[:, i] = theta[:, comb].prod(1)
+
+    return np.dot(theta_transform, skymodel_coeff.T)
+
+
 def exposure_time(program, seeing, transparency, airmass, EBV,
                   moon_frac, moon_sep, moon_alt):
     """Calculate the total exposure time for specified observing conditions.
@@ -304,8 +446,8 @@ def __init__(self, save_history=False):
         for program in desisurvey.tiles.Tiles.PROGRAMS:
             self.TEXP_TOTAL[program] = getattr(config.nominal_exposure_time, program)().to(u.day).value
         # Temporary hardcoded exposure factors for moon-up observing.
-        self.TEXP_TOTAL['GRAY'] *= 1.1
-        self.TEXP_TOTAL['BRIGHT'] *= 1.33
+        #self.TEXP_TOTAL['GRAY'] *= 1.1
+        #self.TEXP_TOTAL['BRIGHT'] *= 1.33
 
         # Initialize model of exposure time dependence on seeing.
         self.seeing_coefs = np.array([12.95475751, -7.10892892, 1.21068726])

py/desisurvey/plan.py

@@ -18,7 +18,7 @@
 import desisurvey.utils
 
 
-def load_design_hourangle(name='surveyinit.fits'):
+def load_design_hourangle(name='surveyinit.fits', bgs_footprint=None):
     """Load design hour-angle assignments from disk.
 
     Reads column 'HA' from binary table HDU 'DESIGN'. This is the format
@@ -46,6 +46,8 @@ def load_design_hourangle(name='surveyinit.fits'):
     else:
         with astropy.io.fits.open(fullname, memmap=False) as hdus:
             HA = hdus['DESIGN'].data['HA'].copy()
+        if bgs_footprint is not None: 
+            HA = HA[tiles._reduced_footprint]
         if HA.shape != (tiles.ntiles,):
             raise ValueError('Read unexpected HA shape.')
     return HA

py/desisurvey/scheduler.py

@@ -51,7 +51,7 @@ class Scheduler(object):
         1D array of design hour angles to use in degrees, or use
         :func:`desisurvey.plan.load_design_hourangle` when None.
     """
-    def __init__(self, restore=None, design_hourangle=None):
+    def __init__(self, restore=None, design_hourangle=None, bgs_footprint=None):
         self.log = desiutil.log.get_logger()
         # Load our configuration.
         config = desisurvey.config.Configuration()
@@ -63,11 +63,11 @@ def __init__(self, restore=None, design_hourangle=None):
         self.max_airmass = desisurvey.utils.cos_zenith_to_airmass(np.sin(config.min_altitude()))
         self.max_ha = config.max_hour_angle().to(u.deg).value
         # Load static tile info.
-        self.tiles = desisurvey.tiles.get_tiles()
+        self.tiles = desisurvey.tiles.get_tiles(bgs_footprint=bgs_footprint)
         ntiles = self.tiles.ntiles
         # Check hourangles.
         if design_hourangle is None:
-            self.design_hourangle = desisurvey.plan.load_design_hourangle()
+            self.design_hourangle = desisurvey.plan.load_design_hourangle(bgs_footprint=bgs_footprint)
         else:
             self.design_hourangle = np.asarray(design_hourangle)
         if self.design_hourangle.shape != (self.tiles.ntiles,):
@@ -264,10 +264,13 @@ def init_night(self, night, use_twilight=False):
         self.in_night_pool[avoid_idx] = False
         # Initialize moon tracking during this night.
         self.moon_DECRA = desisurvey.ephem.get_object_interpolator(self.night_ephem, 'moon', altaz=False)
-        #self.moon_ALTAZ = desisurvey.ephem.get_object_interpolator(self.night_ephem, 'moon', altaz=True)
+        self.moon_ALTAZ = desisurvey.ephem.get_object_interpolator(self.night_ephem, 'moon', altaz=True)
+        # Initialize sun tracking during this night.
+        self.sun_DECRA = desisurvey.ephem.get_object_interpolator(self.night_ephem, 'sun', altaz=False) 
+        self.sun_ALTAZ = desisurvey.ephem.get_object_interpolator(self.night_ephem, 'sun', altaz=True) 
 
-    def next_tile(self, mjd_now, ETC, seeing, transp, skylevel, HA_sigma=15., greediness=0.,
-                  program=None):
+    def next_tile(self, mjd_now, ETC, seeing, transp, skylevel, HA_sigma=15.,
+            greediness=0., use_brightsky=False, program=None):
         """Select the next tile to observe.
 
         The :meth:`init_night` method must be called before calling this
@@ -307,6 +310,10 @@ def next_tile(self, mjd_now, ETC, seeing, transp, skylevel, HA_sigma=15., greedi
             values will depend on the value of ``HA_sigma`` and how exposure
             factors are calculated. Refer to the equation above for details.
             Must be between 0 and 1.
+        use_brightsky : bool 
+            If True use improved bright sky model in next_tile selection to
+            calculate exposure factor for bright sky. If False, exposure factor
+            does not include bright sky model.
         program : string
             PROGRAM of tile to select.  Default of None selects the appropriate
             PROGRAM given current moon/twilight conditions.  Forcing a particular
@@ -369,10 +376,11 @@ def next_tile(self, mjd_now, ETC, seeing, transp, skylevel, HA_sigma=15., greedi
         if not np.any(self.tile_sel):
             # No tiles left to observe after airmass cut.
             return None, None, None, None, None, program, mjd_program_end
+
         # Is the moon up?
         if mjd_now > self.night_ephem['moonrise'] and mjd_now < self.night_ephem['moonset']:
             moon_is_up = True
-            # Calculate the moon (RA,DEC).
+            # calculate the moon (RA,DEC).
             moonDEC, moonRA = self.moon_DECRA(mjd_now)
             # Identify tiles that are too close to the moon to observe now.
             too_close = desisurvey.utils.separation_matrix(
@@ -386,12 +394,19 @@ def next_tile(self, mjd_now, ETC, seeing, transp, skylevel, HA_sigma=15., greedi
                 return None, None, None, None, None, program, mjd_program_end
         else:
             moon_is_up = False
+
         # Estimate exposure factors for all available tiles.
         self.exposure_factor[:] = 1e8
         self.exposure_factor[self.tile_sel] = self.tiles.dust_factor[self.tile_sel]
-        self.exposure_factor[self.tile_sel] *= desisurvey.etc.airmass_exposure_factor(self.airmass[self.tile_sel])
+        if use_brightsky and program == 'BRIGHT': 
+            self.exposure_factor[self.tile_sel] *= \
+                    self.update_exposure_factor(mjd_now, self.tiles.tileID[self.tile_sel])
+        else: 
+            self.exposure_factor[self.tile_sel] *= \
+                    desisurvey.etc.airmass_exposure_factor(self.airmass[self.tile_sel])
         # Apply global weather factors that are the same for all tiles.
         self.exposure_factor[self.tile_sel] /= ETC.weather_factor(seeing, transp)
+
         if not np.any(self.tile_sel):
             return None, None, None, None, None, program, mjd_program_end
         # Calculate (the log of a) Gaussian multiplicative penalty for
@@ -442,3 +457,67 @@ def survey_completed(self):
         """Test if all tiles have been completed.
         """
         return self.completed_by_pass.sum() == self.tiles.ntiles
+
+    def update_exposure_factor(self, mjd, tileid, return_obs_cond=False): 
+        """ get updated exposure factor on this night given mjd, and tile ID.
+        """
+        # get tile index 
+        idx = [] 
+        for _id in np.atleast_1d(tileid): 
+            idx.append(np.where(self.tiles.tileID == _id)[0])
+        idx = np.array(idx).flatten() 
+        assert len(idx) > 0  
+
+        # (RA,DEC) of the moon and sun at mjd
+        moonDEC, moonRA = self.moon_DECRA(mjd)
+        moonALT, moonAZ = self.moon_ALTAZ(mjd) 
+        sunDEC, sunRA = self.sun_DECRA(mjd)
+        sunALT, sunAZ = self.sun_ALTAZ(mjd) 
+
+        # moon illumination 
+        moonILL = self.night_ephem['moon_illum_frac']
+
+        # calculate moon and sun separation 
+        moonSEP = desisurvey.utils.separation_matrix(
+            [moonRA], [moonDEC],
+            self.tiles.tileRA[idx], self.tiles.tileDEC[idx])
+        sunSEP = desisurvey.utils.separation_matrix(
+            [sunRA], [sunDEC],
+            self.tiles.tileRA[idx], self.tiles.tileDEC[idx])
+
+        fexp = desisurvey.etc.bright_exposure_factor(
+                self.airmass[idx], moonILL, moonSEP, moonALT, sunSEP, sunALT)
+        if not return_obs_cond: 
+            return fexp
+        else: 
+            return fexp, moonILL, moonSEP, moonALT, sunSEP, sunALT
+
+    def get_observing_conditions(self, mjd, tileid): 
+        """ get observing conditions this night given mjd, and tile ID.
+        """
+        # get tile index 
+        idx = [] 
+        for _id in np.atleast_1d(tileid): 
+            idx.append(np.where(self.tiles.tileID == _id)[0])
+        idx = np.array(idx).flatten() 
+        assert len(idx) > 0  
+
+        # (RA,DEC) of the moon and sun at mjd
+        moonDEC, moonRA = self.moon_DECRA(mjd)
+        moonALT, moonAZ = self.moon_ALTAZ(mjd) 
+        sunDEC, sunRA = self.sun_DECRA(mjd)
+        sunALT, sunAZ = self.sun_ALTAZ(mjd) 
+
+        # moon illumination 
+        moonILL = self.night_ephem['moon_illum_frac']
+
+        # calculate moon and sun separation 
+        moonSEP = desisurvey.utils.separation_matrix(
+            [moonRA], [moonDEC],
+            self.tiles.tileRA[idx], self.tiles.tileDEC[idx])
+        sunSEP = desisurvey.utils.separation_matrix(
+            [sunRA], [sunDEC],
+            self.tiles.tileRA[idx], self.tiles.tileDEC[idx])
+
+        return moonILL, moonSEP, moonALT, sunSEP, sunALT
+