diff --git a/py/desisurvey/ephem.py b/py/desisurvey/ephem.py index caf47b45..3c4fd3b5 100644 --- a/py/desisurvey/ephem.py +++ b/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', diff --git a/py/desisurvey/etc.py b/py/desisurvey/etc.py index 5bc8bee4..dd5d18bc 100644 --- a/py/desisurvey/etc.py +++ b/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]) diff --git a/py/desisurvey/plan.py b/py/desisurvey/plan.py index 1e4915bb..457ade75 100644 --- a/py/desisurvey/plan.py +++ b/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 diff --git a/py/desisurvey/scheduler.py b/py/desisurvey/scheduler.py index d87e3871..96b52f77 100644 --- a/py/desisurvey/scheduler.py +++ b/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 + diff --git a/py/desisurvey/tiles.py b/py/desisurvey/tiles.py index 3c160847..ce7b8b7b 100644 --- a/py/desisurvey/tiles.py +++ b/py/desisurvey/tiles.py @@ -45,7 +45,7 @@ class Tiles(object): Name of the tiles file to use or None for the default specified in our configuration. """ - def __init__(self, tiles_file=None): + def __init__(self, tiles_file=None, bgs_footprint=None): log = desiutil.log.get_logger() config = desisurvey.config.Configuration() # Read the specified tiles file. @@ -62,6 +62,13 @@ def __init__(self, tiles_file=None): unknown = set(tile_programs) - set(self.PROGRAMS) if unknown and not commissioning: raise RuntimeError('Cannot schedule unknown program(s): {}.'.format(unknown)) + + if bgs_footprint is not None: + # impose reduced footprint for BRIGHT tiles + tiles, _reduced_footprint = self._reduced_bgs_footprint(tiles, bgs_footprint) + self._reduced_footprint = _reduced_footprint.copy() + print(tiles.shape) + # Copy tile arrays. self.tileID = tiles['TILEID'].copy() self.passnum = tiles['PASS'].copy() @@ -271,10 +278,58 @@ def _calculate_overlaps(self): self.tileRA[under_sel], self.tileDEC[under_sel], self.tileRA[over_sel], self.tileDEC[over_sel], tile_diameter) + def _reduced_bgs_footprint(self, tiles, bgs_footprint): + ''' impose reduced footprint for BGS. Depending on the specified number + of square degrees in bgs_footprint, remove BRIGHT tiles close to the + ecliptic while keeping NGC a contiguous piece. + ''' + from astropy.coordinates import SkyCoord + assert bgs_footprint < 14000. + + # bgs tiles + is_bgs = (tiles['PROGRAM'] == 'BRIGHT') + + # get ecliptic coordinates of tiles + tile_coord = SkyCoord(ra=tiles['RA'] * u.deg, dec=tiles['DEC'] * u.deg, frame='icrs') + ecl_lat = tile_coord.barycentrictrueecliptic.lat.to(u.deg).value + + is_ngc = (tile_coord.galactic.b.value >= 0) + is_sgc = (tile_coord.galactic.b.value < 0) + + # reduced BGS footprint based on specified sq. deg. + if bgs_footprint == 13000: + bgsfoot = ( + (is_bgs & is_ngc) | + (is_bgs & is_sgc & (np.abs(tile_coord.barycentrictrueecliptic.lat.to(u.deg).value) > 6.5))) + elif bgs_footprint == 12000: + bgsfoot = ( + (is_bgs & is_ngc & (tile_coord.barycentrictrueecliptic.lat.to(u.deg).value > 8.1)) | + (is_bgs & is_sgc)) + elif bgs_footprint == 11000: + bgsfoot = ( + (is_bgs & is_ngc & (tile_coord.barycentrictrueecliptic.lat.to(u.deg).value > 8.1)) | + (is_bgs & is_sgc & (np.abs(tile_coord.barycentrictrueecliptic.lat.to(u.deg).value) > 6.5))) + elif bgs_footprint == 10000: + bgsfoot = ( + (is_bgs & is_ngc & (tile_coord.barycentrictrueecliptic.lat.to(u.deg).value > 12.6)) | + (is_bgs & is_sgc & (np.abs(tile_coord.barycentrictrueecliptic.lat.to(u.deg).value) > 9.7))) + else: + raise NotImplementedError + + # only keep non BRIGHT tiles or BRIGHT tiles far from ecliptic + reduced_footprint = (~is_bgs) | bgsfoot + + print(' reduced BGS footprint (v2)') + print(' %i of %i BRIGHT tiles removed for reduced footprint' % + (np.sum(is_bgs & ~bgsfoot), np.sum(is_bgs))) + print(' ~%f sq.deg' % (np.sum(bgsfoot)/np.sum(is_bgs) * 14000.)) + + return tiles[reduced_footprint], reduced_footprint + _cached_tiles = {} -def get_tiles(tiles_file=None, use_cache=True, write_cache=True): +def get_tiles(tiles_file=None, use_cache=True, write_cache=True, bgs_footprint=None): """Return a Tiles object with optional caching. You should normally always use the default arguments to ensure @@ -301,13 +356,21 @@ def get_tiles(tiles_file=None, use_cache=True, write_cache=True): if use_cache and tiles_file in _cached_tiles: tiles = _cached_tiles[tiles_file] log.debug('Using cached tiles for "{}".'.format(tiles_file)) + print(' Using cached tiles for "{}".'.format(tiles_file)) + for pname in Tiles.PROGRAMS: + pinfo = [] + for passnum in tiles.program_passes[pname]: + pinfo.append('{}({})'.format(passnum, tiles.pass_ntiles[passnum])) + print('{:6s} passes(tiles): {}.'.format(pname, ', '.join(pinfo))) else: - tiles = Tiles(tiles_file) + tiles = Tiles(tiles_file, bgs_footprint=bgs_footprint) log.info('Initialized tiles from "{}".'.format(tiles_file)) + print(' Initialized tiles from "{}".'.format(tiles_file)) for pname in Tiles.PROGRAMS: pinfo = [] for passnum in tiles.program_passes[pname]: pinfo.append('{}({})'.format(passnum, tiles.pass_ntiles[passnum])) + print('{:6s} passes(tiles): {}.'.format(pname, ', '.join(pinfo))) log.info('{:6s} passes(tiles): {}.'.format(pname, ', '.join(pinfo))) if write_cache: