Merge pull request #21 from desihub/refactorI

Re-added avoidobject, with the option to use the JPL dependent versio…

README.md

@@ -4,16 +4,16 @@ The full path to surveysim/py should be added to PYTHONPATH.
 
 To run, for example:
 
-	>>> from surveysim.wrapper import surveySim
-	>>> surveySim((2016, 12, 27), (2017, 1, 4), 123456)
+	>>> from surveysim.surveysim import surveySim
+	>>> surveySim((2016, 12, 27), (2017, 1, 4), seed=123456, use_jpl=False)
 
 The optional seed for the weather module's random number generator has to
-be an int or array_like convertible to an unsigned 32-bit integer.
+be an int or array_like convertible to an unsigned 32-bit integer; use_jpl specifies which version of avoidobject.py to use.
 
 To run the plotting tool:
 
 	>>> from surveysim.plotsurvey import plotsurvey
-	>>> plotsurvey("your_output_file.fits", plot_type='t')
+	>>> plotsurvey("obslist{_all|YYYYMMDD}.fits", plot_type='t', program='m')
 
-The default filename is ./obslist_all.fits and plot_type is either 'f' (footprint, default), 'h' (histograms), 't' (time evolution) or 'e' (exposure time).
+The default filename is ./obslist_all.fits; plot_type is either 'f' (footprint, default), 'h' (histograms), 't' (time evolution) or 'e' (exposure time); and program us either 'm' (main survey), 'b' (BGS) or 'a' (all).
 

doc/changes.rst

@@ -2,12 +2,13 @@
 surveysim change log
 ====================
 
-1.1.2 (unreleased)
+0.2.0 (2016-11-18)
 ------------------
 
-* No changes yet
+* Modified some file names
+* Moved some functions from one file to another
 
-1.1.1 (2016-11-14)
+0.1.1 (2016-11-14)
 ------------------
 
 * fixed crash at end and data/ install (PR #3)

py/surveysim/_version.py

@@ -1,2 +1,2 @@
-__version__ = '0.1.1'
+__version__ = '0.2.0'
 

py/surveysim/avoidobject.py

@@ -0,0 +1,98 @@
+import ephem
+from datetime import datetime
+import numpy as np
+from surveysim.kpno import mayall
+
+MIN_VENUS_SEP = np.radians(2.0)
+MIN_MARS_SEP = np.radians(2.0)
+MIN_JUPITER_SEP = np.radians(2.0)
+MIN_SATURN_SEP = np.radians(2.0)
+MIN_NEPTUNE_SEP = np.radians(2.0)
+MIN_URANUS_SEP = np.radians(2.0)
+MIN_CERES_SEP = np.radians(2.0)
+
+def avoidObject(datetime, ra0, dec0):
+    """
+    Checks whether all the objects on the list are far enough away from
+    the input coordinates.
+    The current list has: Venus, Mars, Jupiter, Saturn, Neptune, Uranus;
+    the Moon is treated separately.
+
+    Args:
+        datetime: datetime object; should have timezone info
+        ra0: float (apparent or observed, degrees)
+        dec0: float (apparent or observed, degrees)
+
+    Returns:
+        bool, True if all objects on the list are far enough away
+    """
+
+    ra = np.radians(ra0)
+    dec = np.radians(dec0)
+
+    dt = ephem.Date(datetime.datetime)
+    gatech = ephem.Observer()
+    gatech.lon, gatech.lat = np.radians(mayall.west_lon_deg), np.radians(mayall.lat_deg)
+    gatech.date = dt
+    gatech.epoch = dt
+
+    venus = ephem.Venus()
+    venus.compute(gatech)
+    if ephem.separation(venus, (ra, dec)) < MIN_VENUS_SEP:
+        return False
+    mars = ephem.Mars()
+    mars.compute(gatech)
+    if ephem.separation(mars, (ra, dec)) < MIN_MARS_SEP:
+        return False
+    #ceres = ephem.Ceres()
+    #ceres.compute(gatech)
+    #if ephem.separation(ceres, (ra, dec)) < MIN_CERES_SEP:
+    #    return False
+    jupiter = ephem.Jupiter()
+    jupiter.compute(gatech)
+    if ephem.separation(jupiter, (ra, dec)) < MIN_JUPITER_SEP:
+        return False
+    saturn = ephem.Saturn()
+    saturn.compute(gatech)
+    if ephem.separation(saturn, (ra, dec)) < MIN_SATURN_SEP:
+        return False
+    neptune = ephem.Neptune()
+    neptune.compute(gatech)
+    if ephem.separation(neptune, (ra, dec)) < MIN_NEPTUNE_SEP:
+        return False
+    uranus = ephem.Uranus()
+    uranus.compute(gatech)
+    if ephem.separation(uranus, (ra, dec)) < MIN_URANUS_SEP:
+        return False
+
+    # If still here, return True
+    return True
+
+def moonLoc (datetime, ra0, dec0):
+    """
+    Returns the distance to the Moon if RA and DEC as well as alt, az.
+
+    Args:
+        datetime: datetime object; should have timezone info
+        ra0: float (apparent or observed, degrees)
+        dec0: float (apparent or observed, degrees)
+
+    Returns:
+        float, distance from the Moon (degrees)
+        float, Moon altitude (degrees)
+        float, Moon azimuth (degrees) 
+    """
+
+    dt = ephem.Date(datetime.datetime)
+    gatech = ephem.Observer()
+    gatech.lon, gatech.lat = np.radians(mayall.west_lon_deg), np.radians(mayall.lat_deg)
+    gatech.date = dt
+    gatech.epoch = dt
+
+    moon = ephem.Moon()
+    moon.compute(gatech)
+    ra = np.radians(ra0)
+    dec = np.radians(dec0)
+    moondist = ephem.separation(moon, (ra, dec))
+
+    return np.degrees(moondist), np.degrees((moon.alt)), np.degrees((moon.az))

py/surveysim/nextobservation.py

@@ -1,17 +1,15 @@
 import numpy as np
 import astropy.io.fits as pyfits
-#from surveysim.utils import angsep
 from surveysim.exposurecalc import airMassCalculator
 from surveysim.utils import mjd2lst
-from surveysim.observefield import setup_time
 from astropy.time import Time
 from desitarget.targetmask import obsconditions as obsbits
-from desisurvey.avoidobject import avoidObject, moonLoc
+
 MAX_AIRMASS = 10.0 #3.0 This new bound effectively does nothing.
 MIN_MOON_SEP = 90.0
 MIN_MOON_SEP_BGS = 5.0
 
-def nextFieldSelector(obsplan, mjd, conditions, tilesObserved, slew, previous_ra, previous_dec):
+def nextFieldSelector(obsplan, mjd, conditions, tilesObserved, slew, previous_ra, previous_dec, use_jpl=False):
     """
     Returns the first tile for which the current time falls inside
     its assigned LST window and is far enough from the Moon and
@@ -25,6 +23,7 @@ def nextFieldSelector(obsplan, mjd, conditions, tilesObserved, slew, previous_ra
         slew: bool, True if a slew time needs to be taken into account
         previous_ra: float, ra of the previous observed tile (degrees)
         previous_dec: float, dec of the previous observed tile (degrees)
+        use_jpl: bool, True if using jplephem and astropy instead of pyephem
 
     Returns:
         target: dictionnary containing the following keys:
@@ -33,6 +32,10 @@ def nextFieldSelector(obsplan, mjd, conditions, tilesObserved, slew, previous_ra
                 'Exposure', 'obsSN2', 'obsConds'
         overhead: float (seconds)
     """
+    if (use_jpl):
+        from desisurvey.avoidobjectJPL import avoidObject, moonLoc
+    else:
+        from surveysim.avoidobject import avoidObject, moonLoc
 
     hdulist = pyfits.open(obsplan)
     tiledata = hdulist[1].data
@@ -85,3 +88,27 @@ def nextFieldSelector(obsplan, mjd, conditions, tilesObserved, slew, previous_ra
     else:
         target = None
     return target, overhead
+
+def setup_time(slew, dra, ddec):
+    """
+    Computes setup time: slew and focus (assumes readout can proceed during
+    slew.
+
+    Args:
+        slew: bool, True if slew time needs to be taken into account
+        dra: float, difference in RA between previous and current tile (degrees)
+        ddec: float, difference in DEC between previous and current tile (degrees)
+
+    Returns:
+        float, total setup time (seconds)
+    """
+
+    focus_time = 30.0
+    slew_time = 0.0
+    if slew:
+        d = np.maximum(dra, ddec)
+        slew_time = 11.5 + d/0.45
+    overhead = focus_time + slew_time
+    if overhead < 120.0:
+        overhead = 120.0
+    return overhead

py/surveysim/wrapper.py → py/surveysim/nightops.py

@@ -3,11 +3,11 @@
 import numpy as np
 from astropy.time import Time
 from datetime import datetime, timedelta
-from surveysim.weather import weatherModule
+#from surveysim.weather import weatherModule
 from surveysim.nightcal import getCal
 from surveysim.exposurecalc import expTimeEstimator, airMassCalculator
 import astropy.io.fits as pyfits
-from surveysim.afternoonplan import surveyPlan
+#from surveysim.afternoonplan import surveyPlan
 from surveysim.nextobservation import nextFieldSelector
 from surveysim.observefield import observeField
 from astropy.table import Table, vstack
@@ -21,11 +21,11 @@ class obsCount:
     will have its own file with its number as part of the filename.
     """
 
-    def __init__(self):
+    def __init__(self, start_val=0):
         """
         Initialise the counter to zero
         """
-        self.obsNumber = 0
+        self.obsNumber = start_val
 
     def update(self):
         """
@@ -54,7 +54,7 @@ def update(self):
             partFileName = '0000000' + str(obsNumber)
         return partFileName
 
-def nightOps(day_stats, obsplan, w, ocnt, tilesObserved, tableOutput=True):
+def nightOps(day_stats, obsplan, w, ocnt, tilesObserved, tableOutput=True, use_jpl=False):
     """
     Carries out observations during one night and writes the output to disk
 
@@ -98,7 +98,7 @@ def nightOps(day_stats, obsplan, w, ocnt, tilesObserved, tableOutput=True):
             conditions = w.updateValues(conditions, mjd)
 
             lst = mjd2lst(mjd)
-            target, setup_time = nextFieldSelector(obsplan, mjd, conditions, tilesObserved, slew, ra_prev, dec_prev)
+            target, setup_time = nextFieldSelector(obsplan, mjd, conditions, tilesObserved, slew, ra_prev, dec_prev, use_jpl)
             if target != None:
                 # Compute mean to apparent to observed ra and dec???
                 airmass = airMassCalculator(target['RA'], target['DEC'], lst)
@@ -197,60 +197,3 @@ def nightOps(day_stats, obsplan, w, ocnt, tilesObserved, tableOutput=True):
 
     return tilesObserved
 
-def surveySim(sd0, ed0, seed=None, tilesubset=None):
-    """
-    Main driver for survey simulations.
-
-    Args:
-        sd0: tuple of three integers: startyear, startmonth, startday
-        ed0: tuple of three integers: endyear, endmonth, endday
-
-    Optional:
-        seed: integer, to initialise random number generator for weather simulator
-        tilesubset : array of integer tileIDs to use while ignoring others
-            in the DESI footprint
-    """
-
-    # Note 1900 UTC is midday at KPNO
-    (startyear, startmonth, startday) = sd0
-    startdate = datetime(startyear, startmonth, startday, 19, 0, 0)
-    (endyear, endmonth, endday) = ed0
-    enddate = datetime(endyear, endmonth, endday, 19, 0, 0)
-
-    sp = surveyPlan(tilesubset=tilesubset)
-    day0 = Time(datetime(startyear, startmonth, startday, 19, 0, 0))
-    mjd_start = day0.mjd
-    w = weatherModule(startdate, seed)
-    ocnt = obsCount()
-
-    tile_file = 'tiles_observed.fits'
-    if os.path.exists(tile_file):
-        tilesObserved = Table.read(tile_file, format='fits')
-    else:
-        print("The survey will start from scratch.")
-        tilesObserved = Table(names=('TILEID', 'STATUS'), dtype=('i8', 'i4'))
-        tilesObserved.meta['MJDBEGIN'] = mjd_start
-
-    oneday = timedelta(days=1)
-    day = startdate
-    day_monsoon_start = 13
-    month_monsoon_start = 7
-    day_monsoon_end = 27
-    month_monsoon_end = 8
-    survey_done = False
-    while (day <= enddate and survey_done == False):
-        if ( not (day >= datetime(day.year, month_monsoon_start, day_monsoon_start) and
-                  day <= datetime(day.year, month_monsoon_end, day_monsoon_end)) ):
-            day_stats = getCal(day)
-            ntodate = len(tilesObserved)
-            w.resetDome(day)
-            tiles_todo, obsplan = sp.afternoonPlan(day_stats, tilesObserved)
-            tilesObserved = nightOps(day_stats, obsplan, w, ocnt, tilesObserved)
-            t = Time(day, format = 'datetime')
-            ntiles_tonight = len(tilesObserved)-ntodate
-            print ('On the night starting ', t.iso, ', we observed ', ntiles_tonight, ' tiles.')
-            if (tiles_todo-ntiles_tonight) == 0:
-                survey_done = True
-        day += oneday
-
-    tilesObserved.write(tile_file, format='fits', overwrite=True)

py/surveysim/observefield.py

@@ -24,27 +24,3 @@ def observeField(target, exposure):
 
     return status, real_exposure, realSN2
 
-def setup_time(slew, dra, ddec):
-    """
-    Computes setup time: slew and focus (assumes readout can proceed during
-    slew.
-
-    Args:
-        slew: bool, True if slew time needs to be taken into account
-        dra: float, difference in RA between previous and current tile (degrees)
-        ddec: float, difference in DEC between previous and current tile (degrees)
-
-    Returns:
-        float, total setup time (seconds)
-    """
-
-    focus_time = 30.0
-    slew_time = 0.0
-    if slew:
-        d = np.maximum(dra, ddec)
-        slew_time = 11.5 + d/0.45
-    overhead = focus_time + slew_time
-    if overhead < 120.0:
-        overhead = 120.0
-    return overhead
-