Short sample processing fixes

peregrine/acquisition.py

@@ -54,14 +54,16 @@ class Acquisition:
     Array of samples to use for acquisition. Can be `None` but in this case
     `init_samples` *must* be called with an array of samples before any other
     acquisition functions are used.
-  sampling_freq : float, optional
+  sampling_freq : float
     The sampling frequency of the samples in Hz.
-  IF : float, optional
+  IF : float
     The receiver intermediate frequency used when capturing the samples.
-  samples_per_code : float, optional
+  samples_per_code : float
     The number of samples corresponding to one code length.
-  code_length : int, optional
+  code_length : int
     The number of chips in the chipping code.
+  n_codes_integrate : int, optional
+    The number of codes to integrate.
   offsets : int, optional
     Offsets, in units of code length (1ms), to use when performing long
     integrations to avoid clobbering by nav bit edges.

peregrine/defaults.py

@@ -15,6 +15,7 @@
 useCache = True
 cacheDir = 'cache'
 ephemMaxAge = 4 * 3600.0  # Reject an ephemeris entry if older than this
+elevMask = None # Optional elevation mask (degrees), None to disable
 
 # the size of the sample data block processed at a time
 processing_block_size = int(20e6)  # [samples]

peregrine/ephemeris.py

@@ -98,6 +98,27 @@ def read4(f):
         filename, count_healthy)
     return ephem
 
+def rinex_creation_date(filename):
+    """ Read RINEX file creation time
+
+    Returns:
+      datetime object
+    """
+    dt = None
+    with open(filename,'r') as f:
+        while True:
+            line = f.readline()[:-1]
+            if not line:
+                break
+            if "PGM / RUN BY / DATE" not in line:
+                continue
+
+            # Extract creation time
+            timestring = '-'.join(line.split()[2:4])
+            dt = datetime.strptime(timestring,"%Y%m%d-%H%M%S")
+
+    return dt
+
 def obtain_ephemeris(t, settings):
     """
     Finds an appropriate GNSS ephemeris file for a certain time,
@@ -121,17 +142,31 @@ def obtain_ephemeris(t, settings):
     """
     print "Obtaining ephemeris file for ", t
 
-    #TODO: If it's today and more than 1 hr old, check for an update
-
     filename = t.strftime("brdc%j0.%yn")
     filedir = os.path.join(settings.cacheDir, "ephem")
     filepath = os.path.join(filedir, filename)
     url = t.strftime(
         'http://qz-vision.jaxa.jp/USE/archives/ephemeris/%Y/') + filename
+
+    # If not in cache, then fetch it
     if not os.path.isfile(filepath):
         if not os.path.exists(filedir):
             os.makedirs(filedir)
         _, hdrs = urllib.urlretrieve(url, filepath)
+
+    # Otherwise, use cache file if isn't stale
+    else:
+        rinex_gen_date = rinex_creation_date(filepath)
+        if rinex_gen_date is None:
+            # Something wrong with RINEX format, update cache
+            _, hdrs = urllib.urlretrieve(url, filepath)
+        else:
+            # If rinex generated more than an hour then fetch it again
+            if (t-rinex_gen_date) > timedelta(hours=1):
+                _, hdrs = urllib.urlretrieve(url, filepath)
+
+    rinex_gen_date = rinex_creation_date(filepath)
+
     ephem = load_rinex_nav_msg(filepath, t, settings)
     if len(ephem) < 22:
         raise ValueError(

peregrine/gps_time.py

@@ -1,5 +1,78 @@
 import datetime
 
+# NTP timestamps are in units of seconds since the NTP epoch
+# See leap-seconds.list for further details.
+NTP_EPOCH = datetime.datetime(1900,1,1,0,0,0)
+
+def sec_since_NTP_epoch(dt):
+    """Return the number of seconds since the NTP epoch."""
+    return (dt-NTP_EPOCH).total_seconds()
+
+def load_leapsecond_table(f="/usr/share/zoneinfo/leap-seconds.list"):
+    """"
+    Loads leap second table from system.
+
+    The default file location is appropriate for Ubuntu and is provided by
+    the tzdata package.  Refer to the leap-seconds.list file for formatting
+    information.
+
+    Parameters
+    ----------
+    f : string
+      Path to a leap-seconds.list file
+
+    Returns: List of tuples in chronological order each containing an
+    epoch start time and number of leap seconds applicable for that epoch.
+    """
+
+    # Check the expiration date in the file
+    with open(f,'r') as fp:
+        for line in fp:
+            if line[0:2] != "#@":
+                continue
+
+            expiration_time = int(line.split('\t')[-1])
+            now = sec_since_NTP_epoch(datetime.datetime.now())
+
+            if expiration_time<now:
+                raise UserWarning("leap-seconds.list file has expired, update tzdata")
+
+    # Load the load the table
+    with open(f,'r') as fp:
+        table = []
+        for line in fp:
+            if line[0]=="#":
+                continue
+            raw = line.split('\t')
+            table.append( (int(raw[0]), int(raw[1])) )
+
+    # Add the expiration time to the end of the table
+    table.append( (expiration_time, None) )
+
+    return table
+
+def get_gpst_leap_seconds(dt):
+    """Returns the number of leap seconds at the time provided."""
+
+    # Load leap second data from system tzdata file
+    lstable = load_leapsecond_table()
+
+    t = sec_since_NTP_epoch(dt)
+
+    ls = None
+    for epoch_start,leapseconds in lstable:
+        if t>=epoch_start:
+            ls = leapseconds
+
+    # Raise warning if specified time is after expiry date of leap second table
+    if ls==None:
+        raise UserWarning("Specified datetime is after expiry time of leap second table.")
+
+    # GPS leap seconds relative to a Jan 1 1980, where TAI-UTC was 19 seconds.
+    gps_leap_seconds = ls - 19
+
+    return datetime.timedelta(seconds = gps_leap_seconds)
+
 def datetime_to_tow(t, mod1024 = True):
     """
     Convert a Python datetime object to GPS Week and Time Of Week.
@@ -39,12 +112,7 @@ def gpst_to_utc(t_gpst):
     -------
     t_utc : datetime
     """
-    t_utc = t_gpst - datetime.timedelta(seconds = 16)
-    if t_utc <= datetime.datetime(2012, 7, 1) or \
-       t_utc >= datetime.datetime(2015, 7, 1):
-        raise ValueError("Don't know how many leap seconds to use.  " +
-                         "Please implement something better in gpst_to_utc()")
-    return t_utc
+    return t_gpst - get_gpst_leap_seconds(t_gpst)
 
 def utc_to_gpst(t_utc):
     """
@@ -59,9 +127,7 @@ def utc_to_gpst(t_utc):
     -------
     t_gpst : datetime
     """
-    t_gpst = t_utc + datetime.timedelta(seconds = 16)
-    if t_utc <= datetime.datetime(2012, 7, 1) or \
-       t_utc >= datetime.datetime(2015, 7, 1):
-        raise ValueError("Don't know how many leap seconds to use.  " +
-                         "Please implement something better in utc_to_gpst()")
-    return t_gpst
+    # TODO: there may be an unhandled corner case here for conversions
+    # where the GPST-to-UTCT interval spans a leap second.
+
+    return t_utc + get_gpst_leap_seconds(t_utc)