From b35e0acb895a26ffd14c23c59b65c717b6ea2f02 Mon Sep 17 00:00:00 2001 From: changhoonhahn Date: Wed, 24 Feb 2021 06:01:45 -0800 Subject: [PATCH 01/12] added sun to objects we need to avoid --- py/desisurvey/ephem.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/py/desisurvey/ephem.py b/py/desisurvey/ephem.py index caf47b45..99ca98dd 100644 --- a/py/desisurvey/ephem.py +++ b/py/desisurvey/ephem.py @@ -195,7 +195,7 @@ 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: + for name in 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', From 5a4e5842f84a4928e8d2a108542e2662c960ae39 Mon Sep 17 00:00:00 2001 From: changhoonhahn Date: Wed, 24 Feb 2021 06:18:57 -0800 Subject: [PATCH 02/12] `etc.py`: added functions to calculate exposure factor during bright time --- py/desisurvey/etc.py | 196 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 196 insertions(+) diff --git a/py/desisurvey/etc.py b/py/desisurvey/etc.py index d3b43170..4adf4cd6 100644 --- a/py/desisurvey/etc.py +++ b/py/desisurvey/etc.py @@ -12,6 +12,8 @@ import numpy as np +from itertools import chain, combinations_with_replacement + import astropy.units as u import specsim.atmosphere @@ -208,6 +210,200 @@ 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 for bright time exposures. + + This factor is based on fits to sky brightness measurements from DESI SV1, + DESI CMX, and BOSS. + + TODO: + - update twilight component, which is only fit using BOSS twilight + + For details, see the following jupyter notebooks: + - https://github.com/desi-bgs/bgs-cmxsv/blob/55850e44d65570da69de0788c652cff698416834/doc/nb/sv1_sky_model_fit.ipynb + - https://github.com/desi-bgs/bgs-cmxsv/blob/55850e44d65570da69de0788c652cff698416834/doc/nb/sky_model_twilight_fit.ipynb + + Parameters + ---------- + airmass : array + Array of airmass used for observing this tile, must be >= 1. + moon_frac : float + Illuminated fraction of the moon within range [0,1]. + moon_sep : array + Array of separation angles between field center and moon in degrees within the + range [0,180]. + moon_alt : float + Altitude angle of the moon above the horizon in degrees within range [-90,90]. + sun_sep : array + Arry of separations angles between field center and sun in degrees + sun_alt : float + Altitude angle of the sunin degrees + + Returns + ------- + array + Dimensionless factors that exposure time should be increased to + account for increased sky brightness during bright time. + """ + # 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 + # for details see https://github.com/desi-bgs/bgs-cmxsv/blob/55850e44d65570da69de0788c652cff698416834/doc/nb/new_bgs_reference_sky.ipynb + 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): + ''' Calculate sky surface brightness at 5000A during bright time without + twilight. + + Surface brightness is based on a polynomial regression model that 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 + + For details, see jupyter notebook: + https://github.com/desi-bgs/bgs-cmxsv/blob/4c5f124164b649c595cd2dca87d14ba9f3b2c64d/doc/nb/sv1_sky_model_fit.ipynb + + Parameters + ---------- + airmass : array + Array of airmass used for observing this tile, must be >= 1. + moon_frac : array + Array of illuminated fraction of the moon within range [0,1]. + moon_sep : array + Array of separation angles between field center and moon in degrees within the + range [0,180]. + moon_alt : array + Array of altitude angle of the moon above the horizon in degrees within + range [-90,90]. + + Returns + ------- + array + Array of sky surface brightness at 5000A for bright time without + twilight + ''' + # 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): + ''' Calculate twilight contribution to sky surface brightness at 5000A + during bright time. + + Surface brightness is based on a polynomial regression model that was fit + to twilight contributions measured from BOSS exposures with sun altitutde > + -18deg. + + TODO: + * update fit using SV1 twilight exposures. + + + For details, see jupyter notebook: + https://github.com/desi-bgs/bgs-cmxsv/blob/55850e44d65570da69de0788c652cff698416834/doc/nb/sky_model_twilight_fit.ipynb + + Parameters + ---------- + airmass : array + Array of airmass used for observing this tile, must be >= 1. + sun_sep : array + Arry of separations angles between field center and sun in degrees + sun_alt : float + Altitude angle of the sunin degrees + + Returns + ------- + array + Array of twilight contribution to sky surface brightness at 5000A for + bright time. + ''' + 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. From d6a3232ebebfd39b31b01edaec8117c675d3368b Mon Sep 17 00:00:00 2001 From: changhoonhahn Date: Wed, 24 Feb 2021 06:19:58 -0800 Subject: [PATCH 03/12] `etc.py`: removed moon_up_factor for BRIGHT since the exposure factor will be calculated based on observing conditions --- py/desisurvey/etc.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/py/desisurvey/etc.py b/py/desisurvey/etc.py index 4adf4cd6..e2f08f89 100644 --- a/py/desisurvey/etc.py +++ b/py/desisurvey/etc.py @@ -513,12 +513,11 @@ def __init__(self, save_history=False): 'using default exposure time of 1000 s') moon_up_factor = getattr(config, 'moon_up_factor', None) if moon_up_factor: - for cond in ['DARK', 'GRAY', 'BRIGHT']: + for cond in ['DARK', 'GRAY']: self.TEXP_TOTAL[cond] *= getattr(moon_up_factor, cond)() else: # Temporary hardcoded exposure factors for moon-up observing. 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]) From 54e87570ae2a2788004271b15f80c6c7b6d2b627 Mon Sep 17 00:00:00 2001 From: changhoonhahn Date: Wed, 24 Feb 2021 06:22:56 -0800 Subject: [PATCH 04/12] `scheduler.py`: tracks moon_ALTAZ and sun_DECRA, sun_ALTAZ --- py/desisurvey/scheduler.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/py/desisurvey/scheduler.py b/py/desisurvey/scheduler.py index 4ab095fa..fc62a900 100644 --- a/py/desisurvey/scheduler.py +++ b/py/desisurvey/scheduler.py @@ -261,7 +261,10 @@ 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 select_program(self, mjd_now, ETC, verbose=False): """Select program to observe now. From b45d725f8c37e84f40b94c0f29bf3af5e87276a6 Mon Sep 17 00:00:00 2001 From: changhoonhahn Date: Wed, 24 Feb 2021 06:34:13 -0800 Subject: [PATCH 05/12] `scheduler.next_tile`: exposure factor for bright time included in exposure factor calculation for available tiles. --- py/desisurvey/scheduler.py | 26 +++++++++++++++++++++++++- 1 file changed, 25 insertions(+), 1 deletion(-) diff --git a/py/desisurvey/scheduler.py b/py/desisurvey/scheduler.py index fc62a900..abe83c43 100644 --- a/py/desisurvey/scheduler.py +++ b/py/desisurvey/scheduler.py @@ -418,7 +418,31 @@ def next_tile(self, mjd_now, ETC, seeing, transp, skylevel, HA_sigma=15., greedi # 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 program == 'BRIGHT': + # Calculate the moon (RA,DEC). + moonDEC, moonRA = self.moon_DECRA(mjd_now) + # Calculate the moon (ALT,AZ). + moonALT, moonAZ = self.moon_ALTAZ(mjd_now) + # Calculate the sun (RA,DEC). + sunDEC, sunRA = self.sun_DECRA(mjd_now) + # Calculate the sun (ALT,AZ). + sunALT, sunAZ = self.sun_ALTAZ(mjd_now) + # moon illumination + moonILL = self.night_ephem['moon_illum_frac'] + # Calculate moon separation + moonSEP = desisurvey.utils.separation_matrix( + [moonRA], [moonDEC], + self.tiles.tileRA[self.tile_sel], self.tiles.tileDEC[self.tile_sel]) + # Calculate sun separation + sunSEP = desisurvey.utils.separation_matrix( + [sunRA], [sunDEC], + self.tiles.tileRA[self.tile_sel], self.tiles.tileDEC[self.tile_sel]) + + self.exposure_factor[self.tile_sel] *= \ + desisurvey.etc.bright_exposure_factor( + self.airmass[self.tile_sel], moonILL, moonSEP, moonALT, sunSEP, sunALT) + 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): From ab3980168f8261e420cd753e5eeb24197628b01f Mon Sep 17 00:00:00 2001 From: changhoonhahn Date: Wed, 24 Feb 2021 12:55:38 -0800 Subject: [PATCH 06/12] erfa import issue fix --- py/desisurvey/ephem.py | 2 +- py/desisurvey/utils.py | 1 - 2 files changed, 1 insertion(+), 2 deletions(-) diff --git a/py/desisurvey/ephem.py b/py/desisurvey/ephem.py index 99ca98dd..b18bab06 100644 --- a/py/desisurvey/ephem.py +++ b/py/desisurvey/ephem.py @@ -195,7 +195,7 @@ 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 + ['sun']: + 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/utils.py b/py/desisurvey/utils.py index d8a4fd0f..b32272dd 100644 --- a/py/desisurvey/utils.py +++ b/py/desisurvey/utils.py @@ -15,7 +15,6 @@ import astropy.utils.iers import astropy.utils.data import astropy.utils.exceptions -import astropy._erfa.core import astropy.units as u import desiutil.log From ef500a650b069dec33644525e56fd8b2197dd14e Mon Sep 17 00:00:00 2001 From: changhoonhahn Date: Thu, 4 Mar 2021 12:50:01 -0800 Subject: [PATCH 07/12] updated sky model and twilight model --- py/desisurvey/etc.py | 118 ++++++++++++++++++++++++++----------------- 1 file changed, 71 insertions(+), 47 deletions(-) diff --git a/py/desisurvey/etc.py b/py/desisurvey/etc.py index e2f08f89..39c0265c 100644 --- a/py/desisurvey/etc.py +++ b/py/desisurvey/etc.py @@ -210,14 +210,13 @@ 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 for bright time exposures. - - This factor is based on fits to sky brightness measurements from DESI SV1, - DESI CMX, and BOSS. - - TODO: - - update twilight component, which is only fit using BOSS twilight +def sky_exposure_factor(tref, airmass, moon_frac, moon_sep, moon_alt, sun_sep, + sun_alt): + """ Calculate exposure time scaling factor due to sky brightness. This is + the factor that exposure time needs to be scaled by in order to get the + same SNR as an exposure during nominal dark sky. It's calculated based on a + sky brightness model fit to observed sky brightness from DESI SV1, DESI + CMX, and BOSS. For details, see the following jupyter notebooks: - https://github.com/desi-bgs/bgs-cmxsv/blob/55850e44d65570da69de0788c652cff698416834/doc/nb/sv1_sky_model_fit.ipynb @@ -225,6 +224,8 @@ def bright_exposure_factor(airmass, moon_frac, moon_sep, moon_alt, sun_sep, sun_ Parameters ---------- + tnom : array + Array of nominal exposure time airmass : array Array of airmass used for observing this tile, must be >= 1. moon_frac : float @@ -245,10 +246,6 @@ def bright_exposure_factor(airmass, moon_frac, moon_sep, moon_alt, sun_sep, sun_ Dimensionless factors that exposure time should be increased to account for increased sky brightness during bright time. """ - # 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() @@ -270,15 +267,14 @@ def bright_exposure_factor(airmass, moon_frac, moon_sep, moon_alt, sun_sep, sun_ # calculate exposure factor config = desisurvey.config.Configuration() - # sky brightness at 5000A for reference BGS exposure - # for details see https://github.com/desi-bgs/bgs-cmxsv/blob/55850e44d65570da69de0788c652cff698416834/doc/nb/new_bgs_reference_sky.ipynb - Isky5000_ref = 3.6956611461286966 # 1e-17 erg/s/cm^2/A/arcsec^2 + # sky brightness at 5000A for nominal dark sky + Isky5000_ref = 1.1282850428182252 # 1e-17 erg/s/cm^2/A/arcsec^2 - # exposure time for reference BGS exposure + # exposure time for nominal 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( + Isky5000_exp = Isky5000_notwilight_regression( airmass, np.repeat(moon_frac, nexp), moon_sep, @@ -286,7 +282,7 @@ def bright_exposure_factor(airmass, moon_frac, moon_sep, moon_alt, sun_sep, sun_ # add twilight contribution if sun_alt >= -18.: - Isky5000_exp += np.clip(bright_Isky5000_twilight_regression( + Isky5000_exp += np.clip(Isky5000_twilight_regression( airmass, sun_sep, np.repeat(sun_alt, nexp)), 0, None) @@ -310,18 +306,19 @@ def bright_exposure_factor(airmass, moon_frac, moon_sep, moon_alt, sun_sep, sun_ return texp/tref -def bright_Isky5000_notwilight_regression(airmass, moon_frac, moon_sep, moon_alt): +def Isky5000_notwilight_regression(airmass, moon_frac, moon_sep, moon_alt): ''' Calculate sky surface brightness at 5000A during bright time without twilight. - Surface brightness is based on a polynomial regression model that was fit - using observed sky surface brightnesses from - * DESI SV1 bright exposures - * DESI CMX exposures with transparency > 0.9 + Surface brightness is based on a polynomial regression model fit to + log(observed sky brightness). The observed sky brightnesses were compiled + from + * DESI SV1 bright exposures with TRANSP > 0.9 and SUN_ALT < -18 + * DESI CMX exposures with transparency > 0.95 * BOSS exposures with sun alt < -18 For details, see jupyter notebook: - https://github.com/desi-bgs/bgs-cmxsv/blob/4c5f124164b649c595cd2dca87d14ba9f3b2c64d/doc/nb/sv1_sky_model_fit.ipynb + https://github.com/desi-bgs/bgs-cmxsv/blob/521211dccff7cb3b71b07d84522802acade26071/doc/nb/sv1_sky_model_fit.ipynb Parameters ---------- @@ -340,41 +337,66 @@ def bright_Isky5000_notwilight_regression(airmass, moon_frac, moon_sep, moon_alt ------- array Array of sky surface brightness at 5000A for bright time without - twilight + twilight in units of erg/s/cm^2/A/arcsec^2 ''' + # polynomial of order + norder = 5 # 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]) - + coeffs = np.array([ + 1.19728287e-02, 2.68146174e-02, 2.78840941e-02, 4.41723544e-03, + -1.02111676e-01, 7.65110151e-03, 1.63878220e-02, -3.08958375e-02, + 3.85657345e-01, 6.83474920e-03, 3.06045766e-01, 2.55798622e-01, + 1.19306254e-03, -1.04242412e-02, -3.34780689e-03, -1.68835833e-02, + 1.11378109e-02, -1.75611187e-01, -2.38751513e-01, 1.46593863e-02, + -4.36980545e-02, -1.78128533e-01, 5.85014404e-03, 1.37745149e-02, + -1.81392009e-03, -1.27115538e-02, 9.35953790e-02, -1.52422946e-01, + -1.74469624e-04, -6.01814594e-03, -2.13820870e-04, -7.40115962e-05, + -4.43960687e-05, 1.31544250e-04, 5.23119315e-05, -3.07378710e-02, + 3.50689477e-02, 1.85015093e-01, 6.13788213e-02, 1.98622649e-02, + -7.03566704e-02, -2.70010457e-01, -4.13256930e-03, -4.65235557e-03, + 2.08199644e-03, 4.06645589e-03, 2.08415057e-01, 4.87587825e-01, + -7.34996335e-03, 1.37774943e-03, 6.26446800e-03, 2.14264776e-05, + -1.34520777e-05, -1.28680576e-04, -2.98020923e-05, 1.41853618e-03, + -2.22111586e-01, -5.77624465e-01, -1.09989601e-02, -2.77969195e-04, + 2.17332386e-03, 1.47785979e-04, 1.63001051e-04, 3.10888811e-05, + -5.98993375e-05, -8.54879999e-08, 3.29709501e-07, 2.92447174e-08, + -2.79763418e-07, -2.38343236e-08, -1.10769335e-01, 8.24805309e-02, + -2.67759016e-02, -1.09365333e-03, 9.87385635e-03, 3.12620536e-02, + 9.15874867e-02, -2.52374361e-05, -4.92100655e-04, -4.21467759e-04, + 5.28695530e-03, 1.93471362e-01, -5.37647583e-02, -4.34387094e-04, + -2.25656833e-03, -7.15039182e-04, 1.18896366e-05, 2.77215519e-05, + 2.81214300e-05, -2.48811976e-06, 1.29700918e-02, -1.92447646e-02, + -2.44212957e-01, -5.99603226e-03, -2.76089664e-03, 9.14612070e-04, + 6.91296852e-05, 1.02360578e-04, 1.16371717e-05, -2.52971480e-05, + -1.37227312e-07, -3.75516901e-07, 5.48164999e-08, 5.06994020e-07, + 1.93984829e-07, 5.32931890e-02, 5.79701139e-01, 4.34770814e-01, + -6.79348983e-03, -1.31426633e-02, -6.14432600e-05, 1.39891173e-04, + 3.29198340e-04, 1.11812563e-04, -1.17664897e-05, -1.29847018e-06, + -3.91089212e-06, -2.85449951e-06, -3.27679136e-07, 3.48519396e-07, + 3.29795490e-09, 7.65512533e-09, 6.37102357e-09, 5.95568585e-10, + -2.63641980e-09, -1.30434882e-09]) + 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)) + combs = chain.from_iterable(combinations_with_replacement(range(4), i) for i in range(0, norder+1)) 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) + return np.exp(np.dot(theta_transform, coeffs.T)) -def bright_Isky5000_twilight_regression(airmass, sun_sep, sun_alt): +def Isky5000_twilight_regression(airmass, sun_sep, sun_alt): ''' Calculate twilight contribution to sky surface brightness at 5000A during bright time. Surface brightness is based on a polynomial regression model that was fit - to twilight contributions measured from BOSS exposures with sun altitutde > - -18deg. - - TODO: - * update fit using SV1 twilight exposures. - - - For details, see jupyter notebook: - https://github.com/desi-bgs/bgs-cmxsv/blob/55850e44d65570da69de0788c652cff698416834/doc/nb/sky_model_twilight_fit.ipynb + to twilight contributions measured from SV1 exposures with sun alt > -18 + and TRANSP > 0.95 and BOSS exposures with sun altitutde > -18deg. For + details, see jupyter notebook: + https://github.com/desi-bgs/bgs-cmxsv/blob/521211dccff7cb3b71b07d84522802acade26071/doc/nb/sky_model_twilight_fit.ipynb Parameters ---------- @@ -383,16 +405,18 @@ def bright_Isky5000_twilight_regression(airmass, sun_sep, sun_alt): sun_sep : array Arry of separations angles between field center and sun in degrees sun_alt : float - Altitude angle of the sunin degrees + Altitude angle of the sun in degrees Returns ------- array Array of twilight contribution to sky surface brightness at 5000A for - bright time. + bright time in units of erg/s/cm^2/A/arcsec^2 ''' norder = 1 - skymodel_coeff = np.array([2.00474700e+00, 3.19827604e+00, 3.13212960e-01, 2.69673079e-03]) + # coefficients fit in notebook + # https://github.com/desi-bgs/bgs-cmxsv/blob/521211dccff7cb3b71b07d84522802acade26071/doc/nb/sky_model_twilight_fit.ipynb + skymodel_coeff = np.array([3.50300842, 2.63453598, 0.40491282, 0.0096106 ]) theta = np.atleast_2d(np.array([airmass, sun_alt, sun_sep]).T) @@ -401,7 +425,7 @@ def bright_Isky5000_twilight_regression(airmass, sun_sep, sun_alt): for i, comb in enumerate(combs): theta_transform[:, i] = theta[:, comb].prod(1) - return np.dot(theta_transform, skymodel_coeff.T) + return np.clip(np.dot(theta_transform, skymodel_coeff.T), 0, None) def exposure_time(program, seeing, transparency, airmass, EBV, From 015a5fbb1fc33448720c786fc3dce69582eaf149 Mon Sep 17 00:00:00 2001 From: changhoonhahn Date: Fri, 5 Mar 2021 14:43:31 -0800 Subject: [PATCH 08/12] exposure factor from sky brightness is now included for DARK and GRAY times; etc.sky_exposure_factor modified to take nominal exposure time as an input --- py/desisurvey/etc.py | 19 +++++-------- py/desisurvey/scheduler.py | 57 +++++++++++++++++++++----------------- 2 files changed, 39 insertions(+), 37 deletions(-) diff --git a/py/desisurvey/etc.py b/py/desisurvey/etc.py index 39c0265c..2a811ac9 100644 --- a/py/desisurvey/etc.py +++ b/py/desisurvey/etc.py @@ -213,9 +213,9 @@ def moon_exposure_factor(moon_frac, moon_sep, moon_alt, airmass): def sky_exposure_factor(tref, airmass, moon_frac, moon_sep, moon_alt, sun_sep, sun_alt): """ Calculate exposure time scaling factor due to sky brightness. This is - the factor that exposure time needs to be scaled by in order to get the - same SNR as an exposure during nominal dark sky. It's calculated based on a - sky brightness model fit to observed sky brightness from DESI SV1, DESI + the factor that nominal exposure time needs to be scaled by in order to get + the same SNR as an exposure during nominal dark sky. It's calculated based + on a sky brightness model fit to observed sky brightness from DESI SV1, DESI CMX, and BOSS. For details, see the following jupyter notebooks: @@ -225,7 +225,8 @@ def sky_exposure_factor(tref, airmass, moon_frac, moon_sep, moon_alt, sun_sep, Parameters ---------- tnom : array - Array of nominal exposure time + Array of *nominal* exposure time in seconds. This is an input parameter + because different programs have different nominal exposure times. airmass : array Array of airmass used for observing this tile, must be >= 1. moon_frac : float @@ -264,16 +265,10 @@ def sky_exposure_factor(tref, airmass, moon_frac, moon_sep, moon_alt, sun_sep, assert len(moon_sep) == nexp assert len(sun_sep) == nexp - # calculate exposure factor - config = desisurvey.config.Configuration() - # sky brightness at 5000A for nominal dark sky Isky5000_ref = 1.1282850428182252 # 1e-17 erg/s/cm^2/A/arcsec^2 - # exposure time for nominal exposure - tref = getattr(config.nominal_exposure_time, 'BRIGHT')().to(u.s).value - - # sky brightness at 5000A for observing conditions + # sky brightness without twilight at 5000A for observing conditions Isky5000_exp = Isky5000_notwilight_regression( airmass, np.repeat(moon_frac, nexp), @@ -313,7 +308,7 @@ def Isky5000_notwilight_regression(airmass, moon_frac, moon_sep, moon_alt): Surface brightness is based on a polynomial regression model fit to log(observed sky brightness). The observed sky brightnesses were compiled from - * DESI SV1 bright exposures with TRANSP > 0.9 and SUN_ALT < -18 + * DESI SV1 bright exposures with TRANSP > 0.95 and SUN_ALT < -18 * DESI CMX exposures with transparency > 0.95 * BOSS exposures with sun alt < -18 diff --git a/py/desisurvey/scheduler.py b/py/desisurvey/scheduler.py index abe83c43..44449bde 100644 --- a/py/desisurvey/scheduler.py +++ b/py/desisurvey/scheduler.py @@ -133,6 +133,10 @@ def __init__(self, restore=None, design_hourangle=None): self.tile_available = np.zeros(self.tiles.ntiles, bool) self.tile_planned = np.zeros(self.tiles.ntiles, bool) self.tile_priority = np.zeros(self.tiles.ntiles, float) + self.nominal_exposure_time_sec = ( + desisurvey.tiles.get_nominal_program_times(self.tiles.tileprogram)) + self.maxtime = config.maxtime() + # Lookup avoidance cone angles. self.avoid_bodies = {} for body in config.avoid_bodies.keys: @@ -415,36 +419,39 @@ 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 + + # Calculate the moon (RA,DEC). + moonDEC, moonRA = self.moon_DECRA(mjd_now) + # Calculate the moon (ALT,AZ). + moonALT, moonAZ = self.moon_ALTAZ(mjd_now) + # Calculate the sun (RA,DEC). + sunDEC, sunRA = self.sun_DECRA(mjd_now) + # Calculate the sun (ALT,AZ). + sunALT, sunAZ = self.sun_ALTAZ(mjd_now) + # moon illumination + moonILL = self.night_ephem['moon_illum_frac'] + # Calculate moon separation + moonSEP = desisurvey.utils.separation_matrix( + [moonRA], [moonDEC], + self.tiles.tileRA[self.tile_sel], self.tiles.tileDEC[self.tile_sel]) + # Calculate sun separation + sunSEP = desisurvey.utils.separation_matrix( + [sunRA], [sunDEC], + self.tiles.tileRA[self.tile_sel], self.tiles.tileDEC[self.tile_sel]) + # Estimate exposure factors for all available tiles. self.exposure_factor[:] = 1e8 self.exposure_factor[self.tile_sel] = self.tiles.dust_factor[self.tile_sel] - if program == 'BRIGHT': - # Calculate the moon (RA,DEC). - moonDEC, moonRA = self.moon_DECRA(mjd_now) - # Calculate the moon (ALT,AZ). - moonALT, moonAZ = self.moon_ALTAZ(mjd_now) - # Calculate the sun (RA,DEC). - sunDEC, sunRA = self.sun_DECRA(mjd_now) - # Calculate the sun (ALT,AZ). - sunALT, sunAZ = self.sun_ALTAZ(mjd_now) - # moon illumination - moonILL = self.night_ephem['moon_illum_frac'] - # Calculate moon separation - moonSEP = desisurvey.utils.separation_matrix( - [moonRA], [moonDEC], - self.tiles.tileRA[self.tile_sel], self.tiles.tileDEC[self.tile_sel]) - # Calculate sun separation - sunSEP = desisurvey.utils.separation_matrix( - [sunRA], [sunDEC], - self.tiles.tileRA[self.tile_sel], self.tiles.tileDEC[self.tile_sel]) - - self.exposure_factor[self.tile_sel] *= \ - desisurvey.etc.bright_exposure_factor( - self.airmass[self.tile_sel], moonILL, moonSEP, moonALT, sunSEP, sunALT) - else: - self.exposure_factor[self.tile_sel] *= desisurvey.etc.airmass_exposure_factor(self.airmass[self.tile_sel]) + self.exposure_factor[self.tile_sel] *= \ + desisurvey.etc.sky_exposure_factor( + self.nominal_exposure_time_sec, # nominal exposure time + self.airmass[self.tile_sel], + moonILL, moonSEP, moonALT, + sunSEP, sunALT) + 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 From 0e39986e15a6bc9810e0e978852a29c4a9fd759b Mon Sep 17 00:00:00 2001 From: changhoonhahn Date: Wed, 24 Feb 2021 06:01:45 -0800 Subject: [PATCH 09/12] added sun to objects we need to avoid --- py/desisurvey/ephem.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/py/desisurvey/ephem.py b/py/desisurvey/ephem.py index 3873bec5..4492e52b 100644 --- a/py/desisurvey/ephem.py +++ b/py/desisurvey/ephem.py @@ -195,7 +195,7 @@ 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: + for name in 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', From 5fd025dbad6d1b858aa45081fde6640e60ba46c1 Mon Sep 17 00:00:00 2001 From: changhoonhahn Date: Wed, 24 Feb 2021 06:18:57 -0800 Subject: [PATCH 10/12] `etc.py`: added functions to calculate exposure factor during bright time --- py/desisurvey/etc.py | 196 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 196 insertions(+) diff --git a/py/desisurvey/etc.py b/py/desisurvey/etc.py index a449f9d2..9c72c421 100644 --- a/py/desisurvey/etc.py +++ b/py/desisurvey/etc.py @@ -12,6 +12,8 @@ import numpy as np +from itertools import chain, combinations_with_replacement + import astropy.units as u import specsim.atmosphere @@ -214,6 +216,200 @@ 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 for bright time exposures. + + This factor is based on fits to sky brightness measurements from DESI SV1, + DESI CMX, and BOSS. + + TODO: + - update twilight component, which is only fit using BOSS twilight + + For details, see the following jupyter notebooks: + - https://github.com/desi-bgs/bgs-cmxsv/blob/55850e44d65570da69de0788c652cff698416834/doc/nb/sv1_sky_model_fit.ipynb + - https://github.com/desi-bgs/bgs-cmxsv/blob/55850e44d65570da69de0788c652cff698416834/doc/nb/sky_model_twilight_fit.ipynb + + Parameters + ---------- + airmass : array + Array of airmass used for observing this tile, must be >= 1. + moon_frac : float + Illuminated fraction of the moon within range [0,1]. + moon_sep : array + Array of separation angles between field center and moon in degrees within the + range [0,180]. + moon_alt : float + Altitude angle of the moon above the horizon in degrees within range [-90,90]. + sun_sep : array + Arry of separations angles between field center and sun in degrees + sun_alt : float + Altitude angle of the sunin degrees + + Returns + ------- + array + Dimensionless factors that exposure time should be increased to + account for increased sky brightness during bright time. + """ + # 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 + # for details see https://github.com/desi-bgs/bgs-cmxsv/blob/55850e44d65570da69de0788c652cff698416834/doc/nb/new_bgs_reference_sky.ipynb + 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): + ''' Calculate sky surface brightness at 5000A during bright time without + twilight. + + Surface brightness is based on a polynomial regression model that 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 + + For details, see jupyter notebook: + https://github.com/desi-bgs/bgs-cmxsv/blob/4c5f124164b649c595cd2dca87d14ba9f3b2c64d/doc/nb/sv1_sky_model_fit.ipynb + + Parameters + ---------- + airmass : array + Array of airmass used for observing this tile, must be >= 1. + moon_frac : array + Array of illuminated fraction of the moon within range [0,1]. + moon_sep : array + Array of separation angles between field center and moon in degrees within the + range [0,180]. + moon_alt : array + Array of altitude angle of the moon above the horizon in degrees within + range [-90,90]. + + Returns + ------- + array + Array of sky surface brightness at 5000A for bright time without + twilight + ''' + # 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): + ''' Calculate twilight contribution to sky surface brightness at 5000A + during bright time. + + Surface brightness is based on a polynomial regression model that was fit + to twilight contributions measured from BOSS exposures with sun altitutde > + -18deg. + + TODO: + * update fit using SV1 twilight exposures. + + + For details, see jupyter notebook: + https://github.com/desi-bgs/bgs-cmxsv/blob/55850e44d65570da69de0788c652cff698416834/doc/nb/sky_model_twilight_fit.ipynb + + Parameters + ---------- + airmass : array + Array of airmass used for observing this tile, must be >= 1. + sun_sep : array + Arry of separations angles between field center and sun in degrees + sun_alt : float + Altitude angle of the sunin degrees + + Returns + ------- + array + Array of twilight contribution to sky surface brightness at 5000A for + bright time. + ''' + 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. From 467afcfe44f994b4aa1a98f7ec33fb0e05ea9783 Mon Sep 17 00:00:00 2001 From: changhoonhahn Date: Thu, 18 Mar 2021 20:15:42 -0700 Subject: [PATCH 11/12] sky_level function added to etc module. this function is called by get_weather in surveysim.nightops --- py/desisurvey/etc.py | 149 ++++++++++++++++++------------------- py/desisurvey/scheduler.py | 5 +- 2 files changed, 72 insertions(+), 82 deletions(-) diff --git a/py/desisurvey/etc.py b/py/desisurvey/etc.py index ad405973..a7128b9c 100644 --- a/py/desisurvey/etc.py +++ b/py/desisurvey/etc.py @@ -14,7 +14,9 @@ from itertools import chain, combinations_with_replacement +import astropy.time import astropy.units as u +from astropy.coordinates import SkyCoord import specsim.atmosphere @@ -22,6 +24,7 @@ import desisurvey.config import desisurvey.tiles +import desisurvey.ephem def seeing_exposure_factor(seeing): @@ -216,95 +219,85 @@ def moon_exposure_factor(moon_frac, moon_sep, moon_alt, airmass): return _moonCoefficients.dot(X) -def sky_exposure_factor(tref, airmass, moon_frac, moon_sep, moon_alt, sun_sep, - sun_alt): - """ Calculate exposure time scaling factor due to sky brightness. This is - the factor that nominal exposure time needs to be scaled by in order to get - the same SNR as an exposure during nominal dark sky. It's calculated based - on a sky brightness model fit to observed sky brightness from DESI SV1, DESI - CMX, and BOSS. - - For details, see the following jupyter notebooks: - - https://github.com/desi-bgs/bgs-cmxsv/blob/55850e44d65570da69de0788c652cff698416834/doc/nb/sv1_sky_model_fit.ipynb - - https://github.com/desi-bgs/bgs-cmxsv/blob/55850e44d65570da69de0788c652cff698416834/doc/nb/sky_model_twilight_fit.ipynb +def sky_level(mjd, ra, dec): + ''' Calculate the sky surface brightness given MJD, RA, and Dec Parameters ---------- - tnom : array - Array of *nominal* exposure time in seconds. This is an input parameter - because different programs have different nominal exposure times. - airmass : array - Array of airmass used for observing this tile, must be >= 1. - moon_frac : float - Illuminated fraction of the moon within range [0,1]. - moon_sep : array - Array of separation angles between field center and moon in degrees within the - range [0,180]. - moon_alt : float - Altitude angle of the moon above the horizon in degrees within range [-90,90]. - sun_sep : array - Arry of separations angles between field center and sun in degrees - sun_alt : float - Altitude angle of the sunin degrees + mjd : float - Returns - ------- - array - Dimensionless factors that exposure time should be increased to - account for increased sky brightness during bright time. - """ - # 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.') + ra : float + RA in deg - # check size of inputs - nexp = len(airmass) - assert len(moon_sep) == nexp - assert len(sun_sep) == nexp + dec : float + Dec in degrees deg - # sky brightness at 5000A for nominal dark sky - Isky5000_ref = 1.1282850428182252 # 1e-17 erg/s/cm^2/A/arcsec^2 - + Returns + ------- + float + sky surface brightness in units of 1e-17 erg/s/cm^2/A/arcsec^2 + ''' + ephem = desisurvey.ephem.get_ephem() + night = desisurvey.utils.get_date(mjd) + night_ephem = ephem.get_night(night) + + # get moon illumination + moon_ill = night_ephem['moon_illum_frac'] + + if ra is None or dec is None: + # default values when there's no tile position yet + frame = desisurvey.utils.get_observer(astropy.time.Time(mjd, format='mjd')) + # get RA and Dec at the zenith + altaz = SkyCoord(alt=90.*u.deg, az=0*u.deg, + obstime=astropy.time.Time(mjd, format='mjd'), + frame='altaz', + location=frame.location) + ra, dec = altaz.icrs.ra.to(u.deg).value, altaz.icrs.dec.to(u.deg).value + + # calculate airmass + airmass = desisurvey.utils.get_airmass( + astropy.time.Time(mjd, format='mjd'), + ra * u.deg, + dec * u.deg) + + # calculate moon altitude and separation + moon_DECRA = desisurvey.ephem.get_object_interpolator(night_ephem, 'moon', altaz=False) + moon_ALTAZ = desisurvey.ephem.get_object_interpolator(night_ephem, 'moon', altaz=True) + + moon_alt, _ = moon_ALTAZ(mjd) # moon altitude + moon_dec, moon_ra = moon_DECRA(mjd) + moon_sep = desisurvey.utils.separation_matrix( # separation + [moon_ra], [moon_dec], [ra], [dec])[0][0] + + # sun moon altitude and separation + sun_DECRA = desisurvey.ephem.get_object_interpolator(night_ephem, 'sun', altaz=False) + sun_ALTAZ = desisurvey.ephem.get_object_interpolator(night_ephem, 'sun', altaz=True) + + sun_alt, _ = sun_ALTAZ(mjd) # altitude + sun_dec, sun_ra = sun_DECRA(mjd) + sun_sep = desisurvey.utils.separation_matrix([sun_ra], [sun_dec], [ra], + [dec])[0][0] + + print(airmass, moon_ill,moon_sep, moon_alt, sun_sep, sun_alt) + + assert moon_sep > 30. - moon_alt + assert moon_sep < 160. - 1.1 * moon_alt + assert moon_sep < moon_alt + 250. + assert moon_sep > moon_alt - 50. # sky brightness without twilight at 5000A for observing conditions - Isky5000_exp = Isky5000_notwilight_regression( - airmass, - np.repeat(moon_frac, nexp), - moon_sep, - np.repeat(moon_alt, nexp)) + Isky5000_exp = Isky5000_notwilight_regression(airmass, moon_ill, moon_sep, moon_alt) # add twilight contribution if sun_alt >= -18.: - Isky5000_exp += np.clip(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 + Isky5000_exp += Isky5000_twilight_regression(airmass, sun_sep, sun_alt) - # (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 + ## sky fiber flux + #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_exp = fibflux5000_exp * 0.0629735016982807 + #return sky_photon_per_sec_exp + return Isky5000_exp def Isky5000_notwilight_regression(airmass, moon_frac, moon_sep, moon_alt): diff --git a/py/desisurvey/scheduler.py b/py/desisurvey/scheduler.py index 88a0b4b9..c6b9617b 100644 --- a/py/desisurvey/scheduler.py +++ b/py/desisurvey/scheduler.py @@ -158,10 +158,7 @@ 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) - # 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) + #self.moon_ALTAZ = desisurvey.ephem.get_object_interpolator(self.night_ephem, 'moon', altaz=True) def select_program(self, mjd_now, ETC, verbose=False): """Select program to observe now. From cab78618bcaaa534bfbaa71e0d4312e88603fa3e Mon Sep 17 00:00:00 2001 From: changhoonhahn Date: Fri, 7 May 2021 06:07:07 -0700 Subject: [PATCH 12/12] snr2frac accumulation now uses effective sky rather than sky, which more accurately accounts for read noise. I've also added some documentation on sky_level --- py/desisurvey/etc.py | 185 ++++++++++++++++++++++++------------------- 1 file changed, 104 insertions(+), 81 deletions(-) diff --git a/py/desisurvey/etc.py b/py/desisurvey/etc.py index 4c27d2fb..c3ea00d2 100644 --- a/py/desisurvey/etc.py +++ b/py/desisurvey/etc.py @@ -236,30 +236,49 @@ def moon_exposure_factor(moon_frac, moon_sep, moon_alt, airmass): return _moonCoefficients.dot(X) -def sky_level(mjd, ra, dec): - ''' Calculate the sky surface brightness given MJD, RA, and Dec +def sky_level(mjd, ra, dec, moon_ill=None, moon_DECRA=None, moon_ALTAZ=None, sun_DECRA=None, sun_ALTAZ=None): + ''' Calculate the sky level. Sky level is defined as sky surface brightness + for given MJD, RA, and Dec over nominal sky sufrace brightness Parameters ---------- mjd : float - + date ra : float RA in deg - dec : float Dec in degrees deg + moon_ill : float, optional + moon illumination + moon_DECRA : optional + moon Dec and RA interpolator object from desisurvey.ephem + moon_ALTAZ : optional + moon Alt and AZ interpolator object from desisurvey.ephem + sun_DECRA : optional + sun Dec and RA interpolator object from desisurvey.ephem + sun_ALTAZ : optional + sun Alt and AZ interpolator object from desisurvey.ephem Returns ------- - float - sky surface brightness in units of 1e-17 erg/s/cm^2/A/arcsec^2 + sky level : float + current sky brightness over nominal sky brightness ''' - ephem = desisurvey.ephem.get_ephem() - night = desisurvey.utils.get_date(mjd) - night_ephem = ephem.get_night(night) + if moon_DECRA is None or moon_ALTAZ is None or sun_DECRA is None or sun_ALTAZ is None: + ephem = desisurvey.ephem.get_ephem() + night = desisurvey.utils.get_date(mjd) + night_ephem = ephem.get_night(night) + + # calculate moon altitude and separation + moon_DECRA = desisurvey.ephem.get_object_interpolator(night_ephem, 'moon', altaz=False) + moon_ALTAZ = desisurvey.ephem.get_object_interpolator(night_ephem, 'moon', altaz=True) + + # sun moon altitude and separation + sun_DECRA = desisurvey.ephem.get_object_interpolator(night_ephem, 'sun', altaz=False) + sun_ALTAZ = desisurvey.ephem.get_object_interpolator(night_ephem, 'sun', altaz=True) - # get moon illumination - moon_ill = night_ephem['moon_illum_frac'] + # get moon illumination + moon_ill = night_ephem['moon_illum_frac'] if ra is None or dec is None: # default values when there's no tile position yet @@ -276,45 +295,29 @@ def sky_level(mjd, ra, dec): astropy.time.Time(mjd, format='mjd'), ra * u.deg, dec * u.deg) - - # calculate moon altitude and separation - moon_DECRA = desisurvey.ephem.get_object_interpolator(night_ephem, 'moon', altaz=False) - moon_ALTAZ = desisurvey.ephem.get_object_interpolator(night_ephem, 'moon', altaz=True) - + # moon ephem moon_alt, _ = moon_ALTAZ(mjd) # moon altitude moon_dec, moon_ra = moon_DECRA(mjd) moon_sep = desisurvey.utils.separation_matrix( # separation [moon_ra], [moon_dec], [ra], [dec])[0][0] - - # sun moon altitude and separation - sun_DECRA = desisurvey.ephem.get_object_interpolator(night_ephem, 'sun', altaz=False) - sun_ALTAZ = desisurvey.ephem.get_object_interpolator(night_ephem, 'sun', altaz=True) - + # sun ephem sun_alt, _ = sun_ALTAZ(mjd) # altitude sun_dec, sun_ra = sun_DECRA(mjd) sun_sep = desisurvey.utils.separation_matrix([sun_ra], [sun_dec], [ra], [dec])[0][0] - - print(airmass, moon_ill,moon_sep, moon_alt, sun_sep, sun_alt) - - assert moon_sep > 30. - moon_alt - assert moon_sep < 160. - 1.1 * moon_alt - assert moon_sep < moon_alt + 250. - assert moon_sep > moon_alt - 50. # sky brightness without twilight at 5000A for observing conditions Isky5000_exp = Isky5000_notwilight_regression(airmass, moon_ill, moon_sep, moon_alt) - # add twilight contribution if sun_alt >= -18.: Isky5000_exp += Isky5000_twilight_regression(airmass, sun_sep, sun_alt) + + Isky5000_nom = 1.1282850428182252 # 1e-17 erg/s/cm^2/A/arcsec2 - ## sky fiber flux - #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_exp = fibflux5000_exp * 0.0629735016982807 - #return sky_photon_per_sec_exp - return Isky5000_exp + #if moon_sep < 30. - moon_alt: print('moon sep < 30 - moon alt', Isky5000_exp) + #if moon_sep > 160. - 1.1 * moon_alt: print('moon sep > 160 - 1.1 moon alt', Isky5000_exp) + #if moon_sep > moon_alt + 250: print('moon sep > 250. + moon alt', Isky5000_exp) + #if moon_sep < moon_alt - 50: print('moon sep > moon alt - 50', moon_sep, moon_alt, Isky5000_exp) + return Isky5000_exp/Isky5000_nom def Isky5000_notwilight_regression(airmass, moon_frac, moon_sep, moon_alt): @@ -351,42 +354,62 @@ def Isky5000_notwilight_regression(airmass, moon_frac, moon_sep, moon_alt): twilight in units of erg/s/cm^2/A/arcsec^2 ''' # polynomial of order - norder = 5 + norder = 6 # polynomial regression cofficients for estimating exposure time factor during # non-twilight from airmass, moon_frac, moon_sep, moon_alt - coeffs = np.array([ - 1.19728287e-02, 2.68146174e-02, 2.78840941e-02, 4.41723544e-03, - -1.02111676e-01, 7.65110151e-03, 1.63878220e-02, -3.08958375e-02, - 3.85657345e-01, 6.83474920e-03, 3.06045766e-01, 2.55798622e-01, - 1.19306254e-03, -1.04242412e-02, -3.34780689e-03, -1.68835833e-02, - 1.11378109e-02, -1.75611187e-01, -2.38751513e-01, 1.46593863e-02, - -4.36980545e-02, -1.78128533e-01, 5.85014404e-03, 1.37745149e-02, - -1.81392009e-03, -1.27115538e-02, 9.35953790e-02, -1.52422946e-01, - -1.74469624e-04, -6.01814594e-03, -2.13820870e-04, -7.40115962e-05, - -4.43960687e-05, 1.31544250e-04, 5.23119315e-05, -3.07378710e-02, - 3.50689477e-02, 1.85015093e-01, 6.13788213e-02, 1.98622649e-02, - -7.03566704e-02, -2.70010457e-01, -4.13256930e-03, -4.65235557e-03, - 2.08199644e-03, 4.06645589e-03, 2.08415057e-01, 4.87587825e-01, - -7.34996335e-03, 1.37774943e-03, 6.26446800e-03, 2.14264776e-05, - -1.34520777e-05, -1.28680576e-04, -2.98020923e-05, 1.41853618e-03, - -2.22111586e-01, -5.77624465e-01, -1.09989601e-02, -2.77969195e-04, - 2.17332386e-03, 1.47785979e-04, 1.63001051e-04, 3.10888811e-05, - -5.98993375e-05, -8.54879999e-08, 3.29709501e-07, 2.92447174e-08, - -2.79763418e-07, -2.38343236e-08, -1.10769335e-01, 8.24805309e-02, - -2.67759016e-02, -1.09365333e-03, 9.87385635e-03, 3.12620536e-02, - 9.15874867e-02, -2.52374361e-05, -4.92100655e-04, -4.21467759e-04, - 5.28695530e-03, 1.93471362e-01, -5.37647583e-02, -4.34387094e-04, - -2.25656833e-03, -7.15039182e-04, 1.18896366e-05, 2.77215519e-05, - 2.81214300e-05, -2.48811976e-06, 1.29700918e-02, -1.92447646e-02, - -2.44212957e-01, -5.99603226e-03, -2.76089664e-03, 9.14612070e-04, - 6.91296852e-05, 1.02360578e-04, 1.16371717e-05, -2.52971480e-05, - -1.37227312e-07, -3.75516901e-07, 5.48164999e-08, 5.06994020e-07, - 1.93984829e-07, 5.32931890e-02, 5.79701139e-01, 4.34770814e-01, - -6.79348983e-03, -1.31426633e-02, -6.14432600e-05, 1.39891173e-04, - 3.29198340e-04, 1.11812563e-04, -1.17664897e-05, -1.29847018e-06, - -3.91089212e-06, -2.85449951e-06, -3.27679136e-07, 3.48519396e-07, - 3.29795490e-09, 7.65512533e-09, 6.37102357e-09, 5.95568585e-10, - -2.63641980e-09, -1.30434882e-09]) + coeffs = np.array([ 1.37122205e-02, 1.19640701e-02, 8.16266472e-03, 1.39280691e-01, + -3.34998306e-05, 1.03669586e-02, 4.14594997e-03, 4.98021653e-02, + 2.04047595e-03, 2.65405851e-03, 9.64071111e-02, 5.05776641e-02, + -1.83615882e-02, -3.08046695e-02, -6.42666661e-03, 1.85826555e-02, + -1.31684802e-03, -6.31654630e-02, -2.75126470e-02, -6.66878948e-04, + 5.98373278e-02, -2.63741965e-02, 3.76787045e-02, 4.90090134e-02, + 1.88864338e-02, -4.62783691e-03, 2.14026758e-02, -2.84568954e-03, + 2.70850466e-02, 1.81878506e-02, 1.19064002e-03, 2.01685330e-05, + 6.32815029e-05, 2.19722530e-04, -2.84904174e-05, 4.19099793e-02, + -3.35908059e-03, -4.04610260e-02, 5.02943104e-02, -2.51984239e-03, + -1.07739033e-03, -1.30186514e-01, -4.01654696e-02, -3.46066336e-02, + -1.18782966e-02, -1.50282808e-03, 2.07729906e-02, -6.47500952e-02, + -3.04598385e-02, -1.57347975e-02, 2.54821102e-03, 3.64911542e-05, + -4.62796429e-05, -3.23300406e-04, -1.78909655e-04, -1.24732804e-03, + 6.34377564e-02, -2.72372301e-01, -2.71942506e-02, -1.10716987e-02, + 4.76993374e-03, -3.56730604e-04, -1.03304340e-04, -1.68009283e-04, + -4.62115541e-05, 2.50550445e-06, 9.68241078e-07, 8.13914082e-07, + 5.54375178e-07, 1.48456774e-06, 7.75842516e-02, 5.87414817e-04, + 7.64725150e-02, 2.20412705e-02, -3.13542091e-04, -2.85748098e-02, + -5.79352294e-02, 1.68630867e-02, 8.05659390e-03, -9.36412652e-04, + 6.02392576e-03, 6.04139180e-02, -3.02037097e-02, 2.86976403e-02, + -8.55838311e-04, 4.45485436e-03, -5.92335120e-05, 1.30608798e-04, + 2.76992145e-04, 1.73417298e-04, 8.46802477e-05, 1.44803613e-01, + -1.41019852e-01, -8.82311141e-03, 2.40245147e-02, 8.55080444e-03, + 1.65751285e-04, -3.06026118e-04, 2.47806624e-05, -1.04175553e-04, + -9.60076655e-07, 1.25270206e-06, -2.26510053e-07, -8.13034123e-07, + -3.62994514e-07, 8.10275226e-03, 1.30265610e-01, -2.59539811e-02, + 2.81777445e-02, -2.72609442e-02, -6.15755933e-05, -6.96052078e-05, + 4.44190852e-04, 1.88926726e-04, -9.41703486e-05, 3.33875647e-06, + 4.69869178e-06, 3.38184992e-07, 3.36524261e-07, 7.88612964e-07, + -2.53343718e-08, -6.01290739e-08, -3.55412440e-08, -9.68385387e-09, + -2.70149823e-09, -7.00167017e-09, 1.20582063e-01, 1.44232169e-02, + -5.99644863e-02, -5.46588272e-02, 3.79091286e-03, 1.13505888e-01, + 1.40498432e-01, -1.28698148e-03, 3.88926808e-03, 3.27254434e-03, + 1.84300473e-02, 7.48888214e-02, 3.77593511e-02, -1.68928339e-02, + -1.38366855e-02, -7.65439127e-03, 3.79474387e-05, -7.28387298e-05, + -1.86982036e-04, -8.33117928e-05, 2.20980133e-03, 2.16668147e-01, + 1.36526658e-01, 4.16746892e-03, -6.47206461e-03, -2.96650732e-03, + 8.65687911e-05, 4.72609545e-04, 4.72384634e-04, 1.51577611e-04, + -3.88374762e-07, -1.41057120e-06, -1.18013283e-08, 1.54156990e-06, + 4.42399919e-07, -4.41145646e-04, 1.39000704e-01, -4.37408535e-01, + -1.54311899e-02, -1.91755401e-03, 4.36138614e-03, -4.42645315e-05, + -8.34403032e-06, -1.53210208e-04, -3.48132905e-05, 7.41531513e-07, + -2.36866678e-06, -4.74525136e-06, -4.84825494e-06, -1.02298977e-06, + -1.75128388e-09, 1.07399591e-08, 1.26208068e-08, 4.31576944e-09, + -5.13099539e-09, -1.26290889e-09, 1.32565239e-02, 1.36316867e-01, + 4.99691249e-01, -1.74968493e-02, -3.34139883e-03, -2.42772081e-03, + 3.36178419e-04, 2.23896110e-04, 1.73634367e-04, -1.38246185e-06, + -2.40280959e-06, -2.66317367e-06, -2.80539906e-06, -4.19012237e-07, + 4.36408043e-07, -5.34778782e-09, -2.87653770e-08, -2.16845793e-09, + 3.21431780e-08, 2.09245593e-08, 1.66330738e-09, 1.08938813e-10, + 3.53672248e-10, 3.43915997e-10, 7.83100372e-11, -2.92671982e-11, + -5.14346513e-13, 1.39651631e-11]) theta = np.atleast_2d(np.array([airmass, moon_frac, moon_alt, moon_sep]).T) @@ -404,12 +427,8 @@ def Isky5000_twilight_regression(airmass, sun_sep, sun_alt): during bright time. Surface brightness is based on a polynomial regression model that was fit - to twilight contributions measured from BOSS exposures with sun altitutde > - -18deg. - - TODO: - * update fit using SV1 twilight exposures. - + to twilight contributions measured from SV1 and BOSS exposures with sun + altitutde > -18deg. For details, see jupyter notebook: https://github.com/desi-bgs/bgs-cmxsv/blob/55850e44d65570da69de0788c652cff698416834/doc/nb/sky_model_twilight_fit.ipynb @@ -429,10 +448,13 @@ def Isky5000_twilight_regression(airmass, sun_sep, sun_alt): Array of twilight contribution to sky surface brightness at 5000A for bright time in units of erg/s/cm^2/A/arcsec^2 ''' - norder = 1 + norder = 2 # coefficients fit in notebook # https://github.com/desi-bgs/bgs-cmxsv/blob/521211dccff7cb3b71b07d84522802acade26071/doc/nb/sky_model_twilight_fit.ipynb - skymodel_coeff = np.array([3.50300842, 2.63453598, 0.40491282, 0.0096106 ]) + skymodel_coeff = np.array([ + 7.74536884e-01, 1.25866845e+00, 1.06073841e+00, 2.27595902e-01, + 9.12671470e-01, -4.16253157e-02, -2.26899688e-02, 3.52752773e-02, + 6.39485772e-03, -5.03964852e-04]) theta = np.atleast_2d(np.array([airmass, sun_alt, sun_sep]).T) @@ -686,7 +708,8 @@ def start(self, mjd_now, tileid, program, snr2frac, exposure_factor, seeing, tra self.snr2frac_target = snr2frac + (texp_remaining / nexp) / self.texp_total # Initialize signal and background rate factors. self.srate0 = np.sqrt(self.weather_factor(seeing, transp, 1.0)) - self.brate0 = sky + # effective sky / nominal sky based on SurveySpeed calculation + self.brate0 = sky + 0.932/3.373 * (self.texp_total / 0.0115741) + 1.71 / 3.373 self.signal = 0. self.background = 0. self.last_snr2frac = 0. @@ -724,9 +747,9 @@ def update(self, mjd_now, seeing, transp, sky): dt = mjd_now - self.mjd_last self.mjd_last = mjd_now srate = np.sqrt(self.weather_factor(seeing, transp, 1.0)) - brate = sky + # effective sky / nominal sky based on SurveySpeed calculation + brate = sky + 0.932/3.373 * (dt / 0.0115741) + 1.71 / 3.373 self.signal += dt * srate / self.srate0 - #self.background += dt * (srate + brate) / (self.srate0 + self.brate0) self.background += dt * brate / self.brate0 self._snr2frac = self._snr2frac_start + self.signal ** 2 / self.background / self.texp_total if self.save_history: