Repository: hirokawa/cssrlib (language: python)

Issue:
Revised modules using new RINEXv3 signal class



Discussion / hints:

Fix patch:
diff --git a/README.md b/README.md
index cdd299c..2ee28bb 100644
--- a/README.md
+++ b/README.md
@@ -2,7 +2,7 @@
 CSSRlib - Toolkit for PPP-RTK/RTK in Python using Compact SSR
 *************************
 
-# What is CSSRlib? 
+# What is CSSRlib?
 
 CSSRLIB is a open toolkit in Python for high accuracy GNSS positioning. It supports SSR (State-Space Representation) based potitioning such as PPP (Precise Point Positioning) or PPP-RTK (Realtime Kinematic), but also supporting RTK. The goal of the CSSRlib toolkit is to provide an easy-to-understand open implementation to learn PPP/PPP-RTK positioning provided by satellite-based open PPP/PPP-RTK services such as QZSS CLAS, Galileo HAS, BeiDou 3 PPP. It also supports ground based open service by IGS. The code is based on RTKlib.
 
@@ -17,10 +17,9 @@ Prerequisites
 Additional python packages are required as prerequisites and can be installed via the following commands
 
 ```
-pip install bitstruct
-pip install cbitstruct
-pip install galois
-pip install cartopy
+pip install bitstruct, cbitstruct
+pip install galois, cartopy
+pip install notebook, numpy, matplotlib
 ```
 
 If the installation of `cartopy` fails, try installing `libgeos++-dev` first.
@@ -31,19 +30,21 @@ sudo apt-get install libgeos++-dev
 
 Install
 =======
-You can install CSSRlib using pip.
+
+You can install the official version of CSSRlib using pip
 
 ```
 pip install cssrlib
 ```
 
-additional python packages are required as prerequisites and can be installed via the following commands:
+If you want to install the development version from this repository, first clone or download the sources and then run
 
 ```
-pip install notebook, numpy, matplotlib
-pip install bitstruct,cbitstruct,galois,cartopy
+pip install .
 ```
 
+in the root directory, where the ``setup.cfg`` file is located.
+
 Testing
 =======
 
@@ -59,7 +60,7 @@ Run RTK sample.
  python test_rtk.py
 ```
 
-Other samples are also available in a separate repository `cssrlib-data` including :
+Other samples are also available in a separate repository [`cssrlib-data`](https://github.com/hirokawa/cssrlib-data) including :
 
 - Galileo-HAS decoder
 - BDS-PPP decoder
diff --git a/src/cssrlib/ephemeris.py b/src/cssrlib/ephemeris.py
index e68c3f0..c7c9765 100644
--- a/src/cssrlib/ephemeris.py
+++ b/src/cssrlib/ephemeris.py
@@ -132,7 +132,7 @@ def satposs(obs, nav, cs=None):
     for i in range(n):
         sat = obs.sat[i]
         sys, _ = sat2prn(sat)
-        if sys not in nav.gnss_t:
+        if sys not in obs.sig.keys():
             continue
         pr = obs.P[i, 0]
         t = timeadd(obs.t, -pr/rCST.CLIGHT)
diff --git a/src/cssrlib/gnss.py b/src/cssrlib/gnss.py
index 1a53e09..493449b 100644
--- a/src/cssrlib/gnss.py
+++ b/src/cssrlib/gnss.py
@@ -67,17 +67,18 @@ class rCST():
 
 class uGNSS(IntEnum):
     """ class for GNS systems """
+
     NONE = -1
 
     GPS = 0
-    SBS = 1
-    GAL = 2
+    GAL = 1
+    QZS = 2
     BDS = 3
-    QZS = 5
-    GLO = 6
-    IRN = 7
+    GLO = 4
+    SBS = 5
+    IRN = 6
 
-    GNSSMAX = 8
+    GNSSMAX = 7
 
     GPSMAX = 32
     GALMAX = 36
@@ -87,7 +88,15 @@ class uGNSS(IntEnum):
     SBSMAX = 24
     IRNMAX = 10
 
-    MAXSAT = GPSMAX+GLOMAX+GALMAX+BDSMAX+QZSMAX+SBSMAX+IRNMAX
+    GPSMIN = 0
+    GALMIN = GPSMIN+GPSMAX
+    QZSMIN = GALMIN+GALMAX
+    BDSMIN = QZSMIN+QZSMAX
+    GLOMIN = BDSMIN+BDSMAX
+    SBSMIN = GLOMIN+GLOMAX
+    IRNMIN = SBSMIN+SBSMAX
+
+    MAXSAT = GPSMAX+GALMAX+QZSMAX+BDSMAX+GLOMAX+SBSMAX+IRNMAX
 
 
 class uTYP(IntEnum):
@@ -100,8 +109,10 @@ class uTYP(IntEnum):
     D = 3
     S = 4
 
+
 class uSIG(IntEnum):
     """ class for signal band and attribute """
+
     NONE = -1
 
     # GPS   L1  1575.42 MHz
@@ -502,14 +513,13 @@ class Obs():
     """ class to define the observation """
 
     def __init__(self):
-        self.nm = 0
         self.t = gtime_t()
         self.P = []
         self.L = []
         self.S = []
         self.lli = []
-        self.data = []
         self.sat = []
+        self.sig = {}
 
 
 class Eph():
@@ -561,24 +571,23 @@ def __init__(self):
             [0.1167E+06, -0.2294E+06, -0.1311E+06, 0.1049E+07]])
         self.elmin = np.deg2rad(15.0)
         self.tidecorr = False
-        self.nf = 2
+        ######## START OBSOLETE ################################################
+        self.nf = 2  # TODO: make this obsolte if possible
+        ######## END   OBSOLETE ################################################
         self.ne = 0
         self.nc = 0
-        self.excl_sat = []
-        self.freq = [1.57542e9, 1.22760e9,  # L1,L2
-                     1.17645e9, 1.20714e9]  # L5/E5a/B2a,E5b/B2b
+        self.excl_sat = []  # Excluded satellites
         self.rb = [0, 0, 0]  # base station position in ECEF [m]
         self.smode = 0  # position mode 0:NONE,1:std,2:DGPS,4:fix,5:float
-        self.gnss_t = [uGNSS.GPS, uGNSS.GAL, uGNSS.QZS]
-        self.loglevel = 1
+        self.pmode = 1  # 0: static, 1: kinematic
+
+        self.monlevel = 1
         self.cnr_min = 35
         self.maxout = 5  # maximum outage [epoch]
 
         self.sat_ant = None
         self.rcv_ant = None
         self.rcv_ant_b = None
-        self.sig_tab = None
-        self.sig_tab_b = None
 
         # SSR correction placeholder
         self.dorb = np.zeros(uGNSS.MAXSAT)
@@ -588,9 +597,13 @@ def __init__(self):
 
         # satellite observation status
         self.fix = np.zeros((uGNSS.MAXSAT, self.nf), dtype=int)
+        # Measurement outage indicator
         self.outc = np.zeros((uGNSS.MAXSAT, self.nf), dtype=int)
+        # Carrier-phase processed indicator
         self.vsat = np.zeros((uGNSS.MAXSAT, self.nf), dtype=int)
+        # ??? set but not used
         self.lock = np.zeros((uGNSS.MAXSAT, self.nf), dtype=int)
+        # ??? not used
         self.slip = np.zeros((uGNSS.MAXSAT, self.nf), dtype=int)
 
         self.tt = 0
@@ -711,23 +724,23 @@ def time2str(t):
     return "{:04d}-{:02d}-{:02d} {:02d}:{:02d}:{:02d}"\
         .format(e[0], e[1], e[2], e[3], e[4], int(e[5]))
 
+
 def prn2sat(sys, prn):
     """ convert sys+prn to sat """
     if sys == uGNSS.GPS:
         sat = prn
-    elif sys == uGNSS.GLO:
-        sat = prn+uGNSS.GPSMAX
     elif sys == uGNSS.GAL:
-        sat = prn+uGNSS.GPSMAX+uGNSS.GLOMAX
-    elif sys == uGNSS.BDS:
-        sat = prn+uGNSS.GPSMAX+uGNSS.GLOMAX+uGNSS.GALMAX
+        sat = prn+uGNSS.GALMIN
     elif sys == uGNSS.QZS:
-        sat = prn-192+uGNSS.GPSMAX+uGNSS.GLOMAX+uGNSS.GALMAX+uGNSS.BDSMAX
+        sat = prn-192+uGNSS.QZSMIN
+    elif sys == uGNSS.GLO:
+        sat = prn+uGNSS.GLOMIN
+    elif sys == uGNSS.BDS:
+        sat = prn+uGNSS.BDSMIN
     elif sys == uGNSS.SBS:
-        sat = prn-100+uGNSS.GPSMAX+uGNSS.GLOMAX+uGNSS.GALMAX+uGNSS.BDSMAX+uGNSS.QZSMAX
+        sat = prn-100+uGNSS.SBSMIN
     elif sys == uGNSS.IRN:
-        sat = prn+uGNSS.GPSMAX+uGNSS.GLOMAX+uGNSS.GALMAX + \
-            uGNSS.BDSMAX+uGNSS.QZSMAX+uGNSS.SBSMAX
+        sat = prn+uGNSS.IRNMIN
     else:
         sat = 0
     return sat
@@ -735,26 +748,27 @@ def prn2sat(sys, prn):
 
 def sat2prn(sat):
     """ convert sat to sys+prn """
-    if sat > uGNSS.GPSMAX+uGNSS.GLOMAX+uGNSS.GALMAX+uGNSS.BDSMAX+uGNSS.QZSMAX+uGNSS.SBSMAX:
-        prn = sat-(uGNSS.GPSMAX+uGNSS.GLOMAX+uGNSS.GALMAX +
-                   uGNSS.BDSMAX+uGNSS.QZSMAX+uGNSS.SBSMAX)
+    if sat > uGNSS.MAXSAT:
+        prn = 0
+        sys = uGNSS.NONE
+    elif sat > uGNSS.IRNMIN:
+        prn = sat-uGNSS.IRNMIN
         sys = uGNSS.IRN
-    elif sat > uGNSS.GPSMAX+uGNSS.GLOMAX+uGNSS.GALMAX+uGNSS.BDSMAX+uGNSS.QZSMAX:
-        prn = sat-(uGNSS.GPSMAX+uGNSS.GLOMAX+uGNSS.GALMAX +
-                   uGNSS.BDSMAX+uGNSS.QZSMAX)+100
+    elif sat > uGNSS.SBSMIN:
+        prn = sat+100-uGNSS.SBSMIN
         sys = uGNSS.SBS
-    elif sat > uGNSS.GPSMAX+uGNSS.GLOMAX+uGNSS.GALMAX+uGNSS.BDSMAX:
-        prn = sat-(uGNSS.GPSMAX+uGNSS.GLOMAX+uGNSS.GALMAX+uGNSS.BDSMAX)+192
-        sys = uGNSS.QZS
-    elif sat > uGNSS.GPSMAX+uGNSS.GLOMAX+uGNSS.GALMAX:
-        prn = sat-(uGNSS.GPSMAX+uGNSS.GLOMAX+uGNSS.GALMAX)
+    elif sat > uGNSS.BDSMIN:
+        prn = sat-uGNSS.BDSMIN
         sys = uGNSS.BDS
-    elif sat > uGNSS.GPSMAX+uGNSS.GLOMAX:
-        prn = sat-(uGNSS.GPSMAX+uGNSS.GLOMAX)
-        sys = uGNSS.GAL
-    elif sat > uGNSS.GPSMAX:
-        prn = sat-uGNSS.GPSMAX
+    elif sat > uGNSS.GLOMIN:
+        prn = sat-uGNSS.GLOMIN
         sys = uGNSS.GLO
+    elif sat > uGNSS.QZSMIN:
+        prn = sat+192+uGNSS.QZSMIN
+        sys = uGNSS.QZS
+    elif sat > uGNSS.GALMIN:
+        prn = sat-uGNSS.GALMIN
+        sys = uGNSS.GAL
     else:
         prn = sat
         sys = uGNSS.GPS
@@ -789,6 +803,7 @@ def id2sat(id_):
     sat = prn2sat(sys, prn)
     return sat
 
+
 def char2sys(c):
     """ convert character to GNSS """
     gnss_tbl = {'G': uGNSS.GPS, 'R': uGNSS.GLO, 'E': uGNSS.GAL, 'C': uGNSS.BDS,
@@ -801,7 +816,7 @@ def char2sys(c):
 
 
 def sys2char(sys):
-    """ convert character to GNSS """
+    """ convert GNSS to character """
     gnss_tbl = {uGNSS.GPS: 'G', uGNSS.GLO: 'R', uGNSS.GAL: 'E', uGNSS.BDS: 'C',
                 uGNSS.QZS: 'J', uGNSS.SBS: 'S', uGNSS.IRN: 'I'}
 
@@ -811,6 +826,18 @@ def sys2char(sys):
         return gnss_tbl[sys]
 
 
+def sys2str(sys):
+    """ convert GNSS to string """
+    gnss_tbl = {uGNSS.GPS: 'GPS', uGNSS.GLO: 'GLONASS',
+                uGNSS.GAL: 'GALILEO', uGNSS.BDS: 'BEIDOU',
+                uGNSS.QZS: 'QZSS', uGNSS.SBS: 'SBAS', uGNSS.IRN: 'IRNSS'}
+
+    if sys not in gnss_tbl:
+        return "???"
+    else:
+        return gnss_tbl[sys]
+
+
 def vnorm(r):
     """ calculate norm of a vector """
     return r/np.linalg.norm(r)
@@ -903,6 +930,8 @@ def xyz2enu(pos):
 def enu2xyz(pos):
     """ return ENU to ECEF conversion matrix from LLH """
     return np.array(np.matrix(xyz2enu(pos)).I)
+
+
 def ecef2pos(r):
     """  ECEF to LLH position conversion """
     e2 = rCST.FE_WGS84*(2-rCST.FE_WGS84)
@@ -1066,3 +1095,37 @@ def tropmodel(t, pos, el=np.pi/2, humi=0.7):
     trop_hs = 0.0022768*pres/(1.0-0.00266*np.cos(2*pos[0])-0.00028e-3*hgt)
     trop_wet = 0.002277*(1255.0/temp+0.05)*e
     return trop_hs, trop_wet, z
+
+
+def meteoHpf():
+    """
+    Hopfield atmosphere model
+    """
+    pres = 1010.0
+    temp = 291.1
+    e = 10.4
+    return pres, temp, e
+
+
+def tropmapfHpf(el):
+    """
+    Tropospheric mapping functions for Hopfield model
+    """
+
+    mapfh = 1.0/(np.sin(np.sqrt(el**2+(np.pi/72.0)**2)))
+    mapfw = 1.0/(np.sin(np.sqrt(el**2+(np.pi/120.0)**2)))
+
+    return mapfh, mapfw
+
+
+def tropmodelHpf():
+    """
+    Hopfield tropospheric delay model
+    """
+    # standard atmosphere for Hopfield model
+    pres, temp, e = meteoHpf()
+    # Hopfield
+    trop_dry = (77.6e-6 * (-613.3768/temp+148.98)*pres)/5.0
+    trop_wet = (77.6e-6 * 11000.0 * 4810.0 * e/temp**2)/5.0
+
+    return trop_dry, trop_wet
diff --git a/src/cssrlib/peph.py b/src/cssrlib/peph.py
index 423cbdc..9b314be 100644
--- a/src/cssrlib/peph.py
+++ b/src/cssrlib/peph.py
@@ -14,7 +14,7 @@
 from cssrlib.rinex import rnxdec
 import numpy as np
 from math import pow, sin, cos
-import gzip
+
 
 NMAX = 10
 MAXDTE = 900.0
@@ -460,25 +460,27 @@ def readpcv(self, fname):
 def searchpcv(pcvs, name, time):
     """ get satellite or receiver antenna pcv """
 
-    if isinstance(name, int):
+    if isinstance(name, str):
+
         for pcv in pcvs:
-            if pcv.sat != name:
+            if pcv.type != name:
                 continue
             if pcv.ts.time != 0 and timediff(pcv.ts, time) > 0.0:
                 continue
             if pcv.te.time != 0 and timediff(pcv.te, time) < 0.0:
                 continue
             return pcv
+
     else:
+
         for pcv in pcvs:
-            if pcv.type != name:
+            if pcv.sat != name:
                 continue
             if pcv.ts.time != 0 and timediff(pcv.ts, time) > 0.0:
                 continue
             if pcv.te.time != 0 and timediff(pcv.te, time) < 0.0:
                 continue
             return pcv
-        return None
 
     return None
 
@@ -644,10 +646,12 @@ def antModelTx(nav, e, sigs, sat, time, rs):
     # Select satellite antenna
     #
     ant = searchpcv(nav.sat_ant, sat, time)
+    if ant is None:
+        return None
 
     # Zenit angle and zenit angle grid
     #
-    za = np.rad2deg(np.arccos(np.dot(ez,-e)))
+    za = np.rad2deg(np.arccos(np.dot(ez, -e)))
     za_t = np.arange(ant.zen[0], ant.zen[1]+ant.zen[2], ant.zen[2])
 
     # Interpolate PCV and map PCO on line-of-sight vector
diff --git a/src/cssrlib/plot.py b/src/cssrlib/plot.py
index 8c24e4d..192634c 100644
--- a/src/cssrlib/plot.py
+++ b/src/cssrlib/plot.py
@@ -46,7 +46,8 @@ def plot_elv(t, elv, elmask=0, satlist=None):
 
     tmax = t[-1]//60+1
     plt.ylabel('Elevation Angle [deg]')
-    plt.xlabel('Time [min]')    plt.grid()
+    plt.xlabel('Time [min]')
+    plt.grid()
     plt.axis([0, tmax, 0, 90])
     plt.show()
     return nsat
diff --git a/src/cssrlib/pntpos.py b/src/cssrlib/pntpos.py
index 6e0969a..6314bec 100644
--- a/src/cssrlib/pntpos.py
+++ b/src/cssrlib/pntpos.py
@@ -3,7 +3,7 @@
 """
 import numpy as np
 from cssrlib.gnss import rCST, ecef2pos, geodist, satazel, ionmodel,\
-    tropmodel, Nav, tropmapf, kfupdate, sat2prn
+    tropmodel, Nav, tropmapf, kfupdate
 from cssrlib.ephemeris import findeph, satposs
 
 
@@ -38,7 +38,13 @@ def stdinit():
 
 
 def rescode(obs, nav, rs, dts, svh, x):
-    """ calculate code residuals """
+    """
+    calculate code residuals
+
+    NOTE: currently only supports TGD application for GPS L1 P(Y) and GAL E1
+    NOTE: currently only supports iono delay modeling for GPS L1 and GAL E1
+    """
+
     n = obs.sat.shape[0]
     rr = x[0:3]
     dtr = x[nav.na]
@@ -48,34 +54,55 @@ def rescode(obs, nav, rs, dts, svh, x):
     azv = np.zeros(n)
     elv = np.zeros(n)
     nv = 0
+
     for i in range(n):
-        sys, _ = sat2prn(obs.sat[i])
-        if sys not in nav.gnss_t:
-            continue
+
+        # Check for valid satellite positions and health
+        #
         if np.linalg.norm(rs[i, :]) < rCST.RE_WGS84 or svh[i] > 0:
             continue
+
+        # Get geometric distance and line-of-sight vector
+        # (with Earth rotation correction)
+        #
         r, e = geodist(rs[i, :], rr)
         az, el = satazel(pos, e)
         if el < nav.elmin:
             continue
+
         eph = findeph(nav.eph, obs.t, obs.sat[i])
+
+        # Use first pseudorange in the dictionary (must be on L1/E1!)
+        #        
         if obs.P[i, 0] == 0:
             continue
+
+        # Apply group delay correction
+        #        
         P = obs.P[i, 0]-eph.tgd*rCST.CLIGHT
+        
+        # Apply ionospheric delay correction for L1/E1
+        #
         dion = ionmodel(obs.t, pos, az, el, nav.ion)
+
+        # Apply tropospheric correction
+        #        
         trop_hs, trop_wet, _ = tropmodel(obs.t, pos, el)
         mapfh, mapfw = tropmapf(obs.t, pos, el)
         dtrp = mapfh*trop_hs+mapfw*trop_wet
+        
         v[nv] = P-(r+dtr-rCST.CLIGHT*dts[i]+dion+dtrp)
         H[nv, 0:3] = -e
         H[nv, nav.na] = 1
         azv[nv] = az
         elv[nv] = el
         nv += 1
+        
     v = v[0:nv]
     H = H[0:nv, :]
     azv = azv[0:nv]
     elv = elv[0:nv]
+    
     return v, H, nv, azv, elv
 
 
diff --git a/src/cssrlib/ppp.py b/src/cssrlib/ppp.py
index f3f6b43..3987e64 100644
--- a/src/cssrlib/ppp.py
+++ b/src/cssrlib/ppp.py
@@ -1,9 +1,9 @@
 """
 module for PPP processing
 """
-
-import numpy as np
 import cssrlib.gnss as gn
+from cssrlib.peph import gpst2utc
+import numpy as np
 from math import sin, cos, atan2, asin
 
 
@@ -116,7 +116,7 @@ def nut_iau1980(t_, f):
         [-1,   0,   0,   1,   1,  -388.3,       1,    0.0,     0,   0.0],
         [-1,  -1,   0,   2,   1,    35.0,       1,    0.0,     0,   0.0],
         [0,   1,   0,   1,   0,    27.3,       1,    0.0,     0,   0.0]
-        ])
+    ])
     dpsi = 0
     deps = 0
     for i in range(106):
@@ -269,15 +269,27 @@ def shapiro(rsat, rrcv):
     return corr
 
 
-def windupcorr(time, rs, vs, rr, phw):
+def windupcorr(time, rs, vs, rr, phw, full=False):
     """ calculate windup correction """
-    we = np.array([0, 0, gn.rCST.OMGE])
-    # rsun, rmoon, gmst = sunmoonpos(gn.gpst2utc(time))
     ek = gn.vnorm(rr-rs)
-    ezs = gn.vnorm(-rs)
-    ess = gn.vnorm(vs+np.cross(we, rs))
-    eys = gn.vnorm(np.cross(ezs, ess))
-    exs = np.cross(eys, ezs)
+    if full:
+        # Satellite antenna frame unit vectors assuming standard yaw attitude law
+        #
+        rsun, _, _ = sunmoonpos(gpst2utc(time))
+        r = -rs
+        ezs = r/np.linalg.norm(r)
+        r = rsun-rs
+        ess = r/np.linalg.norm(r)
+        r = np.cross(ezs, ess)
+        eys = r/np.linalg.norm(r)
+        exs = np.cross(eys, ezs)
+    else:
+        we = np.array([0, 0, gn.rCST.OMGE])
+        ek = gn.vnorm(rr-rs)
+        ezs = gn.vnorm(-rs)
+        ess = gn.vnorm(vs+np.cross(we, rs))
+        eys = gn.vnorm(np.cross(ezs, ess))
+        exs = np.cross(eys, ezs)
     pos = gn.ecef2pos(rr)
     E = gn.xyz2enu(pos)
     exr = E[0, :]
diff --git a/src/cssrlib/pppigs.py b/src/cssrlib/pppigs.py
index 8c4d6d1..5d9a9de 100644
--- a/src/cssrlib/pppigs.py
+++ b/src/cssrlib/pppigs.py
@@ -1,109 +1,57 @@
 """
 module for standard PPP positioning
 """
-from collections import defaultdict
-from copy import copy
-from datetime import datetime
+
 import numpy as np
 import sys
+
 import cssrlib.gnss as gn
-from cssrlib.gnss import tropmodel, antmodel, uGNSS
-from cssrlib.ephemeris import satposs
-from cssrlib.cssrlib import sSigGPS, sSigGAL, sSigQZS
+from cssrlib.ephemeris import findeph
+from cssrlib.gnss import rCST, sat2id, sat2prn
+from cssrlib.gnss import timeadd, time2str
+from cssrlib.gnss import uTYP
 from cssrlib.ppp import tidedisp, shapiro, windupcorr
-from cssrlib.rtk import IB, ddres, resamb_lambda, valpos, holdamb, initx
-
-MAXITR = 10
-ELMIN = 10
-NX = 4
+from cssrlib.peph import antModelRx, antModelTx
+from cssrlib.rtk import IB, ddcov, resamb_lambda, valpos, holdamb, initx
+from cssrlib.rtk import varerr
 
 
-def logmon(nav, t, sat, cs, iu=None):
-    """ log variables for monitoring """
-    _, tow = gn.time2gpst(t)
-    if iu is None:
-        cpc = cs.cpc
-        prc = cs.prc
-        osr = cs.osr
-    else:
-        cpc = cs.cpc[iu, :]
-        prc = cs.prc[iu, :]
-        osr = cs.osr[iu, :]
-    if nav.loglevel >= 2:
-        n = cpc.shape
-        for i in range(n[0]):
-            if cpc[i, 0] == 0 and cpc[i, 1] == 0:
-                continue
-            sys, prn = gn.sat2prn(sat[i])
-            pb = osr[i, 0:2]
-            cb = osr[i, 2:4]
-            antr = osr[i, 4:6]
-            phw = osr[i, 6:8]
-            trop = osr[i, 8]
-            iono = osr[i, 9]
-            relatv = osr[i, 10]
-            dorb = osr[i, 11]
-            dclk = osr[i, 12]
-            # tow  sys  prn  pb1  pb2  cb1 cb2 trop iono  antr1  antr2  relatv
-            # wup1  wup2  CPC1  CPC2  PRC1  PRC2  dorb  dclk
-            nav.fout.write("%6d\t%2d\t%3d\t%8.3f\t%8.3f\t%8.3f\t%8.3f\t"
-                           % (tow, sys, prn, pb[0], pb[1], cb[0], cb[1]))
-            nav.fout.write("%8.3f\t%8.3f\t%8.3f\t%8.3f\t"
-                           % (trop, iono, antr[0], antr[1]))
-            nav.fout.write("%8.3f\t%8.3f\t%8.3f\t%8.3f\t%8.3f\t"
-                           % (relatv, phw[0], phw[1], cpc[i, 0], cpc[i, 1]))
-            nav.fout.write("%8.3f\t%8.3f\t%8.3f\t%8.3f\n"
-                           % (prc[i, 0], prc[i, 1], dorb, dclk))
-    return 0
+def IT(na):
+    """ return index of zenith tropospheric delay estimate """
+    return na-gn.uGNSS.MAXSAT-1
 
 
-def rtkinit(nav, pos0=np.zeros(3)):
-    """ initialize variables for RTK """
-
-    # Logging level
-    nav.monlevel = 1
-
-    # Number of frequencies
-    nav.nf = 2
+def II(s, na):
+    """ return index of slant ionospheric delay estimate """
+    idx = na-gn.uGNSS.MAXSAT+s-1
+    return idx
 
-    # Positioning mode
-    # 0:static, 1:kinematic
-    nav.pmode = 0
 
-    # Number of tracking channels (1 per satellite)
-    nav.nChan = 15
-    nav.satIdx = defaultdict(list)
-    for n in range(nav.nChan):
-        nav.satIdx[None].append(n)
+def ionoDelay(sig1, sig2, pr1, pr2):
+    """
+    Compute ionospheric delay based on dual-frequency pseudorange difference
+    """
+    f1_2 = sig1.frequency()**2
+    f2_2 = sig2.frequency()**2
+    return (pr1-pr2)/(1-f1_2/f2_2)
 
-    # State index
-    nav.idx_pos = -1
-    nav.idx_vel = -1
-    nav.idx_ztd = -1
-    nav.idx_ion = -1
 
-    # Position (and optional velocity) states
-    nav.na = 3 if nav.pmode == 0 else 6
-    nav.nq = 3 if nav.pmode == 0 else 6
+def rtkinit(nav, pos0=np.zeros(3)):
+    """ initialize variables for PPP """
 
-    # Zenith tropospheric delay state
-    nav.na += 1
-    nav.nq += 1
+    # Number of frequencies (actually signals!)
+    #
+    nav.nf = 2  # TODO: make obsolete if possible
 
-    # Slant ionospheric delay states
-    nav.na += nav.nChan
-    nav.nq += nav.nChan
+    # Position (+ optional velocity), zenith tropo delay and slant ionospheric delay states
+    #
+    nav.na = (4 if nav.pmode == 0 else 7) + gn.uGNSS.MAXSAT
+    nav.nq = (4 if nav.pmode == 0 else 7) + gn.uGNSS.MAXSAT
 
-    nav.ratio = 0
     nav.thresar = [2.0]
 
-    # Index of position, velocity and ztd
-    nav.idx_pos = 0
-    nav.idx_vel = -1 if nav.pmode == 0 else 3
-    nav.idx_ztd = 3 if nav.pmode == 0 else 6
-    nav.idx_ion = 4 if nav.pmode == 0 else 7
-
-    # state vector dimensions
+    # State vector dimensions (inlcuding slat iono delay and ambiguities)
+    #
     nav.nx = nav.na+gn.uGNSS.MAXSAT*nav.nf
 
     nav.x = np.zeros(nav.nx)
@@ -112,80 +60,84 @@ def rtkinit(nav, pos0=np.zeros(3)):
     nav.xa = np.zeros(nav.na)
     nav.Pa = np.zeros((nav.na, nav.na))
 
-    nav.nfix = nav.neb = 0
     nav.phw = np.zeros(gn.uGNSS.MAXSAT)
     nav.el = np.zeros(gn.uGNSS.MAXSAT)
 
-    # parameter for PPP
-
-    # observation noise parameters
-    nav.eratio = [50, 50]
-    nav.err = [0, 0.01, 0.005]/np.sqrt(2)
+    # Parameters for PPP
+    #
+    # Observation noise parameters
+    #
+    nav.eratio = [100, 100]
+    nav.err = [0, 0.001, 0.001]/np.sqrt(2)
+    #nav.eratio = [50, 50]
+    #nav.err = [0, 0.01, 0.005]/np.sqrt(2)
 
-    # initial sigma for state covariance
-    nav.sig_p0 = 30.0
+    # Initial sigma for state covariance
+    #
+    nav.sig_p0 = 100.0
     nav.sig_v0 = 1.0
-    nav.sig_ztd0 = 0.010
-    nav.sig_ion0 = 10.0
+    nav.sig_ztd0 = 0.25
+    nav.sig_ion0 = 1.0
     nav.sig_n0 = 30.0
 
-    # process noise sigma
-    nav.sig_qp = 0.01
-    nav.sig_qv = 1.0
-    nav.sig_qion = 10.0  # [m/sqrt(s)]
-    nav.sig_qztd = 0.1e-3/np.sqrt(30)  # [m/sqrt(s)]
+    # Process noise sigma
+    #
+    if nav.pmode == 0:
+        nav.sig_qp = 100.0/np.sqrt(1)  # [m/sqrt(s)]
+        nav.sig_qv = None
+    else:
+        nav.sig_qp = 0.01/np.sqrt(1)  # [m/sqrt(s)]
+        nav.sig_qv = 1.0/np.sqrt(1)  # [m/s/sqrt(s)]
+    nav.sig_qztd = 0.1/np.sqrt(3600)  # [m/sqrt(s)] -> 1 cm**2/h
+    nav.sig_qion = 1.0  # [m]
 
     nav.tidecorr = True
-    nav.armode = 3  # 1:continuous,2:instantaneous,3:fix-and-hold
-    nav.elmaskar = np.deg2rad(20)  # elevation mask for AR
-    nav.gnss_t = [uGNSS.GPS, uGNSS.GAL]
+    nav.armode = 3  # 0:float-ppp,1:continuous,2:instantaneous,3:fix-and-hold
+    nav.elmaskar = np.deg2rad(20.0)  # elevation mask for AR
+    nav.elmin = np.deg2rad(10.0)
 
     # Initial state vector
+    #
     nav.x[0:3] = pos0
     if nav.pmode >= 1:  # kinematic
         nav.x[3:6] = 0.0  # velocity
 
     # Diagonal elements of covariance matrix
+    #
     dP = np.diag(nav.P)
     dP.flags['WRITEABLE'] = True
 
     dP[0:3] = nav.sig_p0**2
+    # Velocity
     if nav.pmode >= 1:  # kinematic
         dP[3:6] = nav.sig_v0**2
-    dP[nav.idx_ztd] = nav.sig_ztd0**2
-    dP[nav.idx_ion:nav.idx_ion+nav.nChan] = nav.sig_ion0**2
+    # Tropo delay
+    if nav.pmode >= 1:  # kinematic
+        dP[6] = nav.sig_ztd0**2
+    else:
+        dP[3] = nav.sig_ztd0**2
 
     # Process noise
+    #
     nav.q = np.zeros(nav.nq)
+    nav.q[0:3] = nav.sig_qp**2
+    # Velocity
     if nav.pmode >= 1:  # kinematic
-        nav.q[0:3] = nav.sig_qp**2
         nav.q[3:6] = nav.sig_qv**2
+    # Tropo delay
+    if nav.pmode >= 1:  # kinematic
+        nav.q[6] = nav.sig_qztd**2
+    else:
+        nav.q[3] = nav.sig_qztd**2
+    # Iono delay
+    if nav.pmode >= 1:  # kinematic
+        nav.q[7:7+gn.uGNSS.MAXSAT] = nav.sig_qion**2
     else:
-        nav.q[0:3] = nav.sig_qp**2
-    nav.q[nav.idx_ztd] = nav.sig_qztd**2
-    nav.q[nav.idx_ion:nav.idx_ion+nav.nChan] = nav.sig_qion**2
-
-    # obs index
-    i0 = {gn.uGNSS.GPS: 0, gn.uGNSS.GAL: 0, gn.uGNSS.QZS: 0}
-    i1 = {gn.uGNSS.GPS: 1, gn.uGNSS.GAL: 2, gn.uGNSS.QZS: 1}
-    freq0 = {gn.uGNSS.GPS: nav.freq[0], gn.uGNSS.GAL: nav.freq[0],
-             gn.uGNSS.QZS: nav.freq[0]}
-    freq1 = {gn.uGNSS.GPS: nav.freq[1], gn.uGNSS.GAL: nav.freq[2],
-             gn.uGNSS.QZS: nav.freq[1]}
-    nav.obs_idx = [i0, i1]
-    nav.obs_freq = [freq0, freq1]
-    nav.cs_sig_idx = {gn.uGNSS.GPS: [sSigGPS.L1C, sSigGPS.L2W],
-                      gn.uGNSS.GAL: [sSigGAL.L1X, sSigGAL.L5X],
-                      gn.uGNSS.QZS: [sSigQZS.L1C, sSigQZS.L2X]}
-
-    nav.fout = None
-    nav.logfile = 'log.txt'
-    if nav.loglevel >= 2:
-        nav.fout = open(nav.logfile, 'w')
-        nav.fout.write(
-            "# tow\tsys\tprn\tpb1\tpb2\tcb1\tcb2\ttrop\tiono\tantr1\tantr2")
-        nav.fout.write(
-            "# \trelatv\twup1\twup2\tCPC1\tCPC2\tPRC1\tPRC2\tdorb\tdclk\n")
+        nav.q[4:4+gn.uGNSS.MAXSAT] = nav.sig_qion**2
+
+    # Logging level
+    #
+    nav.monlevel = 2
 
 
 def sysidx(satlist, sys_ref):
@@ -198,19 +150,9 @@ def sysidx(satlist, sys_ref):
     return idx
 
 
-def varerr(nav, el, f):
-    """ variation of measurement """
-    s_el = np.sin(el)
-    if s_el <= 0.0:
-        return 0.0
-    fact = nav.eratio[f-nav.nf] if f >= nav.nf else 1
-    a = fact*nav.err[1]
-    b = fact*nav.err[2]
-    return 2.0*(a**2+(b/s_el)**2)
-
-
-def udstate(nav, obs, cs):
+def udstate(nav, obs):
     """ time propagation of states and initialize """
+
     tt = gn.timediff(obs.t, nav.t)
 
     ns = len(obs.sat)
@@ -220,86 +162,158 @@ def udstate(nav, obs, cs):
         sys_i, _ = gn.sat2prn(sat_i)
         sys.append(sys_i)
 
-    # pos,vel,ztd,ion
-    na = nav.na
-    Phi = np.eye(na)
-    dPhi = np.diag(Phi)
-    dPhi.flags['WRITEABLE'] = True
-    dPhi[na-nav.nChan:na] = 0.0
+    # pos,vel,ztd,ion,amb
+    #
+    nx = nav.nx
+    ni = nav.na-gn.uGNSS.MAXSAT
+    Phi = np.eye(nx)
+    Phi[ni:nav.na, ni:nav.na] = np.zeros((gn.uGNSS.MAXSAT, gn.uGNSS.MAXSAT))
     if nav.pmode > 0:
+        nav.x[0:3] += nav.x[3:6]*tt
         Phi[0:3, 3:6] = np.eye(3)*tt
-    nav.x[0:na] = Phi@nav.x[0:na]
-    nav.P[0:na, 0:na] = Phi@nav.P[0:na, 0:na]@Phi.T
+    nav.P[0:nx, 0:nx] = Phi@nav.P[0:nx, 0:nx]@Phi.T
 
     # Process noise
+    #
     dP = np.diag(nav.P)
     dP.flags['WRITEABLE'] = True
     dP[0:nav.nq] += nav.q[0:nav.nq]*tt
 
-    # bias
+    # Carrier-phase ambiguity
+    #
     for f in range(nav.nf):
-        # reset phase-bias if instantaneous AR or
-        # expire obs outage counter
+
+        # Reset phase-ambiguity if instantaneous AR or expire obs outage counter
+        #
         for i in range(gn.uGNSS.MAXSAT):
+
             sat_ = i+1
             nav.outc[i, f] += 1
             reset = (nav.outc[i, f] > nav.maxout)
             sys_i, _ = gn.sat2prn(sat_)
-            if sys_i not in nav.gnss_t:
+            if sys_i not in obs.sig.keys():
                 continue
+
+            # Reset ambiguity estimate
+            #
             j = IB(sat_, f, nav.na)
             if reset and nav.x[j] != 0.0:
                 initx(nav, 0.0, 0.0, j)
                 nav.outc[i, f] = 0
-        # cycle slip check by LLI
+
+                if nav.monlevel > 0:
+                    print("{}  {} - reset ambiguity  {}"
+                          .format(time2str(obs.t), sat2id(sat_),
+                                  obs.sig[sys_i][uTYP.L][f]))
+
+            # Reset slant ionospheric delay estimate
+            #
+            j = II(sat_, nav.na)
+            if reset and nav.x[j] != 0.0:
+                initx(nav, 0.0, 0.0, j)
+
+                if nav.monlevel > 0:
+                    print("{}  {} - reset ionosphere"
+                          .format(time2str(obs.t), sat2id(sat_)))
+
+        # Cycle  slip check by LLI
+        #
         for i in range(ns):
-            if sys[i] not in nav.gnss_t:
-                continue
-            j = nav.obs_idx[f][sys[i]]
-            if obs.lli[i, j] & 1 == 0:
+            if sys[i] not in obs.sig.keys():
                 continue
-            initx(nav, 0.0, 0.0, IB(sat[i], f, nav.na))
-        # reset bias if correction is not available
-        for i in range(ns):
-            if sat[i] in cs.sat_n:
+            if obs.lli[i, f] & 1 == 0:
                 continue
             initx(nav, 0.0, 0.0, IB(sat[i], f, nav.na))
-        # bias
+
+        # Ambiguity
+        #
         bias = np.zeros(ns)
+        ion = np.zeros(ns)
+
         offset = 0
         na = 0
         for i in range(ns):
-            if sys[i] not in nav.gnss_t:
+
+            if sys[i] not in obs.sig.keys():
                 continue
-            j = nav.obs_idx[f][sys[i]]
-            freq = nav.obs_freq[f][sys[i]]
-            cp = obs.L[i, j]
-            pr = obs.P[i, j]
-            if cp == 0.0 or pr == 0.0 or freq == 0.0:
+
+            # Get dual-frequency pseudoranges for this constellation
+            #
+            sig1 = obs.sig[sys[i]][uTYP.C][0]
+            sig2 = obs.sig[sys[i]][uTYP.C][1]
+
+            pr1 = obs.P[i, 0]
+            pr2 = obs.P[i, 1]
+
+            # Skip zero observations
+            #
+            if pr1 == 0.0 or pr2 == 0.0:
+                continue
+
+            # Get iono delay at frequency of first signal
+            #
+            ion[i] = ionoDelay(sig1, sig2, pr1, pr2)
+
+            # Get pseudorange and carrier-phase observation of signal f
+            #
+            sig = obs.sig[sys[i]][uTYP.L][f]
+            lam = sig.wavelength()
+
+            cp = obs.L[i, f]
+            pr = obs.P[i, f]
+            if cp == 0.0 or pr == 0.0 or lam is None:
                 continue
-            bias[i] = cp-pr*freq/gn.rCST.CLIGHT
+
+            bias[i] = cp - pr/lam + 2.0 * \
+                (sig1.frequency()**2/sig.frequency()**2)*ion[i]/lam
+
             amb = nav.x[IB(sat[i], f, nav.na)]
             if amb != 0.0:
-                offset += bias[i]-amb
+                offset += bias[i] - amb
                 na += 1
-        # adjust phase-code coherency
+
+        # Adjust phase-code coherency
+        #
         if na > 0:
             db = offset/na
             for i in range(gn.uGNSS.MAXSAT):
                 if nav.x[IB(i+1, f, nav.na)] != 0.0:
                     nav.x[IB(i+1, f, nav.na)] += db
-        # initialize ambiguity
+
+        # Initialize ambiguity
+        #
         for i in range(ns):
+
+            sys_i, _ = sat2prn(sat[i])
+
             j = IB(sat[i], f, nav.na)
-            if bias[i] == 0.0 or nav.x[j] != 0.0:
-                continue
-            initx(nav, bias[i], nav.sig_n0**2, j)
+            if bias[i] != 0.0 and nav.x[j] == 0.0:
+
+                initx(nav, bias[i], nav.sig_n0**2, j)
+
+                if nav.monlevel > 0:
+                    sig = obs.sig[sys_i][uTYP.L][f]
+                    print("{}  {} - init  ambiguity  {} {:12.3f}"
+                          .format(time2str(obs.t), sat2id(sat[i]), sig, bias[i]))
+
+            j = II(sat[i], nav.na)
+            if ion[i] != 0 and nav.x[j] == 0.0:
+
+                initx(nav, ion[i], nav.sig_ion0**2, j)
+
+                if nav.monlevel > 0:
+                    print("{}  {} - init  ionosphere      {:12.3f}"
+                          .format(time2str(obs.t), sat2id(sat[i]), ion[i]))
+
     return 0
 
 
-def zdres(nav, obs, rs, vs, dts, svh, rr, cs):
+def zdres(nav, obs, bsx, rs, vs, dts, svh, rr):
     """ non-differential residual """
+
     _c = gn.rCST.CLIGHT
+    ns2m = _c*1e-9
+
     nf = nav.nf
     n = len(obs.P)
     y = np.zeros((n, nf*2))
@@ -307,243 +321,289 @@ def zdres(nav, obs, rs, vs, dts, svh, rr, cs):
     e = np.zeros((n, 3))
     rr_ = rr.copy()
 
-    # Tide corrections
+    # Solid Earth tide corrections
+    #
+    # TODO: add solid earth tide displacements
+    #
     if nav.tidecorr:
         pos = gn.ecef2pos(rr_)
         disp = tidedisp(gn.gpst2utc(obs.t), pos)
-        rr_ += disp
+    else:
+        disp = np.zeros(3)
+    rr_ += disp
 
     # Geodetic position
+    #
     pos = gn.ecef2pos(rr_)
-    inet = cs.find_grid_index(pos)
 
-    # Tropospheric dry and wet delays at position
-    trop_hs, trop_wet, _ = tropmodel(obs.t, pos)
+    # Zenith tropospheric dry and wet delays at user position
+    #
+    trop_hs, trop_wet, _ = gn.tropmodel(obs.t, pos)
 
-    cs.cpc = np.zeros((n, nf))
-    cs.prc = np.zeros((n, nf))
-    cs.osr = np.zeros((n, 4*nf+5))
+    cpc = np.zeros((n, nf))
+    prc = np.zeros((n, nf))
 
     for i in range(n):
+
         sat = obs.sat[i]
         sys, _ = gn.sat2prn(sat)
-        if svh[i] > 0 or sys not in nav.gnss_t or sat in nav.excl_sat:
-            continue
-        if sat not in cs.lc[inet].sat_n:
-            continue
-        idx_n = np.where(cs.sat_n == sat)[0][0]  # global
-        idx_l = np.where(cs.lc[inet].sat_n == sat)[0][0]  # local
-        kidx = [-1]*nav.nf
-        nsig = 0
-        for k, sig in enumerate(cs.sig_n[idx_n]):
-            for f in range(nav.nf):
-                if sig == nav.cs_sig_idx[sys][f]:
-                    kidx[f] = k
-                    nsig += 1
-        if nsig < nav.nf:
+
+        if svh[i] > 0 or sys not in obs.sig.keys() or sat in nav.excl_sat:
             continue
-        # check for measurement consistency
+
+        # Check for measurement consistency
+        #
         flg_m = True
         for f in range(nav.nf):
-            k = nav.obs_idx[f][sys]
-            if obs.P[i, k] == 0.0 or obs.L[i, k] == 0.0 \
-                    or obs.lli[i, k] == 1:
+
+            if obs.P[i, f] == 0.0 or obs.L[i, f] == 0.0 or obs.lli[i, f] == 1:
                 flg_m = False
+
+        # Check C/N0 at first frequency
+        #
         if obs.S[i, 0] < nav.cnr_min:
             flg_m = False
+
         if flg_m is False:
             continue
 
+        # Geometric distance corrected for Earth rotation during flight time
+        #
         r, e[i, :] = gn.geodist(rs[i, :], rr_)
         _, el[i] = gn.satazel(pos, e[i, :])
         if el[i] < nav.elmin:
             continue
 
-        freq = np.zeros(nav.nf)
-        lam = np.zeros(nav.nf)
-        for f in range(nav.nf):
-            freq[f] = nav.obs_freq[f][sys]
-            lam[f] = gn.rCST.CLIGHT/freq[f]
-
-        # global/local signal bias
-        cbias = np.zeros(nav.nf)
-        pbias = np.zeros(nav.nf)
-        # t1 = timediff(obs.t, cs.lc[0].t0[sCType.ORBIT])
-        # t2 = timediff(obs.t, cs.lc[inet].t0[sCType.PBIAS])
-
-        if cs.lc[0].cbias is not None:
-            cbias += cs.lc[0].cbias[idx_n][kidx]
-        if cs.lc[0].pbias is not None:
-            pbias += cs.lc[0].pbias[idx_n][kidx]
-        if cs.lc[inet].cbias is not None:
-            cbias += cs.lc[inet].cbias[idx_l][kidx]
-        if cs.lc[inet].pbias is not None:
-            pbias += cs.lc[inet].pbias[idx_l][kidx]
-            # if t1 >= 0 and t1 < 30 and t2 >= 30:
-            #     pbias += nav.dsis[sat]
-
-        # relativity effect
+        # Pseudorange and carrier-phase signals
+        #
+        sigsPR = obs.sig[sys][gn.uTYP.C]
+        sigsCP = obs.sig[sys][gn.uTYP.L]
+
+        # Code and phase signal bias [ns]
+        #
+        cbias = np.array([bsx.getosb(sat, obs.t, s)[0]*ns2m for s in sigsPR])
+        pbias = np.array([bsx.getosb(sat, obs.t, s)[0]*ns2m for s in sigsCP])
+
+        # Shapipo relativistic effect
+        #
         relatv = shapiro(rs[i, :], rr_)
 
-        # tropospheric delay mapping functions
+        # Tropospheric delay mapping functions
+        #
         mapfh, mapfw = gn.tropmapf(obs.t, pos, el[i])
 
-        # tropospheric delay mapping functions
+        # Tropospheric delay
+        #
         trop = mapfh*trop_hs + mapfw*trop_wet
 
-        # phase wind-up effect
+        # Phase wind-up effect
+        #
         nav.phw[sat-1] = windupcorr(obs.t, rs[i, :], vs[i, :], rr_,
-                                    nav.phw[sat-1])
+                                    nav.phw[sat-1], full=True)
+
+        # Wavelength
+        #
+        lam = np.array([s.wavelength() for s in sigsCP])
         phw = lam*nav.phw[sat-1]
-        antr = antmodel(nav, el[i], nav.nf)
-        # range correction
-        prc_c = trop+relatv+antr
-        # prc_c += nav.dorb[sat]-nav.dclk[sat]
-        cs.prc[i, :] = prc_c+cbias
-        cs.cpc[i, :] = prc_c+pbias+phw
-        cs.osr[i, :] = [pbias[0], pbias[1], cbias[0], cbias[1],
-                        antr[0], antr[1], phw[0], phw[1],
-                        trop, None, relatv, nav.dorb[sat], nav.dclk[sat]]
-        r += -_c*dts[i]
 
-        for f in range(nf):
-            k = nav.obs_idx[f][sys]
-            y[i, f] = obs.L[i, k]*lam[f]-(r+cs.cpc[i, f])
-            y[i, f+nf] = obs.P[i, k]-(r+cs.prc[i, f])
+        # Receiver/satellite antenna offset
+        #
+        antrPR = antModelRx(nav, pos, e[i, :], sigsPR)
+        antrCP = antModelRx(nav, pos, e[i, :], sigsCP)
+        antsPR = antModelTx(nav, e[i, :], sigsPR, sat, obs.t, rs[i, :])
+        antsCP = antModelTx(nav, e[i, :], sigsCP, sat, obs.t, rs[i, :])
 
-    return y, e, el
+        # Range correction
+        #
+        prc[i, :] = trop + antsPR + antrPR + cbias
+        cpc[i, :] = trop + antsCP + antrCP + pbias + phw
 
+        r += relatv - _c*dts[i]
 
-def ddcov(nb, n, Ri, Rj, nv):
-    """ DD measurement error covariance """
-    R = np.zeros((nv, nv))
-    k = 0
-    for b in range(n):
-        for i in range(nb[b]):
-            for j in range(nb[b]):
-                R[k+i, k+j] = Ri[k+i]
-                if i == j:
-                    R[k+i, k+j] += Rj[k+i]
-        k += nb[b]
-    return R
+        for f in range(nf):
+            y[i, f] = obs.L[i, f]*lam[f]-(r+cpc[i, f])
+            y[i, f+nf] = obs.P[i, f]-(r+prc[i, f])
+
+    return y, e, el
 
 
-def ddres(nav, t, x, y, e, sat, el, log=False):
+def sdres(nav, obs, x, y, e, sat, el):
     """
-    SD/DD phase/code residuals
-
-        nav :
-        t   :
-        x   :
-        y   :
-        e   :
-        sat :
-        el  :
+    SD phase/code residuals
+
+    Parameters
+    ----------
+
+    nav : Nav()
+        Auxiliary data structure
+    obs : Obs()
+        Data structure with observations
+    x   :
+        State vector elements
+    y   :
+        Un-differenced corrected observations
+    e   :
+        Line-of-sight vectors
+    sat : np.array of int
+        List of satellites
+    el  : np.array of float values
+        Elevation angles
+
+    Returns
+    -------
+    v   : np.array of float values
+        Residuals of single-difference measurements
+    H   : np.array of float values
+        Jacobian matrix with partial derivatives of state variables
+    R   : np.array of float values
+        Covariance matrix of single-difference measurements
     """
-    _c = gn.rCST.CLIGHT
+
     nf = nav.nf
     ns = len(el)
-    mode = 1 if len(y) == ns else 0  # 0:DD,1:SD
-    nb = np.zeros(2*len(nav.gnss_t)*nf, dtype=int)
+
+    nb = np.zeros(2*len(obs.sig.keys())*nf, dtype=int)
+
     Ri = np.zeros(ns*nf*2)
     Rj = np.zeros(ns*nf*2)
+
     nv = 0
     b = 0
+
     H = np.zeros((ns*nf*2, nav.nx))
     v = np.zeros(ns*nf*2)
 
     # Geodetic position
+    #
     pos = gn.ecef2pos(x[0:3])
 
-    idx_f = [0, 1]
-    for sys in nav.gnss_t:
+    # Loop over constellations
+    #
+    for sys in obs.sig.keys():
 
-        for f in range(nf):
-            idx_f[f] = nav.obs_idx[f][sys]
+        # Slant ionospheric delay reference frequency
+        #
+        freq0 = obs.sig[sys][uTYP.L][0].frequency()
 
         # Loop over twice the number of frequencies
         #   first for all carrier-phase observations
         #   second all pseudorange observations
+        #
         for f in range(0, nf*2):
 
             # Select carrier-phase frequency and iono frequency ratio
-            freq = nav.freq[idx_f[f % 2]]
+            #
             if f < nf:  # carrier
-                mu = -(nav.freq[idx_f[0]]/freq)**2
+                sig = obs.sig[sys][uTYP.L][f]
+                mu = -(freq0/sig.frequency())**2
             else:  # code
-                mu = +(nav.freq[idx_f[0]]/freq)**2
+                sig = obs.sig[sys][uTYP.C][f % 2]
+                mu = +(freq0/sig.frequency())**2
 
-            # Select reference satellite
+            # Select satellites from one constellation only
+            #
             idx = sysidx(sat, sys)
 
-            if log and 1 == 0:
-                for i in idx:
-                    print(sys, sat[i])
-
+            # Select reference satellite with highest elevation
+            #
             if len(idx) > 0:
                 i = idx[np.argmax(el[idx])]
 
-            # Loop over satellite
+            # Loop over satellites
+            #
             for j in idx:
 
-                # Skip reference satellite
+                # Skip reference satellite i
+                #
                 if i == j:
                     continue
 
+                # Skip invalid signals
+                #
                 if y[i, f] == 0.0 or y[j, f] == 0.0:
                     continue
 
-                #  DD residual
-                if mode == 0:
-                    if y[i+ns, f] == 0.0 or y[j+ns, f] == 0.0:
-                        continue
-                    v[nv] = (y[i, f]-y[i+ns, f])-(y[j, f]-y[j+ns, f])
-                else:
-                    v[nv] = y[i, f]-y[j, f]
+                #  Single-difference measurement
+                #
+                v[nv] = y[i, f] - y[j, f]
 
                 # SD line-of-sight vectors
-                H[nv, 0:3] = -e[i, :]+e[j, :]
+                #
+                H[nv, 0:3] = -e[i, :] + e[j, :]
 
                 # SD troposphere
-                _, mapfwi = gn.tropmapf(t, pos, el[i])
-                _, mapfwj = gn.tropmapf(t, pos, el[j])
-                H[nv, nav.idx_ztd] = (mapfwi-mapfwj)
-                v[nv] -= (mapfwi-mapfwj)*x[nav.idx_ztd]
+                #
+                _, mapfwi = gn.tropmapf(obs.t, pos, el[i])
+                _, mapfwj = gn.tropmapf(obs.t, pos, el[j])
 
-                # SD ionosphere
-                idx_i = nav.satIdx[sat[i]] + nav.idx_ion
-                idx_j = nav.satIdx[sat[j]] + nav.idx_ion
-                H[nv, idx_i] = mu
-                H[nv, idx_j] = mu
-                v[nv] -= mu*(x[idx_i]-x[idx_j])
-
-                if log:
-                    print("{:%D %T} {:10.3f} {:10.3f} {:10.3f} "
-                          .format(datetime.utcfromtimestamp(t.time),
+                idx_i = IT(nav.na)
+                H[nv, idx_i] = mapfwi - mapfwj
+                v[nv] -= (mapfwi - mapfwj)*x[idx_i]
+
+                if nav.monlevel > 1:
+                    print("{} ztd ({:3d},{:3d}) {:10.3f} {:10.3f} {:10.3f} "
+                          .format(time2str(obs.t), idx_i, idx_i,
                                   (mapfwi-mapfwj),
-                                  nav.x[nav.idx_ztd],
-                                  np.sqrt(nav.P[nav.idx_ztd, nav.idx_ztd])))
+                                  x[IT(nav.na)],
+                                  np.sqrt(nav.P[IT(nav.na), IT(nav.na)])))
+
+                # SD ionosphere
+                #
+                idx_i = II(sat[i], nav.na)
+                idx_j = II(sat[j], nav.na)
+                H[nv, idx_i] = +mu
+                H[nv, idx_j] = -mu
+                v[nv] -= mu*(x[idx_i] - x[idx_j])
+
+                if nav.monlevel > 1:
+                    print("{} {}-{} ion ({:3d},{:3d}) {:10.3f} {:10.3f} {:10.3f} {:10.3f} {:10.3f}"
+                          .format(time2str(obs.t),
+                                  sat2id(sat[i]), sat2id(sat[j]),
+                                  idx_i, idx_j,
+                                  mu,
+                                  x[idx_i], x[idx_j],
+                                  np.sqrt(nav.P[idx_i, idx_i]),
+                                  np.sqrt(nav.P[idx_j, idx_j])))
 
                 # SD ambiguity
+                #
                 if f < nf:  # carrier-phase
+
                     idx_i = IB(sat[i], f, nav.na)
                     idx_j = IB(sat[j], f, nav.na)
-                    lami = _c/freq
+
+                    lami = sig.wavelength()
+
                     v[nv] -= lami*(x[idx_i]-x[idx_j])
-                    H[nv, idx_i] = lami
+
+                    H[nv, idx_i] = +lami
                     H[nv, idx_j] = -lami
+
                     Ri[nv] = varerr(nav, el[i], f)  # measurement variance
                     Rj[nv] = varerr(nav, el[j], f)  # measurement variance
+
                     nav.vsat[sat[i]-1, f] = 1
                     nav.vsat[sat[j]-1, f] = 1
+
+                    if nav.monlevel > 1:
+                        print("{} {}-{} amb ({:3d},{:3d}) {:10.3f} {:10.3f} {:10.3f} {:10.3f} {:10.3f}"
+                              .format(time2str(obs.t),
+                                      sat2id(sat[i]), sat2id(sat[j]),
+                                      idx_i, idx_j,
+                                      lami,
+                                      x[idx_i], x[idx_j],
+                                      np.sqrt(nav.P[idx_i, idx_i]),
+                                      np.sqrt(nav.P[idx_j, idx_j])))
+
                 else:  # pseudorange
+
                     Ri[nv] = varerr(nav, el[i], f)  # measurement variance
                     Rj[nv] = varerr(nav, el[j], f)  # measurement variance
 
                 nb[b] += 1
-                nv += 1
+                nv += 1  # counter for single-difference observations
 
-            b += 1
+            b += 1  # counter for satellite
 
     v = np.resize(v, nv)
     H = np.resize(H, (nv, nav.nx))
@@ -563,91 +623,139 @@ def kfupdate(x, P, H, v, R):
     return x, P, S
 
 
-def pppigspos(nav, obs, cs):
-    """ PPP positioning with IGS files and conventions"""
+def satpreposs(obs, nav, orb):
+    """
+    Calculate pos/vel/clk for observed satellites
+
+    The satellite position, velocity and clock offset are computed at
+    transmission epoch. The signal time-of-flight is computed from a pseudorange
+    measurement corrected by the satellite clock offset, hence the observations
+    are required at this stage. The satellite clock is already corrected for the
+    relativistic effects. The satellite health indicator is extracted from the
+    broadcast navigation message.
+
+    Parameters
+    ----------
+    obs : Obs()
+        contains GNSS measurments
+    nav : Nav()
+        contains coarse satellite orbit and clock offset information
+    obs : peph()
+        contains precise satellite orbit and clock offset information
+
+    Returns
+    -------
+    rs  : np.array() of float
+        satellite position in ECEF [m]
+    vs  : np.array() of float
+        satellite velocities in ECEF [m/s]
+    dts : np.array() of float
+        satellite clock offsets [s]
+    svh : np.array() of int
+        satellite health code [-]
+    """
+
+    n = obs.sat.shape[0]
+    rs = np.zeros((n, 6))
+    dts = np.zeros((n, 2))
+    svh = np.zeros(n, dtype=int)
+
+    for i in range(n):
+
+        sat = obs.sat[i]
+        sys, _ = sat2prn(sat)
+
+        pr = obs.P[i, 0]
+        if pr == 0.0:
+            continue
+
+        t = timeadd(obs.t, -pr/rCST.CLIGHT)
+
+        rs[i, :], dts[i, :], _ = orb.peph2pos(t, sat, nav)
+
+        t = timeadd(t, -dts[i, 0])
+        rs[i, :], dts[i, :], _ = orb.peph2pos(t, sat, nav)
+
+        eph = findeph(nav.eph, t, sat)
+        if eph is None:
+            svh[i] = 1
+            continue
+        svh[i] = eph.svh
+
+    return rs[:, 0:3], rs[:, 3:6], dts[:, 0], svh
+
+
+def pppigspos(nav, obs, orb, bsx):
+    """
+    PPP positioning with IGS files and conventions
+    """
+
+    # Skip empty epochs
+    #
+    if len(obs.sat) == 0:
+        return
 
     # GNSS satellite positions, velocities and clock offsets
-    rs, vs, dts, svh = satposs(obs, nav, cs)
+    #
+    rs, vs, dts, svh = satpreposs(obs, nav, orb)
 
-    # Kalman filter time propagation
-    udstate(nav, obs, cs)
+    # Kalman filter time propagation, initialization of ambiguities and iono
+    #
+    udstate(nav, obs)
 
-    xa = np.zeros(nav.nx)
+    #xa = np.zeros(nav.nx)
     xp = nav.x.copy()
 
-    # Non-differential residuals for rover
-    yu, eu, elu = zdres(nav, obs, rs, vs, dts, svh, xp[0:3], cs)
+    # Non-differential residuals
+    #
+    yu, eu, elu = zdres(nav, obs, bsx, rs, vs, dts, svh, xp[0:3])
 
     # Select satellites above minimum elevation
+    #
     iu = np.where(elu >= nav.elmin)[0]
     sat = obs.sat[iu]
     y = yu[iu, :]
     e = eu[iu, :]
     el = elu[iu]
+
+    # Store reduced satellite list
+    #
     nav.sat = sat
     nav.el[sat-1] = el
     nav.y = y
     ns = len(sat)
 
-    # Update satellite index map
-    #
-    # NOTE: currently iono and ambiguity parameters are not newly initialized
-    #       on constellation changes! Must be added here.
-    #
-    # Remove missing satellites
-    satsIdx_ = copy(nav.satIdx)
-    for s in satsIdx_.keys():
-        # Skip void
-        if s not in sat and s is not None:
-            idx = nav.satIdx.pop(s)
-            nav.satIdx[None].append(idx)
-            print("Removing {} at {}".format(s, idx))
-    # Add new satellites
-    for s in sat:
-        if s not in nav.satIdx.keys():
-            idx = nav.satIdx.pop(None)
-            if len(idx) > 0:
-                nav.satIdx[s] = idx[0]
-                nav.satIdx[None] = idx[1:]
-                print("Adding   {} at {}".format(s, idx[0]))
-            else:
-                print("ERROR: satellite index full!")
-                sys.exit(1)
-
-    print()
-    print("Satellite index")
-    for s, n in nav.satIdx.items():
-        print(s, n)
-    print()
-
-    if nav.loglevel > 1:
-        logmon(nav, obs.t, sat, cs, iu)
-
-    # ???
+    # Check if observations of at least 6 satellites are left over after editing
+    #
     ny = y.shape[0]
     if ny < 6:
         nav.P[np.diag_indices(3)] = 1.0
         nav.smode = 5
         return -1
 
-    # DD residual
-    v, H, R = ddres(nav, obs.t, xp, y, e, sat, el)
+    # SD residuals
+    #
+    # NOTE: where are working on a reduced list of observations from here on
+    #
+    v, H, R = sdres(nav, obs, xp, y, e, sat, el)
     Pp = nav.P.copy()
 
     # Kalman filter measurement update
+    #
     xp, Pp, _ = kfupdate(xp, Pp, H, v, R)
 
-    # non-differential residual for rover after measurement update
-    yu, eu, elu = zdres(nav, obs, rs, vs, dts, svh, xp[0:3], cs)
+    # Non-differential residuals after measurement update
+    #
+    yu, eu, elu = zdres(nav, obs, bsx, rs, vs, dts, svh, xp[0:3])
     y = yu[iu, :]
     e = eu[iu, :]
     ny = y.shape[0]
-    nav.y2 = y
     if ny < 6:
         return -1
 
-    # residual for float solution
-    v, H, R = ddres(nav, obs.t, xp, y, e, sat, el)
+    # Residuals for float solution
+    #
+    v, H, R = sdres(nav, obs, xp, y, e, sat, el)
     if valpos(nav, v, R):
         nav.x = xp
         nav.P = Pp
@@ -664,20 +772,21 @@ def pppigspos(nav, obs, cs):
     else:
         nav.smode = 0
 
-    nav.smode = 5  # 5: float ambiguities
-    try:
-        nb, xa = resamb_lambda(nav, sat)
-    except:
-        nb = 0
+    nav.smode = 5  # 4: fixed ambiguities, 5: float ambiguities
+    nb, xa = resamb_lambda(nav, sat)
     if nb > 0:
-        yu, eu, elu = zdres(nav, obs, rs, vs, dts, svh, xa[0:3], cs)
+        yu, eu, elu = zdres(nav, obs, bsx, rs, vs, dts, svh, xa[0:3])
         y = yu[iu, :]
         e = eu[iu, :]
-        v, H, R = ddres(nav, obs.t, xa, y, e, sat, el)
-        if valpos(nav, v, R):       # R <= Q=H'PH+R  chisq<max_inno[3] (0.5)
+        v, H, R = sdres(nav, obs, xa, y, e, sat, el)
+        # R <= Q=H'PH+R  chisq<max_inno[3] (0.5)
+        if valpos(nav, v, R):
             if nav.armode == 3:     # fix and hold
                 holdamb(nav, xa)    # hold fixed ambiguity
             nav.smode = 4           # fix
 
+    # Store epoch for solution
+    #
     nav.t = obs.t
+
     return 0
diff --git a/src/cssrlib/ppprtk.py b/src/cssrlib/ppprtk.py
index e23a3cb..ed788a4 100644
--- a/src/cssrlib/ppprtk.py
+++ b/src/cssrlib/ppprtk.py
@@ -4,9 +4,10 @@
 
 import numpy as np
 import cssrlib.gnss as gn
-from cssrlib.gnss import tropmodel, antmodel, uGNSS
+from cssrlib.gnss import tropmodel
 from cssrlib.ephemeris import satposs
 from cssrlib.cssrlib import sSigGPS, sSigGAL, sSigQZS
+from cssrlib.peph import antModelRx
 from cssrlib.ppp import tidedisp, shapiro, windupcorr
 from cssrlib.rtk import IB, ddres, resamb_lambda, valpos, holdamb, initx
 
@@ -62,14 +63,14 @@ def rtkinit(nav, pos0=np.zeros(3)):
 
     nav.na = 3 if nav.pmode == 0 else 6
     nav.nq = 3 if nav.pmode == 0 else 6
-    nav.ratio = 0
+
     nav.thresar = [2.0]
     nav.nx = nav.na+gn.uGNSS.MAXSAT*nav.nf
     nav.x = np.zeros(nav.nx)
     nav.P = np.zeros((nav.nx, nav.nx))
     nav.xa = np.zeros(nav.na)
     nav.Pa = np.zeros((nav.na, nav.na))
-    nav.nfix = nav.neb = 0
+
     nav.phw = np.zeros(gn.uGNSS.MAXSAT)
     nav.el = np.zeros(gn.uGNSS.MAXSAT)
 
@@ -84,7 +85,6 @@ def rtkinit(nav, pos0=np.zeros(3)):
     nav.tidecorr = True
     nav.armode = 1  # 1:contunous,2:instantaneous,3:fix-and-hold
     nav.elmaskar = np.deg2rad(20)  # elevation mask for AR
-    nav.gnss_t = [uGNSS.GPS, uGNSS.GAL, uGNSS.QZS]
     nav.x[0:3] = pos0
 
     dP = np.diag(nav.P)
@@ -98,25 +98,19 @@ def rtkinit(nav, pos0=np.zeros(3)):
         nav.q[3:6] = nav.sig_qv**2
     else:
         nav.q[0:3] = nav.sig_qp**2
-    # obs index
-    i0 = {gn.uGNSS.GPS: 0, gn.uGNSS.GAL: 0, gn.uGNSS.QZS: 0}
-    i1 = {gn.uGNSS.GPS: 1, gn.uGNSS.GAL: 2, gn.uGNSS.QZS: 1}
-    freq0 = {gn.uGNSS.GPS: nav.freq[0], gn.uGNSS.GAL: nav.freq[0],
-             gn.uGNSS.QZS: nav.freq[0]}
-    freq1 = {gn.uGNSS.GPS: nav.freq[1], gn.uGNSS.GAL: nav.freq[2],
-             gn.uGNSS.QZS: nav.freq[1]}
-    nav.obs_idx = [i0, i1]
-    nav.obs_freq = [freq0, freq1]
+
     nav.cs_sig_idx = {gn.uGNSS.GPS: [sSigGPS.L1C, sSigGPS.L2W],
                       gn.uGNSS.GAL: [sSigGAL.L1X, sSigGAL.L5X],
                       gn.uGNSS.QZS: [sSigQZS.L1C, sSigQZS.L2X]}
 
     nav.fout = None
     nav.logfile = 'log.txt'
-    if nav.loglevel >= 2:
+    if nav.monlevel >= 2:
         nav.fout = open(nav.logfile, 'w')
-        nav.fout.write("# tow\tsys\tprn\tpb1\tpb2\tcb1\tcb2\ttrop\tiono\tantr1\tantr2")
-        nav.fout.write("# \trelatv\twup1\twup2\tCPC1\tCPC2\tPRC1\tPRC2\tdorb\tdclk\n")
+        nav.fout.write(
+            "# tow\tsys\tprn\tpb1\tpb2\tcb1\tcb2\ttrop\tiono\tantr1\tantr2")
+        nav.fout.write(
+            "# \trelatv\twup1\twup2\tCPC1\tCPC2\tPRC1\tPRC2\tdorb\tdclk\n")
 
 
 def udstate(nav, obs, cs):
@@ -150,7 +144,7 @@ def udstate(nav, obs, cs):
             nav.outc[i, f] += 1
             reset = (nav.outc[i, f] > nav.maxout)
             sys_i, _ = gn.sat2prn(sat_)
-            if sys_i not in nav.gnss_t:
+            if sys_i not in obs.sig.keys():
                 continue
             j = IB(sat_, f, nav.na)
             if reset and nav.x[j] != 0.0:
@@ -158,10 +152,9 @@ def udstate(nav, obs, cs):
                 nav.outc[i, f] = 0
         # cycle slip check by LLI
         for i in range(ns):
-            if sys[i] not in nav.gnss_t:
+            if sys[i] not in obs.sig.keys():
                 continue
-            j = nav.obs_idx[f][sys[i]]
-            if obs.lli[i, j] & 1 == 0:
+            if obs.lli[i, f] & 1 == 0:
                 continue
             initx(nav, 0.0, 0.0, IB(sat[i], f, nav.na))
         # reset bias if correction is not available
@@ -174,15 +167,16 @@ def udstate(nav, obs, cs):
         offset = 0
         na = 0
         for i in range(ns):
-            if sys[i] not in nav.gnss_t:
+            if sys[i] not in obs.sig.keys():
                 continue
-            j = nav.obs_idx[f][sys[i]]
-            freq = nav.obs_freq[f][sys[i]]
-            cp = obs.L[i, j]
-            pr = obs.P[i, j]
-            if cp == 0.0 or pr == 0.0 or freq == 0.0:
+
+            lam = obs.sig[sys[i]][gn.uTYP.L][f].wavelength()
+
+            cp = obs.L[i, f]
+            pr = obs.P[i, f]
+            if cp == 0.0 or pr == 0.0 or lam is None:
                 continue
-            bias[i] = cp-pr*freq/gn.rCST.CLIGHT
+            bias[i] = cp-pr/lam
             amb = nav.x[IB(sat[i], f, nav.na)]
             if amb != 0.0:
                 offset += bias[i]-amb
@@ -235,7 +229,7 @@ def zdres(nav, obs, rs, vs, dts, svh, rr, cs):
     for i in range(n):
         sat = obs.sat[i]
         sys, _ = gn.sat2prn(sat)
-        if svh[i] > 0 or sys not in nav.gnss_t or sat in nav.excl_sat:
+        if svh[i] > 0 or sys not in obs.sig.keys() or sat in nav.excl_sat:
             continue
         if sat not in cs.lc[inet].sat_n:
             continue
@@ -253,9 +247,8 @@ def zdres(nav, obs, rs, vs, dts, svh, rr, cs):
         # check for measurement consistency
         flg_m = True
         for f in range(nav.nf):
-            k = nav.obs_idx[f][sys]
-            if obs.P[i, k] == 0.0 or obs.L[i, k] == 0.0 \
-                    or obs.lli[i, k] == 1:
+            if obs.P[i, f] == 0.0 or obs.L[i, f] == 0.0 \
+                    or obs.lli[i, f] == 1:
                 flg_m = False
         if obs.S[i, 0] < nav.cnr_min:
             flg_m = False
@@ -267,14 +260,17 @@ def zdres(nav, obs, rs, vs, dts, svh, rr, cs):
         if el[i] < nav.elmin:
             continue
 
-        freq = np.zeros(nav.nf)
-        lam = np.zeros(nav.nf)
-        iono = np.zeros(nav.nf)
-        for f in range(nav.nf):
-            freq[f] = nav.obs_freq[f][sys]
-            lam[f] = gn.rCST.CLIGHT/freq[f]
-            iono[f] = 40.3e16/(freq[f]*freq[f])*stec[idx_l]
-        iono_ = 40.3e16/(freq[0]*freq[1])*stec[idx_l]
+        sigsPR = obs.sig[sys][gn.uTYP.C]
+        sigsCP = obs.sig[sys][gn.uTYP.L]
+
+        frq = [s.frequency() for s in sigsPR]
+        ionoPR = np.array([40.3e16/(f*f)*stec[idx_l] for f in frq])
+
+        frq = np.array([s.frequency() for s in sigsCP])
+        lam = np.array([s.wavelength() for s in sigsCP])
+        ionoCP = np.array([40.3e16/(f*f)*stec[idx_l] for f in frq])
+
+        iono_ = 40.3e16/(frq[0]*frq[1])*stec[idx_l]
 
         # global/local signal bias
         cbias = np.zeros(nav.nf)
@@ -299,25 +295,29 @@ def zdres(nav, obs, rs, vs, dts, svh, rr, cs):
         # tropospheric delay
         mapfh, mapfw = gn.tropmapf(obs.t, pos, el[i])
         trop = mapfh*trph*r_hs+mapfw*trpw*r_wet
+
         # phase wind-up effect
         nav.phw[sat-1] = windupcorr(obs.t, rs[i, :], vs[i, :], rr_,
                                     nav.phw[sat-1])
         phw = lam*nav.phw[sat-1]
-        antr = antmodel(nav, el[i], nav.nf)
+
+        antrPR = antModelRx(nav, pos, e[i, :], sigsPR)
+        antrCP = antModelRx(nav, pos, e[i, :], sigsCP)
+
         # range correction
-        prc_c = trop+relatv+antr
+        prc_c = trop+relatv
+
         # prc_c += nav.dorb[sat]-nav.dclk[sat]
-        cs.prc[i, :] = prc_c+iono+cbias
-        cs.cpc[i, :] = prc_c-iono+pbias+phw
+        cs.prc[i, :] = prc_c+antrPR+ionoPR+cbias
+        cs.cpc[i, :] = prc_c+antrCP-ionoCP+pbias+phw
         cs.osr[i, :] = [pbias[0], pbias[1], cbias[0], cbias[1],
-                        antr[0], antr[1], phw[0], phw[1],
+                        antrPR[0], antrPR[1], phw[0], phw[1],
                         trop, iono_, relatv, nav.dorb[sat], nav.dclk[sat]]
         r += -_c*dts[i]
 
         for f in range(nf):
-            k = nav.obs_idx[f][sys]
-            y[i, f] = obs.L[i, k]*lam[f]-(r+cs.cpc[i, f])
-            y[i, f+nf] = obs.P[i, k]-(r+cs.prc[i, f])
+            y[i, f] = obs.L[i, f]*lam[f]-(r+cs.cpc[i, f])
+            y[i, f+nf] = obs.P[i, f]-(r+cs.prc[i, f])
 
     return y, e, el
 
@@ -338,7 +338,7 @@ def ppprtkpos(nav, obs, cs):
     rs, vs, dts, svh = satposs(obs, nav, cs)
     # Kalman filter time propagation
     udstate(nav, obs, cs)
-    xa = np.zeros(nav.nx)
+    #xa = np.zeros(nav.nx)
     xp = nav.x.copy()
 
     # non-differential residuals for rover
@@ -354,7 +354,7 @@ def ppprtkpos(nav, obs, cs):
     nav.y = y
     ns = len(sat)
 
-    if nav.loglevel > 1:
+    if nav.monlevel > 1:
         logmon(nav, obs.t, sat, cs, iu)
 
     ny = y.shape[0]
@@ -364,7 +364,7 @@ def ppprtkpos(nav, obs, cs):
         return -1
 
     # DD residual
-    v, H, R = ddres(nav, xp, y, e, sat, el)
+    v, H, R = ddres(nav, obs, xp, y, e, sat, el)
     Pp = nav.P.copy()
     # Kalman filter measurement update
     xp, Pp, _ = kfupdate(xp, Pp, H, v, R)
@@ -379,7 +379,7 @@ def ppprtkpos(nav, obs, cs):
         return -1
 
     # residual for float solution
-    v, H, R = ddres(nav, xp, y, e, sat, el)
+    v, H, R = ddres(nav, obs, xp, y, e, sat, el)
     if valpos(nav, v, R):
         nav.x = xp
         nav.P = Pp
@@ -402,7 +402,7 @@ def ppprtkpos(nav, obs, cs):
         yu, eu, elu = zdres(nav, obs, rs, vs, dts, svh, xa[0:3], cs)
         y = yu[iu, :]
         e = eu[iu, :]
-        v, H, R = ddres(nav, xa, y, e, sat, el)
+        v, H, R = ddres(nav, obs, xa, y, e, sat, el)
         if valpos(nav, v, R):  # R <= Q=H'PH+R  chisq<max_inno[3] (0.5)
             if nav.armode == 3:  # fix and hold
                 holdamb(nav, xa)  # hold fixed ambiguity
diff --git a/src/cssrlib/rinex.py b/src/cssrlib/rinex.py
index b2b254f..acde733 100644
--- a/src/cssrlib/rinex.py
+++ b/src/cssrlib/rinex.py
@@ -8,8 +8,9 @@
 from cssrlib.gnss import prn2sat, char2sys
 from cssrlib.gnss import Eph, Obs
 
+
 class pclk_t:
-    def __init__(self, time = None):
+    def __init__(self, time=None):
         if time is not None:
             self.time = time
         else:
@@ -17,16 +18,19 @@ def __init__(self, time = None):
         self.clk = np.zeros(uGNSS.MAXSAT)
         self.std = np.zeros(uGNSS.MAXSAT)
 
+
 class rnxdec:
     """ class for RINEX decoder """
-    MAXSAT = uGNSS.GPSMAX+uGNSS.GLOMAX+uGNSS.GALMAX+uGNSS.BDSMAX+uGNSS.QZSMAX
 
     def __init__(self):
+
         self.ver = -1.0
         self.fobs = None
-        
-        self.sig_map = {}  # signal code mapping to columns in data section
-        self.sig_tab = {}  # signal selection for internal data structure
+
+        # signal code mapping from RINEX header to columns in data section
+        self.sig_map = {}
+        # signal selection for internal data structure
+        self.sig_tab = {}
         self.nsig = {uTYP.C: 0, uTYP.L: 0, uTYP.D: 0, uTYP.S: 0}
 
         self.pos = np.array([0, 0, 0])
@@ -93,13 +97,14 @@ def decode_nav(self, navfile, nav):
                             nav.ion[1, k] = self.flt(line[5+k*12:5+(k+1)*12])
 
             for line in fnav:
+
                 sys = char2sys(line[0])
-                
+
                 # Skip undesired constellations
-                #                
+                #
                 if sys not in (uGNSS.GPS, uGNSS.GAL, uGNSS.QZS):
                     continue
-                
+
                 prn = int(line[1:3])
                 if sys == uGNSS.QZS:
                     prn += 192
@@ -164,19 +169,20 @@ def decode_nav(self, navfile, nav):
         return nav
 
     def decode_clk(self, clkfile, nav):
-        """decode RINEX Navigation message from file """
+        """decode Clock-RINEX data from file """
         nav.pclk = []
         with open(clkfile, 'rt') as fnav:
             for line in fnav:
+
                 if line[0:2] != 'AS':
                     continue
-                if line[3] not in self.gnss_tbl:
-                    continue
-                sys = self.gnss_tbl[line[3]]
+
+                sys = char2sys(line[3])
                 prn = int(line[4:7])
                 if sys == uGNSS.QZS:
                     prn += 192
                 sat = prn2sat(sys, prn)
+
                 t = self.decode_time(line, 8, 9)
                 if nav.nc <= 0 or abs(timediff(nav.pclk[-1].time, t)) > 1e-9:
                     nav.nc += 1
@@ -215,6 +221,7 @@ def decode_obsh(self, obsfile):
                                      float(line[28:42]),  # North
                                      float(line[0:14])])  # Up
             elif line[60:79] == 'SYS / # / OBS TYPES':
+
                 gns = char2sys(line[0])
                 nsig = int(line[3:6])
 
@@ -236,16 +243,21 @@ def decode_obsh(self, obsfile):
                 self.ts = epoch2time([float(v) for v in line[0:44].split()])
             elif 'TIME OF LAST OBS' in line[60:]:
                 self.te = epoch2time([float(v) for v in line[0:44].split()])
+
         return 0
 
     def decode_obs(self):
-        """decode RINEX Observation message from file """
+        """ decode RINEX Observation message from file """
+
         obs = Obs()
 
         for line in self.fobs:
+
             if line[0] != '>':
                 continue
+
             nsat = int(line[32:35])
+
             year = int(line[2:6])
             month = int(line[7:9])
             day = int(line[10:12])
@@ -253,13 +265,18 @@ def decode_obs(self):
             minute = int(line[16:18])
             sec = float(line[19:29])
             obs.t = epoch2time([year, month, day, hour, minute, sec])
+
+            # Initialize data structures
+            #
             obs.P = np.empty((0, self.nsig[uTYP.C]), dtype=np.float64)
             obs.L = np.empty((0, self.nsig[uTYP.L]), dtype=np.float64)
             obs.S = np.empty((0, self.nsig[uTYP.S]), dtype=np.float64)
             obs.lli = np.empty((0, self.nsig[uTYP.L]), dtype=np.int32)
             obs.sat = np.empty(0, dtype=np.int32)
+            obs.sig = self.sig_tab
 
             for _ in range(nsat):
+
                 line = self.fobs.readline()
                 sys = char2sys(line[0])
 
@@ -278,6 +295,8 @@ def decode_obs(self):
                 prn = int(line[1:3])
                 if sys == uGNSS.QZS:
                     prn += 192
+                elif sys == uGNSS.SBS:
+                    prn += 100
                 sat = prn2sat(sys, prn)
 
                 pr = np.zeros(len(self.getSignals(sys, uTYP.C)),
@@ -297,23 +316,27 @@ def decode_obs(self):
                             sig not in self.sig_tab[sys][sig.typ]:
                         continue
 
-                    j = self.sig_tab[sys][sig.typ].index(sig)
-                    obs_ = line[16*i+4:16*i+17].strip()
-                    
-                    # Replace empty string with zero
-                    # TBD: consider using None-type instead?
+                    # Get string representation of measurement value
                     #
-                    if not obs_:
-                        obs_ = '0.0'
+                    sval = line[16*i+3:16*i+17].strip()
+                    slli = line[16*i+17] if len(line) > 16*i+17 else ''
+
+                    # Convert from string to numerical value
+                    #
+                    val = 0.0 if not sval else float(sval)
+                    lli = 1 if slli == '1' else 0
+
+                    # Signal index in data structure
+                    #
+                    j = self.sig_tab[sys][sig.typ].index(sig)
 
                     if sig.typ == uTYP.C:
-                        pr[j] = float(obs_)
+                        pr[j] = val
                     elif sig.typ == uTYP.L:
-                        cp[j] = float(obs_)
-                        if line[16*i+17] == '1':
-                            ll[j] = 1
+                        cp[j] = val
+                        ll[j] = lli
                     elif sig.typ == uTYP.S:
-                        cn[j] = float(obs_)
+                        cn[j] = val
                     else:
                         continue
 
@@ -329,7 +352,9 @@ def decode_obs(self):
             obs.L = obs.L.reshape(len(obs.sat), self.nsig[uTYP.L])
             obs.S = obs.S.reshape(len(obs.sat), self.nsig[uTYP.S])
             obs.lli = obs.lli.reshape(len(obs.sat), self.nsig[uTYP.L])
+
             break
+
         return obs
 
 
diff --git a/src/cssrlib/rtk.py b/src/cssrlib/rtk.py
index 8018606..41d88ad 100644
--- a/src/cssrlib/rtk.py
+++ b/src/cssrlib/rtk.py
@@ -4,13 +4,14 @@
 """
 
 import numpy as np
+
 import cssrlib.gnss as gn
 from cssrlib.ephemeris import satposs
+from cssrlib.gnss import uTYP
+from cssrlib.peph import antModelRx
 from cssrlib.ppp import tidedisp
 from cssrlib.mlambda import mlambda
 
-from cssrlib.peph import searchpcv, antModelRx, antModelTx
-
 VAR_HOLDAMB = 0.001
 
 
@@ -22,7 +23,7 @@ def rtkinit(nav, pos0=np.zeros(3)):
 
     nav.na = 3 if nav.pmode == 0 else 6
     nav.nq = 3 if nav.pmode == 0 else 6
-    nav.ratio = 0
+
     nav.thresar = [2]
     nav.nx = nav.na+gn.uGNSS.MAXSAT*nav.nf
     nav.x = np.zeros(nav.nx)
@@ -30,9 +31,8 @@ def rtkinit(nav, pos0=np.zeros(3)):
     nav.xa = np.zeros(nav.na)
     nav.Pa = np.zeros((nav.na, nav.na))
     nav.el = np.zeros(gn.uGNSS.MAXSAT)
-    nav.nfix = nav.neb = 0
 
-    # parameter for RTK    
+    # parameter for RTK
     nav.eratio = [50, 50]
     nav.err = [0, 0.01, 0.005]/np.sqrt(2)
     nav.sig_p0 = 30.0
@@ -40,10 +40,9 @@ def rtkinit(nav, pos0=np.zeros(3)):
     nav.sig_n0 = 30.0
     nav.sig_qp = 0.01
     nav.sig_qv = 1.0
-    
+
     nav.armode = 1  # 1:contunous,2:instantaneous,3:fix-and-hold
     nav.elmaskar = np.deg2rad(20)  # elevation mask for AR
-    nav.gnss_t = [gn.uGNSS.GPS, gn.uGNSS.GAL, gn.uGNSS.QZS]
     nav.x[0:3] = pos0
     nav.x[3:6] = 0.0
 
@@ -58,16 +57,6 @@ def rtkinit(nav, pos0=np.zeros(3)):
     else:
         nav.q[0:3] = nav.sig_qp**2
 
-    # obs index
-    i0 = {gn.uGNSS.GPS: 0, gn.uGNSS.GAL: 0, gn.uGNSS.QZS: 0}
-    i1 = {gn.uGNSS.GPS: 1, gn.uGNSS.GAL: 2, gn.uGNSS.QZS: 1}
-    freq0 = {gn.uGNSS.GPS: nav.freq[0], gn.uGNSS.GAL: nav.freq[0],
-             gn.uGNSS.QZS: nav.freq[0]}
-    freq1 = {gn.uGNSS.GPS: nav.freq[1], gn.uGNSS.GAL: nav.freq[2],
-             gn.uGNSS.QZS: nav.freq[1]}
-    nav.obs_idx = [i0, i1]
-    nav.obs_freq = [freq0, freq1]
-
 
 def zdres(nav, obs, rs, dts, svh, rr, rtype=1):
     """ non-differencial residual """
@@ -85,7 +74,7 @@ def zdres(nav, obs, rs, dts, svh, rr, rtype=1):
     pos = gn.ecef2pos(rr_)
     for i in range(n):
         sys, _ = gn.sat2prn(obs.sat[i])
-        if svh[i] > 0 or sys not in nav.gnss_t or obs.sat[i] in nav.excl_sat:
+        if svh[i] > 0 or sys not in obs.sig.keys() or obs.sat[i] in nav.excl_sat:
             continue
         r, e[i, :] = gn.geodist(rs[i, :], rr_)
         _, el[i] = gn.satazel(pos, e[i, :])
@@ -96,24 +85,23 @@ def zdres(nav, obs, rs, dts, svh, rr, rtype=1):
         mapfh, _ = gn.tropmapf(obs.t, pos, el[i])
         r += mapfh*zhd
 
-        #dant = gn.antmodel(nav, el[i], nav.nf, rtype)
-        if rtype == 1:
-            sigs = nav.sig_tab[sys][gn.uTYP.C]
-        else:
-            sigs = nav.sig_tab_b[sys][gn.uTYP.C]
-        
-        dant = antModelRx(nav, rr_, e[i, :], sigs, rtype)
-        
+        sigPR = obs.sig[sys][gn.uTYP.C]
+        sigCP = obs.sig[sys][gn.uTYP.L]
+
+        dantPR = antModelRx(nav, pos, e[i, :], sigPR, rtype)
+        dantCP = antModelRx(nav, pos, e[i, :], sigCP, rtype)
+
+        lam = [s.wavelength() for s in sigCP]
+
         for f in range(nf):
-            j = nav.obs_idx[f][sys]
             if obs.L[i, f] == 0.0:
                 y[i, f] = 0.0
             else:
-                y[i, f] = obs.L[i, f]*_c/nav.freq[j]-r-dant[f]
+                y[i, f] = obs.L[i, f]*lam[f]-r-dantCP[f]
             if obs.P[i, f] == 0.0:
                 y[i, f+nf] = 0.0
             else:
-                y[i, f+nf] = obs.P[i, f]-r-dant[f]
+                y[i, f+nf] = obs.P[i, f]-r-dantPR[f]
     return y, e, el
 
 
@@ -158,13 +146,13 @@ def varerr(nav, el, f):
     return 2.0*(a**2+(b/s_el)**2)
 
 
-def ddres(nav, x, y, e, sat, el):
+def ddres(nav, obs, x, y, e, sat, el):
     """ DD phase/code residual """
     _c = gn.rCST.CLIGHT
     nf = nav.nf
     ns = len(el)
     mode = 1 if len(y) == ns else 0  # 0:DD,1:SD
-    nb = np.zeros(2*len(nav.gnss_t)*nf, dtype=int)
+    nb = np.zeros(2*len(obs.sig.keys())*nf, dtype=int)
     Ri = np.zeros(ns*nf*2)
     Rj = np.zeros(ns*nf*2)
 
@@ -172,13 +160,9 @@ def ddres(nav, x, y, e, sat, el):
     b = 0
     H = np.zeros((ns*nf*2, nav.nx))
     v = np.zeros(ns*nf*2)
-    idx_f = [0, 1]
-    for sys in nav.gnss_t:
-        for f in range(nf):
-            idx_f[f] = nav.obs_idx[f][sys]
+
+    for sys in obs.sig.keys():
         for f in range(0, nf*2):
-            if f < nf:
-                freq = nav.freq[idx_f[f]]
             # reference satellite
             idx = sysidx(sat, sys)
             if len(idx) > 0:
@@ -199,7 +183,7 @@ def ddres(nav, x, y, e, sat, el):
                 if f < nf:  # carrier
                     idx_i = IB(sat[i], f, nav.na)
                     idx_j = IB(sat[j], f, nav.na)
-                    lami = _c/freq
+                    lami = obs.sig[sys][gn.uTYP.L][f].wavelength()
                     v[nv] -= lami*(x[idx_i]-x[idx_j])
                     H[nv, idx_i] = lami
                     H[nv, idx_j] = -lami
@@ -221,7 +205,7 @@ def ddres(nav, x, y, e, sat, el):
 
 
 def valpos(nav, v, R, thres=4.0):
-    """ post-file residual test """
+    """ post-fit residual test """
     nv = len(v)
     fact = thres**2
     for i in range(nv):
@@ -245,7 +229,7 @@ def ddidx(nav, sat):
             for i in range(k, k+n):
                 sat_i = i-k+1
                 sys, _ = gn.sat2prn(sat_i)
-                if (sys != m) or sys not in nav.gnss_t:
+                if (sys != m):  # or sys not in nav.gnss_t:
                     continue
                 if sat_i not in sat or nav.x[i] == 0.0 \
                    or nav.vsat[sat_i-1, f] == 0:
@@ -258,7 +242,7 @@ def ddidx(nav, sat):
             for j in range(k, k+n):
                 sat_j = j-k+1
                 sys, _ = gn.sat2prn(sat_j)
-                if (sys != m) or sys not in nav.gnss_t:
+                if (sys != m):  # or sys not in nav.gnss_t:
                     continue
                 if i == j or sat_j not in sat or nav.x[j] == 0.0 \
                    or nav.vsat[sat_j-1, f] == 0:
@@ -284,7 +268,7 @@ def restamb(nav, bias, nb):
             index = []
             for i in range(gn.uGNSS.MAXSAT):
                 sys, _ = gn.sat2prn(i+1)
-                if sys != m or (sys not in nav.gnss_t) or nav.fix[i, f] != 2:
+                if sys != m or nav.fix[i, f] != 2:
                     continue
                 index.append(IB(i+1, f, nav.na))
                 n += 1
@@ -380,7 +364,7 @@ def udstate(nav, obs, obsb, iu, ir):
             nav.outc[i, f] += 1
             reset = (nav.outc[i, f] > nav.maxout)
             sys_i, _ = gn.sat2prn(i+1)
-            if sys_i not in nav.gnss_t:
+            if sys_i not in obs.sig.keys():
                 continue
             j = IB(i+1, f, nav.na)
             if reset and nav.x[j] != 0.0:
@@ -388,9 +372,8 @@ def udstate(nav, obs, obsb, iu, ir):
                 nav.outc[i, f] = 0
         # cycle slip check by LLI
         for i in range(ns):
-            if sys[i] not in nav.gnss_t:
+            if sys[i] not in obs.sig.keys():
                 continue
-            j = nav.obs_idx[f][sys[i]]
             if obsb.lli[ir[i], f] & 1 == 0 and obs.lli[iu[i], f] & 1 == 0:
                 continue
             initx(nav, 0.0, 0.0, IB(sat[i], f, nav.na))
@@ -399,13 +382,12 @@ def udstate(nav, obs, obsb, iu, ir):
         offset = 0
         na = 0
         for i in range(ns):
-            if sys[i] not in nav.gnss_t:
+            if sys[i] not in obs.sig.keys():
                 continue
-            j = nav.obs_idx[f][sys[i]]
-            freq = nav.obs_freq[f][sys[i]]
+            lam = obs.sig[sys[i]][uTYP.C][f].wavelength()
             cp = obs.L[iu[i], f]-obsb.L[ir[i], f]
             pr = obs.P[iu[i], f]-obsb.P[ir[i], f]
-            bias[i] = cp-pr*freq/gn.rCST.CLIGHT
+            bias[i] = cp-pr/lam
             amb = nav.x[IB(sat[i], f, nav.na)]
             if amb != 0.0:
                 offset += bias[i]-amb
@@ -429,10 +411,7 @@ def selsat(nav, obs, obsb, elb):
     """ select common satellite between rover and base station """
     # exclude satellite with missing observation and cycle slip for rover
     idx_u = []
-    for k, sat in enumerate(obs.sat):
-        sys, _ = gn.sat2prn(sat)
-        j0 = nav.obs_idx[0][sys]
-        j1 = nav.obs_idx[1][sys]
+    for k, _ in enumerate(obs.sat):
         if obs.P[k, 0] == 0.0 or obs.P[k, 1] == 0.0 or \
            obs.L[k, 0] == 0.0 or obs.L[k, 1] == 0.0 or \
            obs.lli[k, 0] > 0 or obs.lli[k, 1] > 0:
@@ -441,10 +420,7 @@ def selsat(nav, obs, obsb, elb):
 
     # exclude satellite with missing observation and cycle slip for base
     idx_r = []
-    for k, sat in enumerate(obsb.sat):
-        sys, _ = gn.sat2prn(sat)
-        j0 = nav.obs_idx[0][sys]
-        j1 = nav.obs_idx[1][sys]
+    for k, _ in enumerate(obsb.sat):
         if obsb.P[k, 0] == 0.0 or obsb.P[k, 1] == 0.0 or \
            obsb.L[k, 0] == 0.0 or obsb.L[k, 1] == 0.0 or \
            obsb.lli[k, 0] > 0 or obsb.lli[k, 1] > 0 or \
@@ -514,7 +490,7 @@ def relpos(nav, obs, obsb):
     # Kalman filter time propagation
     udstate(nav, obs, obsb, iu, ir)
 
-    xa = np.zeros(nav.nx)
+    #xa = np.zeros(nav.nx)
     xp = nav.x
 
     # non-differencial residual for rover
@@ -526,7 +502,7 @@ def relpos(nav, obs, obsb):
     sat = obs.sat[iu]
     nav.el[sat-1] = el
     # DD residual
-    v, H, R = ddres(nav, xp, y, e, sat, el)
+    v, H, R = ddres(nav, obs, xp, y, e, sat, el)
     Pp = nav.P
 
     # Kalman filter measurement update
@@ -537,7 +513,7 @@ def relpos(nav, obs, obsb):
     y[:ns, :] = yu[iu, :]
     e[:ns, :] = eu[iu, :]
     # residual for float solution
-    v, H, R = ddres(nav, xp, y, e, sat, el)
+    v, H, R = ddres(nav, obs, xp, y, e, sat, el)
     if valpos(nav, v, R):
         nav.x = xp
         nav.P = Pp
@@ -550,7 +526,7 @@ def relpos(nav, obs, obsb):
         yu, eu, _ = zdres(nav, obs, rs, dts, svh, xa[0:3])
         y[:ns, :] = yu[iu, :]
         e[:ns, :] = eu[iu, :]
-        v, H, R = ddres(nav, xa, y, e, sat, el)
+        v, H, R = ddres(nav, obs, xa, y, e, sat, el)
         if valpos(nav, v, R):
             if nav.armode == 3:
                 holdamb(nav, xa)


Test patch:
diff --git a/src/cssrlib/test/test_gnss.py b/src/cssrlib/test/test_gnss.py
index ccd44d1..dc79674 100644
--- a/src/cssrlib/test/test_gnss.py
+++ b/src/cssrlib/test/test_gnss.py
@@ -1,74 +1,57 @@
-from cssrlib.gnss import uGNSS, uSIG, uTYP, rSigRnx, sys2char
-
-signals = {}
-signals.update({uGNSS.GPS: [uSIG.L1C, uSIG.L2W, uSIG.L5Q]})
-signals.update({uGNSS.GLO: [uSIG.L1P, uSIG.L2C, uSIG.L3X, uSIG.L4A, uSIG.L5B]})
-signals.update({uGNSS.GAL: [uSIG.L1C, uSIG.L5Q, uSIG.L6C, uSIG.L7Q, uSIG.L8X]})
-signals.update({uGNSS.BDS: [uSIG.L1C, uSIG.L2I, uSIG.L5P, uSIG.L6C, uSIG.L7I,
-                            uSIG.L8Q]})
-signals.update({uGNSS.QZS: [uSIG.L1C, uSIG.L2L, uSIG.L5Q, uSIG.L6S]})
-signals.update({uGNSS.SBS: [uSIG.L1C, uSIG.L5Q]})
-signals.update({uGNSS.IRN: [uSIG.L5A, uSIG.L9B]})
-
-for sys, sigs in signals.items():
-    for sig in sigs:
-        rnxSig = rSigRnx(sys, uTYP.C, sig)
-        print("{} {} {:9.4f} MHz  {:5.2f} cm".format(
-              sys2char(sys), rnxSig.str(),
-              rnxSig.frequency(k=-6)*1e-6,
-              rnxSig.wavelength(k=-6)*1e2))
-    print()
-
-print()
-
-print("String representation and modification functions")
-print()
-
-sig = rSigRnx(uGNSS.GPS, 'L1X')
-print("<{:s}>".format(sig.__repr__()))
-
-sig = sig.toTyp(uTYP.C)
-print("<{:s}>".format(sig.__repr__()))
-
-sig = sig.toAtt()
-print("<{:s}>".format(sig.__repr__()))
-
-sig = sig.toTyp(uTYP.D).toAtt('C')
-print("<{:s}>".format(sig.__repr__()))
-
-print()
-
-print("Test constructors")
-print()
-
-try:
-    sig = rSigRnx()
-    print(sys2char(sig.sys), sig.str())
-except ValueError:
-    print('Value Error')
-    
-try:
-    sig = rSigRnx("EC5Q")
-    print(sys2char(sig.sys), sig.str())
-except ValueError:
-    print('Value Error')
-    
-try:
-    sig = rSigRnx(uGNSS.GPS, "D1X")
-    print(sys2char(sig.sys), sig.str())
-except ValueError:
-    print('Value Error')
-    
-try:
-    sig = rSigRnx(uGNSS.IRN, uTYP.S, uSIG.L1X)
-    print(sys2char(sig.sys), sig.str())
-except ValueError:
-    print('Value Error')
-    
-try:
-     sig = rSigRnx(uGNSS.GPS, "D1X", "what??")
-     print(sys2char(sig.sys), sig.str())
-except ValueError:
-    print('Value Error')
-
-
+from cssrlib.gnss import uGNSS, uSIG, uTYP, rSigRnx, sys2char
+
+signals = {}
+signals.update({uGNSS.GPS: [uSIG.L1C, uSIG.L2W, uSIG.L5Q]})
+signals.update({uGNSS.GLO: [uSIG.L1P, uSIG.L2C, uSIG.L3X, uSIG.L4A, uSIG.L5B]})
+signals.update({uGNSS.GAL: [uSIG.L1C, uSIG.L5Q, uSIG.L6C, uSIG.L7Q, uSIG.L8X]})
+signals.update({uGNSS.BDS: [uSIG.L1C, uSIG.L2I, uSIG.L5P, uSIG.L6C, uSIG.L7I,
+                            uSIG.L8Q]})
+signals.update({uGNSS.QZS: [uSIG.L1C, uSIG.L2L, uSIG.L5Q, uSIG.L6S]})
+signals.update({uGNSS.SBS: [uSIG.L1C, uSIG.L5Q]})
+signals.update({uGNSS.IRN: [uSIG.L5A, uSIG.L9B]})
+
+for sys, sigs in signals.items():
+    for sig in sigs:
+        rnxSig = rSigRnx(sys, uTYP.C, sig)
+        print("{} {} {:9.4f} MHz  {:5.2f} cm".format(
+              sys2char(sys), rnxSig.str(),
+              rnxSig.frequency(k=-6)*1e-6,
+              rnxSig.wavelength(k=-6)*1e2))
+    print()
+
+print()
+
+print("String representation and modification functions")
+print()
+
+sig = rSigRnx(uGNSS.GPS, 'L1X')
+print("<{:s}>".format(sig.__repr__()))
+
+sig = sig.toTyp(uTYP.C)
+print("<{:s}>".format(sig.__repr__()))
+
+sig = sig.toAtt()
+print("<{:s}>".format(sig.__repr__()))
+
+sig = sig.toTyp(uTYP.D).toAtt('C')
+print("<{:s}>".format(sig.__repr__()))
+
+print()
+
+print("Test constructors")
+print()
+
+sig = rSigRnx()
+print(sys2char(sig.sys), sig.str())
+
+sig = rSigRnx("EC5Q")
+print(sys2char(sig.sys), sig.str())
+
+sig = rSigRnx(uGNSS.GPS, "D1X")
+print(sys2char(sig.sys), sig.str())
+
+sig = rSigRnx(uGNSS.IRN, uTYP.S, uSIG.L1X)
+print(sys2char(sig.sys), sig.str())
+
+sig = rSigRnx(uGNSS.GPS, "D1X", "what??")
+print(sys2char(sig.sys), sig.str())
diff --git a/src/cssrlib/test/test_peph.py b/src/cssrlib/test/test_peph.py
index 1ddac09..6d598a1 100644
--- a/src/cssrlib/test/test_peph.py
+++ b/src/cssrlib/test/test_peph.py
@@ -16,13 +16,9 @@
 
 bdir = expanduser('~/GNSS_DAT/')
 atxfile = bdir+"IGS/ANTEX/igs14.atx"
-#orbfile = bdir+"COD0IGSRAP/2021/COD0IGSRAP_20210780000_01D_15M_ORB.SP3"
-#clkfile = bdir+"COD0IGSRAP/2021/COD0IGSRAP_20210780000_01D_30S_CLK.CLK"
-#dcbfile = bdir+"COD0IGSRAP/2021/COD0IGSRAP_20210780000_01D_01D_OSB.BIA"
-
-orbfile = bdir+"COD0IGSRAP/2021/COD0MGXFIN_20210780000_01D_05M_ORB.SP3"
-clkfile = bdir+"COD0IGSRAP/2021/COD0MGXFIN_20210780000_01D_30S_CLK.CLK"
-dcbfile = bdir+"COD0IGSRAP/2021/COD0MGXFIN_20210780000_01D_01D_OSB.BIA"
+orbfile = bdir+"COD0IGSRAP/2021/COD0IGSRAP_20210780000_01D_15M_ORB.SP3"
+clkfile = bdir+"COD0IGSRAP/2021/COD0IGSRAP_20210780000_01D_30S_CLK.CLK"
+dcbfile = bdir+"COD0IGSRAP/2021/COD0IGSRAP_20210780000_01D_01D_OSB.BIA"
 
 time = epoch2time([2021, 3, 19, 12, 0, 0])
 sat = id2sat("G01")
@@ -174,7 +170,7 @@
 
     print()
 
-if True:
+if False:
 
     print("Test Bias-SINEX module")
     print()
@@ -186,16 +182,16 @@
     sig = rSigRnx("GC1W")
 
     bias, std, = bd.getosb(sat, time, sig)
-    #assert bias == 7.6934
-    #assert std == 0.0
+    assert bias == 7.6934
+    assert std == 0.0
 
     print("{:s} {:s} {:8.5f} {:6.4f}"
           .format(sat2id(sat), sig.str(), bias, std))
 
     sig = rSigRnx("GL1W")
     bias, std, = bd.getosb(sat, time, sig)
-    #assert bias == 0.00038
-    #assert std == 0.0
+    assert bias == 0.00038
+    assert std == 0.0
 
     print("{:s} {:s} {:8.5f} {:6.4f}"
           .format(sat2id(sat), sig.str(), bias, std))
diff --git a/src/cssrlib/test/test_rnx.py b/src/cssrlib/test/test_rnx.py
index f7774a2..0b40fb9 100644
--- a/src/cssrlib/test/test_rnx.py
+++ b/src/cssrlib/test/test_rnx.py
@@ -2,44 +2,52 @@
  test of RINEX decoder
 """
 
+from datetime import datetime
+
 from cssrlib.rinex import rnxdec
 from cssrlib.gnss import uTYP, rSigRnx
 from cssrlib.gnss import sat2id, sat2prn
 
-from os.path import expanduser
-
-obsfile = '~/GNSS_OBS/IGS/DAILY/2021/078/CHOF00JPN_S_20210780000_01D_30S_MO.rnx'
-#obsfile = '../data/SEPT078M.21O'
+obsfile = '../data/SEPT078M1.21O'
 
 sigs = [rSigRnx("GC1C"), rSigRnx("EC1X"),
         rSigRnx("GC2W"), rSigRnx("EC5X"),
         rSigRnx("GL1C"), rSigRnx("EL1X"),
-        rSigRnx("GL2W"), rSigRnx("EL5X")]
+        rSigRnx("GL2W"), rSigRnx("EL5X"),
+        rSigRnx("GS1C"), rSigRnx("ES1X"),
+        rSigRnx("GS2W"), rSigRnx("ES5X")]
 
 dec = rnxdec()
 dec.setSignals(sigs)
 
-nep = 1
-if dec.decode_obsh(expanduser(obsfile)) >= 0:
+nep = 2
+if dec.decode_obsh(obsfile) >= 0:
 
     for ne in range(nep):
 
         obs = dec.decode_obs()
 
+        print("{:%Y-%m-%d %T}".format(datetime.utcfromtimestamp(obs.t.time)))
+
         for i, sat in enumerate(obs.sat):
 
             txt = "{} ".format(sat2id(sat))
 
             sys, _ = sat2prn(sat)
-            sigs = dec.getSignals(sys, uTYP.C)
+            sigs = obs.sig[sys][uTYP.C]
             for j, sig in enumerate(sigs):
                 txt += "{} {:13.3f}  ".format(sig.str(), obs.P[i, j])
-            sigs = dec.getSignals(sys, uTYP.L)
+
+            sigs = obs.sig[sys][uTYP.L]
             for j, sig in enumerate(sigs):
                 txt += "{} {:13.3f}  ".format(sig.str(), obs.L[i, j])
 
+            sigs = obs.sig[sys][uTYP.S]
+            for j, sig in enumerate(sigs):
+                txt += "{} {:7.3f}  ".format(sig.str(), obs.S[i, j])
+
             print(txt)
 
-    print()
+        print()
 
     dec.fobs.close()
diff --git a/src/cssrlib/test/test_rtk.py b/src/cssrlib/test/test_rtk.py
index 02b4399..11831b0 100644
--- a/src/cssrlib/test/test_rtk.py
+++ b/src/cssrlib/test/test_rtk.py
@@ -4,31 +4,40 @@
 
 import matplotlib.pyplot as plt
 import numpy as np
+from os.path import expanduser
+import sys
+
 import cssrlib.rinex as rn
 import cssrlib.gnss as gn
 from cssrlib.rtk import rtkinit, relpos
-from cssrlib.peph import searchpcv, atxdec
-from os.path import expanduser
+from cssrlib.gnss import rSigRnx
+from cssrlib.peph import atxdec, searchpcv
 
 bdir = '../data/'
 navfile = bdir+'SEPT078M.21P'
 obsfile = bdir+'SEPT078M1.21O'
 basefile = bdir+'3034078M1.21O'
 
+bdir = expanduser('~/GNSS_DAT/')
+atxfile = bdir+"IGS/ANTEX/igs14.atx"
+
 xyz_ref = [-3962108.673,   3381309.574,   3668678.638]
 pos_ref = gn.ecef2pos(xyz_ref)
 
-sigs_b = [rn.rSigRnx("GC1C"), rn.rSigRnx("EC1X"),
-          rn.rSigRnx("GC2W"), rn.rSigRnx("EC5X"),
-          rn.rSigRnx("GL1C"), rn.rSigRnx("EL1X"),
-          rn.rSigRnx("GL2W"), rn.rSigRnx("EL5X")]
+# Define signals to be processed
+#
+sigs = [rSigRnx("GC1C"), rSigRnx("GC2W"),
+        rSigRnx("EC1C"), rSigRnx("EC5Q"),
+        rSigRnx("GL1C"), rSigRnx("GL2W"),
+        rSigRnx("EL1C"), rSigRnx("EL5Q")]
 
-sigs = [rn.rSigRnx("GC1C"), rn.rSigRnx("EC1C"),
-        rn.rSigRnx("GC2W"), rn.rSigRnx("EC5Q"),
-        rn.rSigRnx("GL1C"), rn.rSigRnx("EL1C"),
-        rn.rSigRnx("GL2W"), rn.rSigRnx("EL5Q")]
+sigsb = [rSigRnx("GC1C"), rSigRnx("GC2W"),
+         rSigRnx("EC1X"), rSigRnx("EC5X"),
+         rSigRnx("GL1C"), rSigRnx("GL2W"),
+         rSigRnx("EL1X"), rSigRnx("EL5X")]
 
 # rover
+#
 dec = rn.rnxdec()
 dec.setSignals(sigs)
 
@@ -36,36 +45,38 @@
 dec.decode_nav(navfile, nav)
 
 # base
+#
 decb = rn.rnxdec()
-decb.setSignals(sigs_b)
+decb.setSignals(sigsb)
 
 decb.decode_obsh(basefile)
 dec.decode_obsh(obsfile)
 
-nep = 60
-#nep = 45
+#
+# Set rover and base antenna
+
+dec.ant = "{:16s}{:4s}".format("JAVRINGANT_DM", "SCIS")
+decb.ant = "{:16s}{:4s}".format("TRM59800.80", "NONE")
 
-time = gn.epoch2time([2021, 3, 19, 12, 0, 0])
+atx = atxdec()
+atx.readpcv(atxfile)
 
-if True:
+# Set antenna PCO/PCV data
+#
+nav.sat_pcv = atx.pcvs
 
-    bdir = expanduser('~/GNSS_DAT/')
-    atxfile = bdir+"IGS/ANTEX/igs14.atx"
-    
-    atx = atxdec()
-    atx.readpcv(atxfile)
-    
-    antr = "{:16s}{:4s}".format("JAVRINGANT_DM", "SCIS")
-    antb = "{:16s}{:4s}".format("TRM59800.80", "NONE")
+nav.rcv_ant = searchpcv(atx.pcvr, dec.ant,  dec.ts)
+nav.rcv_ant_b = searchpcv(atx.pcvr, decb.ant, decb.ts)
 
-    # Store PCV/PCO information for satellites and receiver
-    #
-    nav.sat_ant = atx.pcvs
-    nav.rcv_ant = searchpcv(atx.pcvr, antr, time)
-    nav.rcv_ant_b = searchpcv(atx.pcvr, antb, time)
-    nav.sig_tab = dec.sig_tab
-    nav.sig_tab_b = decb.sig_tab
+if nav.rcv_ant is None:
+    print("ERROR: no rover antenna found!")
+    sys.exit(1)
 
+if nav.rcv_ant_b is None:
+    print("ERROR: no base antenna found!")
+    sys.exit(1)
+
+nep = 60
 
 # GSI 3034 fujisawa
 nav.rb = [-3959400.631, 3385704.533, 3667523.111]