Default iono delay (#138)

* default iono delay

* no self in slant delay

* helper function

* need this

laika/astro_dog.py

@@ -9,7 +9,7 @@
 from .downloader import download_orbits_gps, download_orbits_russia_src, download_nav, download_ionex, download_dcb, download_prediction_orbits_russia_src
 from .downloader import download_cors_station
 from .trop import saast
-from .iono import IonexMap, parse_ionex
+from .iono import IonexMap, parse_ionex, get_slant_delay
 from .dcb import DCB, parse_dcbs
 from .gps_time import GPSTime
 from .dgps import get_closest_station_names, parse_dgps
@@ -345,7 +345,11 @@ def get_delay(self, prn, time, rcv_pos, no_dgps=False, signal='C1C', freq=None):
     # When using internet we expect all data or return None
     if self.auto_update and (ionex is None or dcb is None or freq is None):
       return None
-    iono_delay = ionex.get_delay(rcv_pos, az, el, sat_pos, time, freq) if ionex is not None else 0.
+    if ionex is not None:
+      iono_delay = ionex.get_delay(rcv_pos, az, el, sat_pos, time, freq)
+    else:
+      # 5m vertical delay is a good default
+      iono_delay = get_slant_delay(rcv_pos, az, el, sat_pos, time, freq, vertical_delay=5.0)
     trop_delay = saast(rcv_pos, el)
     code_bias = dcb.get_delay(signal) if dcb is not None else 0.
     return iono_delay + trop_delay + code_bias

laika/iono.py

@@ -9,6 +9,17 @@
 # Altitude of Ionospheric-pierce-point
 IPP_ALT = 6821000
 
+def get_alpha_beta(rcv_pos, el):
+    geocentric_alt = np.linalg.norm(rcv_pos)
+    alpha = np.pi/2 + el
+    beta = np.arcsin(geocentric_alt*np.sin(alpha)/IPP_ALT)
+    return alpha, beta
+
+def get_slant_delay(rcv_pos, az, el, sat_pos, time, freq, vertical_delay):
+    alpha, beta = get_alpha_beta(rcv_pos, el)
+    slant_delay = vertical_delay * ((1 - ((EARTH_RADIUS * np.sin(beta)) /
+                                          (EARTH_RADIUS + 3.5e5))**2)**(-0.5))
+    return slant_delay
 
 def closest_in_list(lst, val, num=2):
   """
@@ -134,18 +145,16 @@ def get_delay(self, rcv_pos, az, el, sat_pos, time, freq):
     # To get a delay from a TEC map, we need to calculate
     # the ionospheric pierce point, geometry described here
     # https://en.wikipedia.org/wiki/Ionospheric_pierce_point
+    alpha, beta = get_alpha_beta(rcv_pos, el)
     conv = LocalCoord.from_ecef(rcv_pos)
-    geocentric_alt = np.linalg.norm(rcv_pos)
-    alpha = np.pi/2 + el
-    beta = np.arcsin(geocentric_alt*np.sin(alpha)/IPP_ALT)
     gamma = np.pi - alpha - beta
+    geocentric_alt = np.linalg.norm(rcv_pos)
     ipp_dist = geocentric_alt*np.sin(gamma)/np.sin(beta)
     ipp_ned = conv.ecef2ned(sat_pos)*(ipp_dist)/np.linalg.norm(sat_pos)
     ipp_geo = conv.ned2geodetic(ipp_ned)
     factor = 40.30E16 / (freq**2) * 10**(self.exp)
     vertical_delay = self.get_TEC(ipp_geo, time) * factor
-    slant_delay = vertical_delay * ((1 - ((EARTH_RADIUS * np.sin(beta)) /
-                                          (EARTH_RADIUS + 3.5e5))**2)**(-0.5))
+    slant_delay = get_slant_delay(rcv_pos, az, el, sat_pos, time, freq, vertical_delay)
     return slant_delay
 
   @staticmethod