Repository: hirokawa/cssrlib (language: python)

Issue:
Merge latest development changes

Revised modules using new RINEXv3 signal class



Discussion / hints:

Fix patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
new file mode 100644
index 0000000..dccd3b2
--- /dev/null
+++ b/CHANGELOG.md
@@ -0,0 +1,74 @@
+# Changelog
+All notable changes to this project will be documented in this file.
+
+The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
+
+
+# [1.0.0] 2024-01-01
+
+### Added
+New integrated class structure (PPPOS) for PPP/PPP-RTK/RTK processing
+
+Support for RTCM3 (Galileo HAS IDD)
+
+Decoder for Septentrio receiver
+
+Parser for RTCM3, L5 SBAS
+
+Experimental support for PPP via SouthPAN (PVS)
+
+New solid Earth tides model using PySolid (2010 IERS Conventions) 
+
+### Changed
+
+Improved documentation
+
+Sign of SSR satellite signal code/phase bias align with RTCM 3 convention
+
+PPPIGS was integrated into PPPOS and removed
+
+### Fixed
+
+Link for cssrlib-data
+
+### Deprecated
+### Removed
+Function based PPP/PPP-RTK/RTK processing
+### Security
+
+# [0.8.0] 2023-09-09
+
+### Added
+New signal structure
+
+Support for open PPP services: Galileo HAS (SIS), BDS PPP, QZSS MADOCA-PPP
+
+Support for IGS (SP3+BIAS)
+
+Parser for SP3, ANTEX, BIAS files
+
+Support for PPP-AR
+
+Jupyter notebook with examples
+
+### Changed
+
+Improved documentation
+
+Added link for Google Colab
+
+### Fixed
+### Deprecated
+### Removed
+### Security
+
+# [0.3.0] 2022-03-01
+
+### Added
+Initial version for PPP-RTK (QZSS CLAS) and RTK
+
+### Changed
+### Fixed
+### Deprecated
+### Removed
+### Security
diff --git a/README.md b/README.md
index d0e2640..7f0e89f 100644
--- a/README.md
+++ b/README.md
@@ -15,7 +15,6 @@ It supports the following open format:
 ************************************************
 Click this button for a quick demo in Google Colab
 
-
 <a target="_blank" href="https://colab.research.google.com/github/hirokawa/cssrlib/blob/main/tutorials/cssrlib.ipynb">
   <img src="https://colab.research.google.com/assets/colab-badge.svg" alt="Open In Colab"/>
 </a>
@@ -27,14 +26,14 @@ Prerequisites
 Additional python packages are required as prerequisites and can be installed via the following commands
 
 ```
-pip install bitstruct galois crccheck
+pip install bitstruct galois crccheck pysolid
 pip install notebook numpy matplotlib
 ```
 
-Optionally, on linux, users can install the `cartopy` package
+Optionally, on linux, users can install the `cartopy` package and the `PySolid` package
 
 ```
-pip install cartopy
+pip install cartopy pysolid
 ```
 
 If the installation of `cartopy` fails, try installing `libgeos++-dev` first.
@@ -75,9 +74,14 @@ Run RTK sample.
  python test_rtk.py
 ```
 
-Other samples are also available in a separate repository [`cssrlib-data`](https://github.com/hirokawa/cssrlib-data) including :
+Other samples with dataset are also available in a separate repository [`cssrlib-data`](https://github.com/hirokawa/cssrlib-data) including :
+
+- Galileo-HAS PPP (SIS/IDD) processing example
+- BDS-PPP processing example
+- QZSS MADOCA-PPP processing example
+- QZSS CLAS PPP-RTK processing example
+- IGS (SP3+BIAS) PPP-AR processing example
+- (experimental) PPP via SouthPAN (PVS) processing example
+- Binary message parser for Septentrio receiver (Mosaic/PolaRX5)
 
-- Galileo-HAS decoder
-- BDS-PPP decoder
-- Jupyter-Notebook including samples for RTK/PPP-RTK.
 
diff --git a/setup.cfg b/setup.cfg
index a449dfc..9a7bd04 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -1,9 +1,8 @@
 [metadata]
 name = cssrlib
-version = 0.8.0
 author = Rui Hirokawa
 author_email = rui.hirokawa@gmail.com
-description = Compact SSR Library for PPP/PPP-RTK
+description = Python Toolkit for PPP/PPP-RTK positioning
 long_description = file: README.md
 long_description_content_type = text/markdown
 url = https://github.com/hirokawa/cssrlib
diff --git a/src/cssrlib/cssr_bds.py b/src/cssrlib/cssr_bds.py
index 630f929..a2d8856 100644
--- a/src/cssrlib/cssr_bds.py
+++ b/src/cssrlib/cssr_bds.py
@@ -156,24 +156,24 @@ def decode_cssr_mask(self, msg, i):
             self.add_gnss(mask_glo, 37, sGNSS.GLO)
 
             inet = 0
-            self.lc[inet].dclk = np.ones(self.nsat_n)*np.nan
-            self.lc[inet].dorb = np.ones((self.nsat_n, 3))*np.nan
-            self.lc[inet].iode = np.zeros(self.nsat_n, dtype=int)
-            self.lc[inet].iodc = np.zeros(self.nsat_n, dtype=int)
-            self.lc[inet].iodc_c = np.zeros(self.nsat_n, dtype=int)
-            self.lc[inet].cbias = np.ones((self.nsat_n, self.nsig_max))*np.nan
+            self.lc[inet].dclk = {}
+            self.lc[inet].dorb = {}
+            self.lc[inet].iode = {}
+            self.lc[inet].iodc = {}
+            self.lc[inet].iodc_c = {}
+            self.lc[inet].cbias = {}
             self.nsig_n = np.ones(self.nsat_n, dtype=int)*self.nsig_max
-            self.sig_n = -1*np.ones((self.nsat_n, self.nsig_max), dtype=int)
-            self.ura = np.zeros(self.nsat_n)
+            self.sig_n = {}
+            self.ura = {}
 
             # fallback for inconsistent clock update
-            self.lc[inet].dclk_p = np.ones(self.nsat_n)*np.nan
-            self.lc[inet].iodc_c_p = np.zeros(self.nsat_n, dtype=int)
+            self.lc[inet].dclk_p = {}
+            self.lc[inet].iodc_c_p = {}
 
         self.iodssr = head['iodssr']
 
         self.lc[0].cstat |= (1 << sCType.MASK)
-        self.lc[0].t0[sCType.MASK] = self.time
+        # self.lc[0].t0[0][sCType.MASK] = self.time
         return i
 
     def decode_cssr_orb_sat(self, msg, i, inet, sat_n):
@@ -186,17 +186,24 @@ def decode_cssr_orb_sat(self, msg, i, inet, sat_n):
         sat = prn2sat(sys, prn)
         if sat not in sat_n:
             return i
-        idx = np.where(sat == sat_n)[0][0]
 
         if (sys == uGNSS.GPS) or (sys == uGNSS.BDS):
             iodn = iodn & 0xff  # IODC -> IODE
 
-        self.lc[inet].iode[idx] = iodn
-        self.lc[inet].iodc[idx] = iodc
-        self.lc[inet].dorb[idx, 0] = self.sval(dx, 15, self.dorb_scl[0])
-        self.lc[inet].dorb[idx, 1] = self.sval(dy, 13, self.dorb_scl[1])
-        self.lc[inet].dorb[idx, 2] = self.sval(dz, 13, self.dorb_scl[2])
-        self.ura[idx] = self.quality_idx(ucls, uval)
+        dorb = np.zeros(3)
+        dorb[0] = self.sval(dx, 15, self.dorb_scl[0])
+        dorb[1] = self.sval(dy, 13, self.dorb_scl[1])
+        dorb[2] = self.sval(dz, 13, self.dorb_scl[2])
+
+        self.lc[inet].iode[sat] = iodn
+        self.lc[inet].iodc[sat] = iodc
+        self.lc[inet].dorb[sat] = dorb
+        self.ura[sat] = self.quality_idx(ucls, uval)
+
+        if sat not in self.lc[inet].t0:
+            self.lc[inet].t0[sat] = {}
+
+        self.lc[inet].t0[sat][sCType.ORBIT] = self.time
         return i
 
     def decode_cssr_orb(self, msg, i, inet=0):
@@ -212,7 +219,7 @@ def decode_cssr_orb(self, msg, i, inet=0):
 
         self.iodssr_c[sCType.ORBIT] = head['iodssr']
         self.lc[inet].cstat |= (1 << sCType.ORBIT)
-        self.lc[inet].t0[sCType.ORBIT] = self.time
+
         i += 19
         return i
 
@@ -233,30 +240,36 @@ def decode_cssr_cbias(self, msg, i, inet=0):
             i += 13
             sys, prn = self.slot2prn(slot)
             sat = prn2sat(sys, prn)
-            idx = np.where(sat == sat_n)[0]
-            if len(idx) == 0:
+            if sat not in sat_n:
                 continue
+            if sat not in self.lc[inet].t0:
+                self.lc[inet].t0[sat] = {}
+            self.lc[inet].t0[sat][sCType.CBIAS] = self.time
+
+            self.lc[inet].cbias[sat] = {}
             for j in range(0, nsig):
                 sig, cb = bs.unpack_from('u4s12', msg, i)
                 i += 16
-                self.sig_n[idx, j] = sig
-                self.lc[inet].cbias[idx, j] = self.sval(cb, 12, self.cb_scl)
+                sig_ = self.ssig2rsig(sys, uTYP.C, sig)
+                self.lc[inet].cbias[sat][sig_] = self.sval(cb, 12, self.cb_scl)
 
         self.iodssr_c[sCType.CBIAS] = head['iodssr']
         self.lc[inet].cstat |= (1 << sCType.CBIAS)
-        self.lc[inet].t0[sCType.CBIAS] = self.time
+
         return i
 
-    def decode_cssr_clk_sat(self, msg, i, inet, idx):
+    def decode_cssr_clk_sat(self, msg, i, inet, sat):
         """ decode clock correction for satellite """
         iodc, dclk = bs.unpack_from('u3s15', msg, i)
         i += 18
-        if iodc != self.lc[inet].iodc_c[idx]:
-            self.lc[inet].iodc_c_p[idx] = self.lc[inet].iodc_c[idx]
-            self.lc[inet].dclk_p[idx] = self.lc[inet].dclk[idx]
-        self.lc[inet].iodc_c[idx] = iodc
+
+        if sat in self.lc[inet].iodc_c.keys() and \
+                iodc != self.lc[inet].iodc_c[sat]:
+            self.lc[inet].iodc_c_p[sat] = self.lc[inet].iodc_c[sat]
+            self.lc[inet].dclk_p[sat] = self.lc[inet].dclk[sat]
+        self.lc[inet].iodc_c[sat] = iodc
         # note: the sign of the clock correction reversed
-        self.lc[inet].dclk[idx] = -self.sval(dclk, 15, self.dclk_scl)
+        self.lc[inet].dclk[sat] = -self.sval(dclk, 15, self.dclk_scl)
         return i
 
     def decode_cssr_clk(self, msg, i, inet=0):
@@ -273,11 +286,15 @@ def decode_cssr_clk(self, msg, i, inet=0):
         for k in range(23):
             idx = st1*23+k
             if idx < self.nsat_n:
-                i = self.decode_cssr_clk_sat(msg, i, inet, idx)
+                sat = self.sat_n[idx]
+                i = self.decode_cssr_clk_sat(msg, i, inet, sat)
+                if sat not in self.lc[inet].t0:
+                    self.lc[inet].t0[sat] = {}
+                self.lc[inet].t0[sat][sCType.CLOCK] = self.time
 
         self.iodssr_c[sCType.CLOCK] = head['iodssr']
         self.lc[inet].cstat |= (1 << sCType.CLOCK)
-        self.lc[inet].t0[sCType.CLOCK] = self.time
+
         i += 10
         return i
 
@@ -295,10 +312,11 @@ def decode_cssr_ura(self, msg, i):
             i += 6
             idx = st2*70+k
             if idx < self.nsat_n:
-                self.ura[idx] = self.quality_idx(v['class'], v['val'])
+                sat = self.sat_n[idx]
+                self.ura[sat] = self.quality_idx(v['class'], v['val'])
 
         self.lc[0].cstat |= (1 << sCType.URA)
-        self.lc[0].t0[sCType.URA] = self.time
+        # self.lc[0].t0[0][sCType.URA] = self.time
         return i
 
     def decode_cssr_comb1(self, msg, i, inet=0):
@@ -312,7 +330,8 @@ def decode_cssr_comb1(self, msg, i, inet=0):
             i += 36
             for k in range(numc):
                 idx = slot_s+k
-                i = self.decode_cssr_clk_sat(msg, i, inet, idx)
+                sat = self.sat_n[idx]
+                i = self.decode_cssr_clk_sat(msg, i, inet, sat)
 
         if numo > 0:
             tod, _, iodssr = bs.unpack_from('u17u4u2', msg, i)
@@ -335,8 +354,7 @@ def decode_cssr_comb2(self, msg, i, inet=0):
                 i += 9
                 sys, prn = self.slot2prn(slot)
                 sat = prn2sat(sys, prn)
-                idx = np.where(sat == sat_n)[0]
-                i = self.decode_cssr_clk_sat(msg, i, inet, idx)
+                i = self.decode_cssr_clk_sat(msg, i, inet, sat)
 
         if numo > 0:
             tod, _, iodssr = bs.unpack_from('u17u4u2', msg, i)
diff --git a/src/cssrlib/cssr_pvs.py b/src/cssrlib/cssr_pvs.py
new file mode 100644
index 0000000..3ecb043
--- /dev/null
+++ b/src/cssrlib/cssr_pvs.py
@@ -0,0 +1,176 @@
+"""
+PVS (PPP via SouthPAN) correction data decoder
+
+[1] Service Definition Document for Open Services, SBAS-STN-0001,
+    Revision 02, December 2022
+
+[2] Service Definition Document for Data Access Services, SBAS-STN-0002,
+    Revision 01, November 2023
+
+"""
+
+import numpy as np
+import bitstruct as bs
+from cssrlib.cssrlib import cssr, sCSSR, sCSSRTYPE, prn2sat, sCType
+from cssrlib.cssrlib import sat2id
+from cssrlib.gnss import uGNSS, rCST, gpst2time, timediff, timeadd, time2gpst
+from binascii import unhexlify
+
+
+def decode_sinca_line(line):
+    """ SINCA (SISNeT compression algorithm) decoder  """
+    if line[0:4] != '*MSG':
+        return None
+    v = line.split(',')
+    s = v[3]
+    week, tow = int(v[1]), int(v[2])
+    t = gpst2time(week, tow)
+    for key in '|/':
+        if key in s:
+            while True:
+                pos = s.find(key)
+                if pos < 0:
+                    break
+                c = s[pos-1]
+                k = 2 if key == '|' else 3
+                n = int(s[pos+1:pos+k], 16)
+                s = s.replace(s[pos-1:pos+k], c*n, 1)
+
+    sb = unhexlify(s.split('*')[0])
+    return t, sb
+
+
+class cssr_pvs(cssr):
+    def __init__(self, foutname=None):
+        super().__init__(foutname)
+        self.MAXNET = 1
+        self.cssrmode = sCSSRTYPE.PVS_PPP
+        self.nsig_max = 0
+        self.sat_n = []
+
+        self.lc[0].dclk = {}
+        self.lc[0].ddft = {}
+        self.lc[0].dorb = {}
+        self.lc[0].dvel = {}
+        self.lc[0].iode = {}
+        self.lc[0].t0 = {}
+
+        self.tmax = {sCType.CLOCK: 120.0, sCType.ORBIT: 120.0}
+
+        # GPS LNAV, Galileo F/NAV
+        self.nav_mode = {uGNSS.GPS: 0, uGNSS.GAL: 1}
+
+    def slot2sat(self, slot):
+        sat = 0
+        if slot >= 1 and slot <= 32:
+            prn = slot
+            sat = prn2sat(uGNSS.GPS, prn)
+        elif slot >= 38 and slot <= 69:
+            prn = slot-37
+            sat = prn2sat(uGNSS.GLO, prn)
+        elif slot >= 75 and slot <= 110:
+            prn = slot-74
+            sat = prn2sat(uGNSS.GAL, prn)
+        elif slot >= 120 and slot <= 158:
+            prn = slot
+            sat = prn2sat(uGNSS.SBS, prn)
+        elif slot >= 159 and slot <= 195:
+            prn = slot-158
+            sat = prn2sat(uGNSS.BDS, prn)
+        return sat
+
+    def adjust_time_week(self, time, time0):
+        dt = timediff(time, time0)
+        if dt > rCST.HALFWEEK_SEC:
+            time = timeadd(time, -rCST.WEEK_SEC)
+        elif dt < -rCST.HALFWEEK_SEC:
+            time = timeadd(time,  rCST.WEEK_SEC)
+        return time
+
+    def check_validity(self, time):
+        for sat in self.sat_n:
+            if timediff(time, self.lc[0].t0[sat][sCType.CLOCK]) > \
+                    self.tmax[sCType.CLOCK]:
+                self.lc[0].dclk[sat] = 0.0
+            if timediff(time, self.lc[0].t0[sat][sCType.ORBIT]) > \
+                    self.tmax[sCType.ORBIT]:
+                self.lc[0].iode[sat] = -1
+                self.lc[0].dorb[sat] = np.zeros(3)
+
+    def decode_cssr_orb(self, msg, i, inet=0):
+        """ Types 32 clock-ephemeris correction and covariance matrix """
+        slot, iodn = bs.unpack_from('u8u10', msg, i)
+        i += 18
+        dx, dy, dz, db = bs.unpack_from('s11s11s11s12', msg, i)
+        i += 45
+        dxd, dyd, dzd, dbd, t0 = bs.unpack_from('s8s8s8s9u13', msg, i)
+        i += 46
+
+        # i, C = self.decode_cov(msg, i)
+        i += 99
+        # dfrei, dRcorr = bs.unpack_from('u4u4', msg, i)
+        i += 8
+
+        sat = self.slot2sat(slot)
+
+        if sat == 0 or t0 % 2 == 0:  # for DFMC SBAS only
+            return
+
+        if sat not in self.sat_n:
+            self.sat_n.append(sat)
+
+        dorb = np.zeros(3)
+        dvel = np.zeros(3)
+
+        dorb[0] = dx*0.0625
+        dorb[1] = dy*0.0625
+        dorb[2] = dz*0.0625
+
+        dclk = db*0.03125
+
+        dvel[0] = dxd*0.00048828
+        dvel[1] = dyd*0.00048828
+        dvel[2] = dzd*0.00048828
+
+        ddft = dbd*0.00024414
+
+        tow = self.tow0 + t0*16.0
+        self.time = self.adjust_time_week(
+            gpst2time(self.week, tow), self.time0)
+
+        self.lc[0].iode[sat] = iodn
+        self.lc[0].dorb[sat] = -dorb
+        self.lc[0].dclk[sat] = dclk
+        self.lc[0].dvel[sat] = -dvel
+        self.lc[0].ddft[sat] = ddft
+
+        self.iodssr = 0
+        self.lc[0].cstat |= (1 << sCType.CLOCK) | (1 << sCType.ORBIT)
+        self.lc[0].t0[sat] = {
+            sCType.CLOCK: self.time0, sCType.ORBIT: self.time0}
+
+        self.iodssr_c[sCType.CLOCK] = self.iodssr
+        self.iodssr_c[sCType.ORBIT] = self.iodssr
+
+        if self.monlevel > 0:
+            self.fh.write("{:.1f} {:3s} {:3d}\n"
+                          .format(tow, sat2id(sat), iodn))
+
+        return i
+
+    def decode_cssr(self, msg, i=0):
+
+        _, mt = bs.unpack_from('u4u6', msg, i)
+        i += 10
+
+        if mt == 32:
+            self.subtype = sCSSR.ORBIT
+            i = self.decode_cssr_orb(msg, i)
+
+        if self.monlevel > 3:
+            if self.time != -1:
+                _, tow = time2gpst(self.time)
+            else:
+                tow = -1
+            self.fh.write("mt={:2d} tow={:6.1f}\n"
+                          .format(mt, tow))
diff --git a/src/cssrlib/cssrlib.py b/src/cssrlib/cssrlib.py
index ca1bb0f..8decf52 100644
--- a/src/cssrlib/cssrlib.py
+++ b/src/cssrlib/cssrlib.py
@@ -6,17 +6,18 @@
 import numpy as np
 from enum import IntEnum
 from cssrlib.gnss import gpst2time, rCST, prn2sat, uGNSS, gtime_t, rSigRnx
-from cssrlib.gnss import uSIG, uTYP, sat2prn, time2str, sat2id
+from cssrlib.gnss import uSIG, uTYP, sat2prn, time2str, sat2id, timediff
 
 
 class sCSSRTYPE(IntEnum):
-    QZS_CLAS = 0
-    QZS_MADOCA = 1
+    QZS_CLAS = 0     # QZS CLAS PPP-RTK
+    QZS_MADOCA = 1   # MADOCA-PPP
     GAL_HAS_SIS = 2  # Galileo HAS Signal-In-Space
     GAL_HAS_IDD = 3  # Galileo HAS Internet Data Distribution
-    BDS_PPP = 4
+    BDS_PPP = 4      # BDS PPP
     IGS_SSR = 5
     RTCM3_SSR = 6
+    PVS_PPP = 7      # PPP via SouthPAN
 
 
 class sGNSS(IntEnum):
@@ -285,6 +286,7 @@ def __init__(self, foutname=None):
             self.lc[inet].flg_trop = 0
             self.lc[inet].flg_stec = 0
             self.lc[inet].nsat_n = 0
+            self.lc[inet].t0 = {}
 
         self.dorb_scl = [0.0016, 0.0064, 0.0064]
         self.dclk_scl = 0.0016
@@ -298,9 +300,14 @@ def __init__(self, foutname=None):
         self.iodssr_c = np.ones(16, dtype=np.int32)*-1
         self.sig_n_p = []
 
+        # maximum validty time for correction
+        self.tmax = {sCType.CLOCK: 30.0, sCType.ORBIT: 120.0,
+                     sCType.CBIAS: 120.0, sCType.PBIAS: 120.0,
+                     sCType.TROP: 120.0, sCType.STEC: 120.0}
+
         # default navigation message mode: 0:LNAV/INAV, 1: CNAV/CNAV1
         self.nav_mode = {uGNSS.GPS: 0, uGNSS.QZS: 0,
-                         uGNSS.GAL: 0, uGNSS.BDS: 1}
+                         uGNSS.GAL: 0, uGNSS.BDS: 1, uGNSS.GLO: 0}
 
         self.fh = None
         if foutname is not None:
@@ -432,6 +439,12 @@ def isset(self, mask, nbit, k):
             return True
         return False
 
+    def set_t0(self, inet=0, sat=0, ctype=0, t=gtime_t()):
+        """ set reference time for correcion to check validity time """
+        if sat not in self.lc[inet].t0:
+            self.lc[inet].t0[sat] = {}
+        self.lc[inet].t0[sat][ctype] = t
+
     def quality_idx(self, cl, val):
         """ calculate quality index """
         if cl == 7 and val == 7:
@@ -509,18 +522,18 @@ def decode_cssr_mask(self, msg, i):
         self.gnss_n = []
         self.sat_n = []
         self.nsig_total = 0
-        self.sig_n = []
+        self.sig_n = {}
         self.nsig_max = 0
         self.nm_idx = np.zeros(self.SYSMAX, dtype=int)
 
         self.dcm = np.ones(self.SYSMAX)  # delta clock multipliers for HAS
-        self.gnss_idx = np.zeros(self.ngnss, dtype=int)
+        self.gnss_idx = np.zeros(self.ngnss, dtype=uGNSS)
         self.nsat_g = np.zeros(self.SYSMAX, dtype=int)
 
         for j in range(self.ngnss):
             v = bs.unpack_from_dict('u4u40u16u1', ['gnssid', 'svmask',
                                                    'sigmask', 'cma'], msg, i)
-            self.gnss_idx[j] = v['gnssid']
+            self.gnss_idx[j] = self.gnss2sys(v['gnssid'])
             sys = self.gnss2sys(v['gnssid'])
             i += 61
             prn, nsat = self.decode_mask(v['svmask'], 40)
@@ -533,6 +546,8 @@ def decode_cssr_mask(self, msg, i):
 
             self.nsig_max = max(self.nsig_max, nsig)
 
+            sig_r = [self.ssig2rsig(sys, uTYP.L, s) for s in sig]
+
             for k in range(0, nsat):
                 if sys == uGNSS.QZS:
                     prn[k] += 192
@@ -542,14 +557,14 @@ def decode_cssr_mask(self, msg, i):
                 self.gnss_n.append(v['gnssid'])
                 if v['cma'] == 1:
                     sig_s, nsig_s = self.decode_mask(vc[k], nsig, 0)
-                    sig_n = [sig[i] for i in sig_s]
+                    sig_n = [sig_r[i] for i in sig_s]
                     self.nsig_n.append(nsig_s)
                     self.nsig_total = self.nsig_total+nsig_s
-                    self.sig_n.append(sig_n)
+                    self.sig_n[sat] = sig_n
                 else:
                     self.nsig_n.append(nsig)
                     self.nsig_total = self.nsig_total+nsig
-                    self.sig_n.append(sig)
+                    self.sig_n[sat] = sig_r
 
             if self.cssrmode == sCSSRTYPE.GAL_HAS_SIS:  # HAS only
                 self.nm_idx[v['gnssid']] = bs.unpack_from('u3', msg, i)[0]
@@ -559,58 +574,58 @@ def decode_cssr_mask(self, msg, i):
             i += 6
 
         self.lc[0].cstat |= (1 << sCType.MASK)
-        self.lc[0].t0[sCType.MASK] = self.time
+        self.set_t0(ctype=sCType.MASK, t=self.time)
         return i
 
-    def decode_orb_sat(self, msg, i, k, sys=uGNSS.NONE, inet=0):
+    def decode_orb_sat(self, msg, i, sat, sys=uGNSS.NONE, inet=0):
         """ decoder orbit correction of cssr """
         n = 10 if sys == uGNSS.GAL else 8
         fmt = 'u{:d}s{:d}s{:d}s{:d}'\
             .format(n, self.dorb_blen[0], self.dorb_blen[1], self.dorb_blen[2])
-        v = bs.unpack_from_dict(fmt, ['iode', 'dx', 'dy', 'dz'], msg, i)
+        iode, dx, dy, dz = bs.unpack_from(fmt, msg, i)
 
-        self.lc[inet].iode[k] = v['iode']
-        self.lc[inet].dorb[k, 0] = \
-            self.sval(v['dx'], self.dorb_blen[0], self.dorb_scl[0])
-        self.lc[inet].dorb[k, 1] = \
-            self.sval(v['dy'], self.dorb_blen[1], self.dorb_scl[1])
-        self.lc[inet].dorb[k, 2] = \
-            self.sval(v['dz'], self.dorb_blen[2], self.dorb_scl[2])
+        dorb = np.zeros(3)
+        dorb[0] = self.sval(dx, self.dorb_blen[0], self.dorb_scl[0])
+        dorb[1] = self.sval(dy, self.dorb_blen[1], self.dorb_scl[1])
+        dorb[2] = self.sval(dz, self.dorb_blen[2], self.dorb_scl[2])
+
+        self.lc[inet].iode[sat] = iode
+        self.lc[inet].dorb[sat] = dorb
 
         if self.cssrmode == sCSSRTYPE.GAL_HAS_SIS:  # HAS SIS
-            self.lc[inet].dorb[k, :] *= -1.0
+            self.lc[inet].dorb[sat] *= -1.0
 
         i += n + self.dorb_blen[0]+self.dorb_blen[1]+self.dorb_blen[2]
         return i
 
-    def decode_clk_sat(self, msg, i, k, inet=0):
+    def decode_clk_sat(self, msg, i, sat, inet=0):
         """ decoder clock correction of cssr """
-        v = bs.unpack_from_dict('s'+str(self.dclk_blen), ['dclk'], msg, i)
+        dclk = bs.unpack_from('s'+str(self.dclk_blen), msg, i)[0]
 
-        self.lc[inet].dclk[k] = \
-            self.sval(v['dclk'], self.dclk_blen, self.dclk_scl)
+        self.lc[inet].dclk[sat] = self.sval(
+            dclk, self.dclk_blen, self.dclk_scl)
 
         if self.cssrmode == sCSSRTYPE.GAL_HAS_SIS:  # HAS SIS
-            self.lc[inet].dclk[k] *= self.dcm[self.gnss_n[k]]
+            sys, _ = sat2prn(sat)
+            self.lc[inet].dclk[sat] *= self.dcm[sys]
 
         i += self.dclk_blen
         return i
 
-    def decode_cbias_sat(self, msg, i, k, j, inet=0):
+    def decode_cbias_sat(self, msg, i, sat, rsig, inet=0):
         """ decoder code bias correction of cssr """
-        v = bs.unpack_from_dict('s'+str(self.cb_blen), ['cbias'], msg, i)
-        self.lc[inet].cbias[k, j] = \
-            self.sval(v['cbias'], self.cb_blen, self.cb_scl)
+        cb = bs.unpack_from('s'+str(self.cb_blen), msg, i)[0]
+        self.lc[inet].cbias[sat][rsig] = self.sval(
+            cb, self.cb_blen, self.cb_scl)
         i += self.cb_blen
         return i
 
-    def decode_pbias_sat(self, msg, i, k, j, inet=0):
+    def decode_pbias_sat(self, msg, i, sat, rsig, inet=0):
         """ decoder phase bias correction of cssr """
-        v = bs.unpack_from_dict('s'+str(self.pb_blen) +
-                                'u2', ['pbias', 'di'], msg, i)
-        self.lc[inet].pbias[k, j] = \
-            self.sval(v['pbias'], self.pb_blen, self.pb_scl)
-        self.lc[inet].di[k, j] = v['di']
+        pb, di = bs.unpack_from('s'+str(self.pb_blen) + 'u2', msg, i)
+        self.lc[inet].pbias[sat][rsig] = self.sval(
+            pb, self.pb_blen, self.pb_scl)
+        self.lc[inet].di[sat][rsig] = di
         i += self.pb_blen + 2
         return i
 
@@ -618,54 +633,45 @@ def decode_cssr_orb(self, msg, i, inet=0):
         """decode MT4073,2 Orbit Correction message """
         head, i = self.decode_head(msg, i)
         self.flg_net = False
-        #if self.iodssr != head['iodssr']:
-        #    return -1
-        dorb_p = self.lc[inet].dorb
-        self.lc[inet].dorb = np.zeros((self.nsat_n, 3))
-        self.lc[inet].iode = np.zeros(self.nsat_n, dtype=int)
-        self.lc[inet].dorb_d = np.ones((self.nsat_n, 3))*np.nan
-        for k in range(0, self.nsat_n):
-            i = self.decode_orb_sat(msg, i, k, self.sys_n[k], inet)
-            if self.sat_n[k] in self.sat_n_p:
-                j = self.sat_n_p.index(self.sat_n[k])
-                #self.lc[inet].dorb_d[k, :] = self.lc[inet].dorb[k, :] \
-                #    - dorb_p[j, :]
+        self.lc[inet].dorb = {}
+        self.lc[inet].iode = {}
+
+        for k in range(self.nsat_n):
+            i = self.decode_orb_sat(msg, i, self.sat_n[k], self.sys_n[k], inet)
+            self.set_t0(inet, self.sat_n[k], sCType.ORBIT, self.time)
 
         self.iodssr_c[sCType.ORBIT] = head['iodssr']
         self.lc[inet].cstat |= (1 << sCType.ORBIT)
-        self.lc[inet].t0[sCType.ORBIT] = self.time
+
         return i
 
     def decode_cssr_clk(self, msg, i, inet=0):
         """decode MT4073,3 Clock Correction message """
         head, i = self.decode_head(msg, i)
         self.flg_net = False
-        #if self.iodssr != head['iodssr']:
+        # if self.iodssr != head['iodssr']:
         #    return -1
 
         if (self.lc[0].cstat & (1 << sCType.MASK)) != (1 << sCType.MASK):
             return -1
 
         if self.cssrmode == sCSSRTYPE.GAL_HAS_SIS:  # HAS only
+            self.dcm = {}
             for k in range(self.ngnss):
-                self.dcm[self.gnss_idx[k]] = \
-                    bs.unpack_from('u2', msg, i)[0]+1.0
+                self.dcm[self.gnss_idx[k]] = bs.unpack_from('u2', msg, i)[
+                    0]+1.0
                 i += 2
 
-        dclk_p = self.lc[inet].dclk
-        self.lc[inet].dclk = np.zeros(self.nsat_n)
-        self.lc[inet].dclk_d = np.ones(self.nsat_n)*np.nan
-        for k in range(0, self.nsat_n):
-            i = self.decode_clk_sat(msg, i, k, inet)
-            # if self.sat_n[k] in self.sat_n_p:
-            #    j = self.sat_n_p.index(self.sat_n[k])
-            #    self.lc[inet].dclk_d[k] = self.lc[inet].dclk[k]-dclk_p[j]
+        self.lc[inet].dclk = {}
+        for k in range(self.nsat_n):
+            i = self.decode_clk_sat(msg, i, self.sat_n[k], inet)
+            self.set_t0(inet, self.sat_n[k], sCType.CLOCK, self.time)
 
         if self.cssrmode == sCSSRTYPE.GAL_HAS_SIS:  # HAS only
             self.sat_n_p = self.sat_n
         self.iodssr_c[sCType.CLOCK] = head['iodssr']
         self.lc[inet].cstat |= (1 << sCType.CLOCK)
-        self.lc[inet].t0[sCType.CLOCK] = self.time
+
         return i
 
     def decode_cssr_cbias(self, msg, i, inet=0):
@@ -673,16 +679,19 @@ def decode_cssr_cbias(self, msg, i, inet=0):
         head, i = self.decode_head(msg, i)
         nsat = self.nsat_n
         self.flg_net = False
-        #if self.iodssr != head['iodssr']:
-        #    return -1
-        self.lc[inet].cbias = np.zeros((nsat, self.nsig_max))
+        self.lc[inet].cbias = {}
         for k in range(nsat):
-            for j in range(0, self.nsig_n[k]):
-                i = self.decode_cbias_sat(msg, i, k, j, inet)
+            sat = self.sat_n[k]
+            sys, _ = sat2prn(sat)
+            self.lc[inet].cbias[sat] = {}
+            for j in range(len(self.sig_n[sat])):
+                rsig = self.sig_n[sat][j].toTyp(uTYP.C)
+                i = self.decode_cbias_sat(msg, i, sat, rsig, inet)
+                self.set_t0(inet, self.sat_n[k], sCType.CBIAS, self.time)
 
         self.iodssr_c[sCType.CBIAS] = head['iodssr']
         self.lc[inet].cstat |= (1 << sCType.CBIAS)
-        self.lc[inet].t0[sCType.CBIAS] = self.time
+
         return i
 
     def decode_cssr_pbias(self, msg, i, inet=0):
@@ -690,24 +699,28 @@ def decode_cssr_pbias(self, msg, i, inet=0):
         head, i = self.decode_head(msg, i)
         nsat = self.nsat_n
         self.flg_net = False
-        #if self.iodssr != head['iodssr']:
-        #    return -1
-        self.lc[inet].pbias = np.zeros((nsat, self.nsig_max))
-        self.lc[inet].di = np.zeros((nsat, self.nsig_max), dtype=int)
+        self.lc[inet].pbias = {}
+        self.lc[inet].di = {}
         for k in range(nsat):
-            for j in range(0, self.nsig_n[k]):
-                i = self.decode_pbias_sat(msg, i, k, j, inet)
+            sat = self.sat_n[k]
+            sys, _ = sat2prn(sat)
+            self.lc[inet].pbias[sat] = {}
+            self.lc[inet].di[sat] = {}
+            for j in range(len(self.sig_n[sat])):
+                rsig = self.sig_n[sat][j].toTyp(uTYP.L)
+                i = self.decode_pbias_sat(msg, i, sat, rsig, inet)
+                self.set_t0(inet, self.sat_n[k], sCType.PBIAS, self.time)
 
         self.iodssr_c[sCType.PBIAS] = head['iodssr']
         self.lc[inet].cstat |= (1 << sCType.PBIAS)
-        self.lc[inet].t0[sCType.PBIAS] = self.time
+
         return i
 
     def decode_cssr_bias(self, msg, i, inet=0):
         """decode MT4073,6 Bias Correction message """
         nsat = self.nsat_n
         head, i = self.decode_head(msg, i)
-        #if self.iodssr != head['iodssr']:
+        # if self.iodssr != head['iodssr']:
         #    return -1
         dfm = bs.unpack_from_dict('b1b1b1', ['cb', 'pb', 'net'], msg, i)
         self.flg_net = dfm['net']
@@ -721,29 +734,40 @@ def decode_cssr_bias(self, msg, i, inet=0):
             self.lc[inet].nsat_n = nsat = len(self.lc[inet].sat_n)
 
         if dfm['cb']:
-            self.lc[inet].cbias = np.zeros((nsat, self.nsig_max))
+            self.lc[inet].cbias = {}
         if dfm['pb']:
-            self.lc[inet].pbias = np.zeros((nsat, self.nsig_max))
-            self.lc[inet].di = np.zeros((nsat, self.nsig_max), dtype=int)
+            self.lc[inet].pbias = {}
+            self.lc[inet].di = {}
         ki = 0
         for k in range(self.nsat_n):
             if not self.isset(v['svmaskn'], self.nsat_n, k):
                 continue
+            sat = self.sat_n[k]
+            if dfm['cb']:
+                self.lc[inet].cbias[sat] = {}
+            if dfm['pb']:
+                self.lc[inet].pbias[sat] = {}
+                self.lc[inet].di[sat] = {}
             for j in range(self.nsig_n[k]):
                 if dfm['cb']:
-                    i = self.decode_cbias_sat(msg, i, ki, j, inet)
+                    rsig = self.sig_n[sat][j].toTyp(uTYP.C)
+                    i = self.decode_cbias_sat(msg, i, sat, rsig, inet)
+                    self.set_t0(inet, self.sat_n[k], sCType.CBIAS, self.time)
+
                 if dfm['pb']:
-                    i = self.decode_pbias_sat(msg, i, ki, j, inet)
+                    rsig = self.sig_n[sat][j].toTyp(uTYP.L)
+                    i = self.decode_pbias_sat(msg, i, sat, rsig, inet)
+                    self.set_t0(inet, self.sat_n[k], sCType.PBIAS, self.time)
+
             ki += 1
 
         if dfm['cb']:
             self.iodssr_c[sCType.CBIAS] = head['iodssr']
             self.lc[inet].cstat |= (1 << sCType.CBIAS)
-            self.lc[inet].t0[sCType.CBIAS] = self.time
         if dfm['pb']:
             self.iodssr_c[sCType.PBIAS] = head['iodssr']
             self.lc[inet].cstat |= (1 << sCType.PBIAS)
-            self.lc[inet].t0[sCType.PBIAS] = self.time
+
         return i
 
     def decode_cssr_ura(self, msg, i):
@@ -751,13 +775,15 @@ def decode_cssr_ura(self, msg, i):
         head, i = self.decode_head(msg, i)
         if self.iodssr != head['iodssr']:
             return -1
-        self.ura = np.zeros(self.nsat_n)
-        for k in range(0, self.nsat_n):
+        self.ura = {}
+        for k in range(self.nsat_n):
+            sat = self.sat_n[k]
             v = bs.unpack_from_dict('u3u3', ['class', 'val'], msg, i)
-            self.ura[k] = self.quality_idx(v['class'], v['val'])
+            self.ura[sat] = self.quality_idx(v['class'], v['val'])
             i += 6
         self.lc[0].cstat |= (1 << sCType.URA)
-        self.lc[0].t0[sCType.URA] = self.time
+        self.set_t0(ctype=sCType.URA, t=self.time)
+
         return i
 
     def decode_cssr_stec_coeff(self, msg, stype, i):
@@ -795,17 +821,20 @@ def decode_cssr_stec(self, msg, i):
         self.lc[inet].sat_n = self.decode_local_sat(netmask)
         self.lc[inet].nsat_n = nsat = len(self.lc[inet].sat_n)
         i += 7+self.nsat_n
-        self.lc[inet].stec_quality = np.zeros(nsat)
-        self.lc[inet].ci = np.zeros((nsat, 6))
+        self.lc[inet].stec_quality = {}
+        self.lc[inet].ci = {}
         for k in range(nsat):
+            sat = self.lc[inet].sat_n[k]
             v = bs.unpack_from_dict('u3u3', ['class', 'val'], msg, i)
-            self.lc[inet].stec_quality[k] = self.quality_idx(v['class'],
-                                                             v['val'])
+            self.lc[inet].stec_quality[sat] = self.quality_idx(v['class'],
+                                                               v['val'])
             i += 6
             ci, i = self.decode_cssr_stec_coeff(msg, dfm['stype'], i)
-            self.lc[inet].ci[k, :] = ci
+            self.lc[inet].ci[sat] = ci
+            self.set_t0(inet, sat, sCType.STEC, self.time)
+
         self.lc[inet].cstat |= (1 << sCType.STEC)
-        self.lc[inet].t0[sCType.STEC] = self.time
+
         return i
 
     def decode_cssr_grid(self, msg, i):
@@ -815,7 +844,7 @@ def decode_cssr_grid(self, msg, i):
             return -1
         dfm = bs.unpack_from_dict('u2u1u5u'+str(self.nsat_n)+'u3u3u6',
                                   ['ttype', 'range', 'inet', 'svmaskn',
-                                   'class', 'value', 'ng'], msg, i)
+                                  'class', 'value', 'ng'], msg, i)
         self.flg_net = True
         inet = dfm['inet']
         self.netmask[inet] = netmask = dfm['svmaskn']
@@ -828,11 +857,12 @@ def decode_cssr_grid(self, msg, i):
         i += 20+self.nsat_n
         sz = 7 if dfm['range'] == 0 else 16
         fmt = 's'+str(sz)
-        self.lc[inet].stec = np.zeros((ng, nsat))
+        self.lc[inet].stec = np.zeros(ng, dtype=list)
         self.lc[inet].dtd = np.zeros(ng)
         self.lc[inet].dtw = np.zeros(ng)
 
-        for j in range(0, ng):
+        for j in range(ng):
+            self.lc[inet].stec[j] = {}
             if dfm['ttype'] > 0:
                 vd = bs.unpack_from_dict('s9s8', ['dtd', 'dtw'], msg, i)
                 i += 17
@@ -840,11 +870,14 @@ def decode_cssr_grid(self, msg, i):
                 self.lc[inet].dtw[j] = self.sval(vd['dtw'], 8, 0.004)
 
             for k in range(nsat):
+                sat = self.lc[inet].sat_n[k]
                 dstec = bs.unpack_from(fmt, msg, i)[0]
                 i += sz
-                self.lc[inet].stec[j, k] = self.sval(dstec, sz, 0.04)
+                self.lc[inet].stec[j][sat] = self.sval(dstec, sz, 0.04)
         self.lc[inet].cstat |= (1 << sCType.TROP)
-        self.lc[inet].t0[sCType.TROP] = self.time
+        for k in range(nsat):
+            self.set_t0(inet, self.lc[inet].sat_n[k], sCType.TROP, self.time)
+
         return i
 
     def parse_sinfo(self):
@@ -871,7 +904,7 @@ def decode_cssr_sinfo(self, msg, i):
     def decode_cssr_comb(self, msg, i, inet=0):
         """decode MT4073,11 Orbit,Clock Combined Correction message """
         head, i = self.decode_head(msg, i)
-        #if self.iodssr != head['iodssr']:
+        # if self.iodssr != head['iodssr']:
         #    return -1
         dfm = bs.unpack_from_dict('b1b1b1', ['orb', 'clk', 'net'], msg, i)
         i += 3
@@ -884,26 +917,31 @@ def decode_cssr_comb(self, msg, i, inet=0):
             i += 5+self.nsat_n
 
         if dfm['orb']:
-            self.lc[inet].dorb = np.zeros((self.nsat_n, 3))
-            self.lc[inet].iode = np.zeros(self.nsat_n, dtype=int)
+            self.lc[inet].dorb = {}
+            self.lc[inet].iode = {}
         if dfm['clk']:
-            self.lc[inet].dclk = np.zeros(self.nsat_n)
+            self.lc[inet].dclk = {}
 
         for k in range(self.nsat_n):
             if self.flg_net and not self.isset(svmask, self.nsat_n, k):
                 continue
+            sat = self.sat_n[k]
             if dfm['orb']:
-                i = self.decode_orb_sat(msg, i, k, self.sys_n[k], inet)
+                i = self.decode_orb_sat(msg, i, sat, self.sys_n[k], inet)
+                self.set_t0(inet, sat, sCType.ORBIT, self.time)
+
             if dfm['clk']:
-                i = self.decode_clk_sat(msg, i, k, inet)
+                i = self.decode_clk_sat(msg, i, sat, inet)
+                self.set_t0(inet, sat, sCType.CLOCK, self.time)
+
         if dfm['clk']:
             # self.iodssr_c[sCType.CLOCK] = head['iodssr']
             self.lc[inet].cstat |= (1 << sCType.CLOCK)
-            self.lc[inet].t0[sCType.CLOCK] = self.time
+
         if dfm['orb']:
             # self.iodssr_c[sCType.ORBIT] = head['iodssr']
             self.lc[inet].cstat |= (1 << sCType.ORBIT)
-            self.lc[inet].t0[sCType.ORBIT] = self.time
+
         return i
 
     def decode_cssr_atmos(self, msg, i):
@@ -959,25 +997,29 @@ def decode_cssr_atmos(self, msg, i):
         self.lc[inet].netmask = netmask
         self.lc[inet].sat_n = self.decode_local_sat(netmask)
         self.lc[inet].nsat_n = nsat = len(self.lc[inet].sat_n)
-        self.lc[inet].stec_quality = np.zeros(nsat)
+        self.lc[inet].stec_quality = {}
         if dfm['stec'] & 2 > 0:
-            self.lc[inet].ci = np.zeros((nsat, 6))
-            self.lc[inet].stype = np.zeros(nsat, dtype=int)
+            self.lc[inet].ci = {}
+            self.lc[inet].stype = {}
         if dfm['stec'] & 1 > 0:
-            self.lc[inet].dstec = np.zeros((nsat, ng))
+            self.lc[inet].dstec = {}
 
         for k in range(nsat):
+            sat = self.lc[inet].sat_n[k]
+            self.set_t0(inet, sat, sCType.STEC, self.time)
+
             if dfm['stec'] > 0:
                 v = bs.unpack_from_dict('u3u3', ['class', 'value'], msg, i)
                 i += 6
-                self.lc[inet].stec_quality[k] = self.quality_idx(v['class'],
-                                                                 v['value'])
+                self.lc[inet].stec_quality[sat] = self.quality_idx(v['class'],
+                                                                   v['value'])
             if dfm['stec'] & 2 > 0:  # functional term
-                self.lc[inet].stype[k] = bs.unpack_from('u2', msg, i)[0]
+                self.lc[inet].stype[sat] = bs.unpack_from('u2', msg, i)[0]
                 i += 2
                 ci, i = self.decode_cssr_stec_coeff(msg,
-                                                    self.lc[inet].stype[k], i)
-                self.lc[inet].ci[k, :] = ci
+                                                    self.lc[inet].stype[sat],
+                                                    i)
+                self.lc[inet].ci[sat] = ci
 
             if dfm['stec'] & 1 > 0:  # residual term
                 sz_idx = bs.unpack_from('u2', msg, i)[0]
@@ -986,15 +1028,17 @@ def decode_cssr_atmos(self, msg, i):
                 scl = self.stec_scl_t[sz_idx]
                 v = bs.unpack_from(('s'+str(sz))*ng, msg, i)
                 i += sz*ng
+                self.lc[inet].dstec[sat] = np.zeros(ng)
                 for j in range(ng):
-                    self.lc[inet].dstec[k, j] = self.sval(v[j], sz, scl)
+                    self.lc[inet].dstec[sat][j] = self.sval(v[j], sz, scl)
 
         if dfm['trop'] > 0:
             self.lc[inet].cstat |= (1 << sCType.TROP)
-            self.lc[inet].t0[sCType.TROP] = self.time
+            self.set_t0(inet, 0, sCType.TROP, self.time)
+
         if dfm['stec'] > 0:
             self.lc[inet].cstat |= (1 << sCType.STEC)
-            self.lc[inet].t0[sCType.STEC] = self.time
+
         return i
 
     def decode_cssr(self, msg, i=0):
@@ -1055,6 +1099,17 @@ def chk_stat(self):
             return False
         return True
 
+    def check_validity(self, time):
+        """ check validity time of correction """
+        for sat in self.sat_n:
+            if timediff(time, self.lc[0].t0[sat][sCType.CLOCK]) > \
+                    self.tmax[sCType.CLOCK]:
+                self.lc[0].dclk[sat] = 0.0
+            if timediff(time, self.lc[0].t0[sat][sCType.ORBIT]) > \
+                    self.tmax[sCType.ORBIT]:
+                self.lc[0].iode[sat] = -1
+                self.lc[0].dorb[sat] = np.zeros(3)
+
     def out_log(self):
 
         if self.time == -1:
@@ -1067,70 +1122,64 @@ def out_log(self):
 
         if self.subtype == sCSSR.CLOCK:
             self.fh.write(" {:s}\t{:s}\n".format("SatID", "dclk [m]"))
-            for k, sat_ in enumerate(self.sat_n):
-                if np.isnan(self.lc[0].dclk[k]):
+            for k, sat_ in enumerate(self.lc[0].dclk.keys()):
+                if np.isnan(self.lc[0].dclk[sat_]):
                     continue
                 self.fh.write(" {:s}\t{:8.4f}\n"
                               .format(sat2id(sat_),
-                                      self.lc[0].dclk[k]))
+                                      self.lc[0].dclk[sat_]))
 
         elif self.subtype == sCSSR.ORBIT:
             self.fh.write(" {:s}\t{:s}\t{:s}\t{:s}\t{:s}\n"
                           .format("SatID", "IODE", "Radial[m]",
                                   "Along[m]", "Cross[m]"))
-            for k, sat_ in enumerate(self.sat_n):
-                if np.isnan(self.lc[0].dorb[k][0]):
+            for k, sat_ in enumerate(self.lc[0].dorb.keys()):
+                if np.isnan(self.lc[0].dorb[sat_][0]):
                     continue
                 self.fh.write(" {:s}\t{:3d}\t{:6.3f}\t{:6.3f}\t{:6.3f}\n"
                               .format(sat2id(sat_),
-                                      self.lc[0].iode[k],
-                                      self.lc[0].dorb[k][0],
-                                      self.lc[0].dorb[k][1],
-                                      self.lc[0].dorb[k][2]))
+                                      self.lc[0].iode[sat_],
+                                      self.lc[0].dorb[sat_][0],
+                                      self.lc[0].dorb[sat_][1],
+                                      self.lc[0].dorb[sat_][2]))
 
         elif self.subtype == sCSSR.COMBINED:
             self.fh.write(" {:s}\t{:s}\t{:s}\t{:s}\t{:s}\t{:s}\n"
                           .format("SatID", "IODE", "Radial[m]",
                                   "Along[m]", "Cross[m]", "dclk[m]"))
-            for k, sat_ in enumerate(self.sat_n):
-                if np.isnan(self.lc[0].dorb[k][0]) or \
-                   np.isnan(self.lc[0].dclk[k]):
+            for k, sat_ in enumerate(self.lc[0].dorb.keys()):
+                if np.isnan(self.lc[0].dorb[sat_][0]) or \
+                   np.isnan(self.lc[0].dclk[sat_]):
                     continue
                 self.fh.write(
                     " {:s}\t{:3d}\t{:6.3f}\t{:6.3f}\t{:6.3f}\t{:6.3f}\n"
                     .format(sat2id(sat_),
-                            self.lc[0].iode[k],
-                            self.lc[0].dorb[k][0],
-                            self.lc[0].dorb[k][1],
-                            self.lc[0].dorb[k][2],
-                            self.lc[0].dclk[k]))
+                            self.lc[0].iode[sat_],
+                            self.lc[0].dorb[sat_][0],
+                            self.lc[0].dorb[sat_][1],
+                            self.lc[0].dorb[sat_][2],
+                            self.lc[0].dclk[sat_]))
 
         elif self.subtype == sCSSR.CBIAS:
             self.fh.write(" {:s}\t{:s}\t{:s}\t{:s}\n"
                           .format("SatID", "SigID", "Code Bias[m]", "..."))
-            for k, sat_ in enumerate(self.sat_n):
-                sys_, _ = sat2prn(sat_)
+            for k, sat_ in enumerate(self.lc[0].cbias.keys()):
                 self.fh.write(" {:s}\t".format(sat2id(sat_)))
-                for j in range(self.nsig_n[k]):
-                    if self.sig_n[k][j] == -1:
-                        continue
-                    sig_ = self.ssig2rsig(sys_, uTYP.C, self.sig_n[k][j])
+                for sig_ in range(self.lc[0].cbias[sat_].keys()):
                     self.fh.write("{:s}\t{:5.2f}\t"
-                                  .format(sig_.str(), self.lc[0].cbias[k, j]))
+                                  .format(sig_.str(),
+                                          self.lc[0].cbias[sat_][sig_]))
                 self.fh.write("\n")
 
         elif self.subtype == sCSSR.PBIAS:
             self.fh.write(" {:s}\t{:s}\t{:s}\t{:s}\n"
                           .format("SatID", "SigID", "Phase Bias[m]", "..."))
-            for k, sat_ in enumerate(self.sat_n):
-                sys_, _ = sat2prn(sat_)
+            for k, sat_ in enumerate(self.lc[0].pbias.keys()):
                 self.fh.write(" {:s}\t".format(sat2id(sat_)))
-                for j in range(self.nsig_n[k]):
-                    if self.sig_n[k][j] == -1:
-                        continue
-                    sig_ = self.ssig2rsig(sys_, uTYP.L, self.sig_n[k][j])
+                for sig_ in range(self.lc[0].pbias[sat_].keys()):
                     self.fh.write("{:s}\t{:5.2f}\t"
-                                  .format(sig_.str(), self.lc[0].pbias[k, j]))
+                                  .format(sig_.str(),
+                                          self.lc[0].pbias[sat_][sig_]))
                 self.fh.write("\n")
 
         self.fh.flush()
@@ -1213,12 +1262,13 @@ def get_stec(self, dlat=0.0, dlon=0.0):
         nsat = self.lc[inet].nsat_n
         stec = np.zeros(nsat)
         for i in range(nsat):
+            sat = self.lc[inet].sat_n[i]
             if self.lc[inet].flg_stec & 2:
-                ci = self.lc[inet].ci[i, :]
+                ci = self.lc[inet].ci[sat]
                 stec[i] = [1, dlat, dlon, dlat*dlon, dlat**2, dlon**2]@ci
             if self.lc[inet].flg_stec & 1:
-                dstec = self.lc[inet].dstec[i,
-                                            self.grid_index-1]@self.grid_weight
+                dstec = self.lc[inet].dstec[sat][self.grid_index-1] \
+                    @ self.grid_weight
                 stec[i] += dstec
         return stec
 
diff --git a/src/cssrlib/ephemeris.py b/src/cssrlib/ephemeris.py
index 805e224..4294c89 100644
--- a/src/cssrlib/ephemeris.py
+++ b/src/cssrlib/ephemeris.py
@@ -4,6 +4,7 @@
 
 import numpy as np
 from cssrlib.gnss import uGNSS, rCST, sat2prn, timediff, timeadd, vnorm
+from cssrlib.gnss import gtime_t, Geph
 from cssrlib.cssrlib import sCSSRTYPE as sc
 from cssrlib.cssrlib import sCType
 
@@ -19,6 +20,8 @@ def findeph(nav, t, sat, iode=-1, mode=0):
         if eph_.sat != sat:
             continue
         dt = timediff(t, eph_.toe)
+        if sat == 43 and eph_.mode == 1:
+            None
         if (iode < 0 or eph_.iode == iode) and eph_.mode == mode and \
                 abs(dt) < dt_p:
             eph = eph_
@@ -36,6 +39,76 @@ def dtadjust(t1, t2, tw=604800):
     return dt
 
 
+def deq(x, acc):
+    xdot = np.zeros(6)
+
+    r2 = x[0:3]@x[0:3]
+    r3 = r2*np.sqrt(r2)
+    omg2 = rCST.OMGE_GLO**2
+
+    if r2 <= 0.0:
+        return xdot
+
+    a = 1.5*rCST.J2_GLO*rCST.MU_GLO*rCST.RE_GLO**2/r2/r3
+    b = 5.0*x[2]**2/r2
+    c = -rCST.MU_GLO/r3-a*(1.0-b)
+
+    xdot[0:3] = x[3:6]
+    xdot[3] = (c+omg2)*x[0]+2.0*rCST.OMGE_GLO*x[4]
+    xdot[4] = (c+omg2)*x[1]-2.0*rCST.OMGE_GLO*x[3]
+    xdot[5] = (c-2.0*a)*x[2]
+    xdot[3:6] += acc
+    return xdot
+
+
+def glorbit(t, x, acc):
+    k1 = deq(x, acc)
+    w = x + k1*t/2.0
+    k2 = deq(w, acc)
+    w = x + k2*t/2.0
+    k3 = deq(w, acc)
+    w = x + k3*t
+    k4 = deq(w, acc)
+    x += (k1+2.0*k2+2.0*k3+k4)*t/6.0
+    return x
+
+
+def geph2pos(time: gtime_t, geph: Geph, flg_v=False, TSTEP=1.0):
+    """ calculate GLONASS satellite position based on ephemeris """
+    t = timediff(time, geph.toe)
+    dts = -geph.taun+geph.gamn*t
+    x = np.zeros(6)
+    x[0:3] = geph.pos
+    x[3:6] = geph.vel
+
+    tt = -TSTEP if t < 0.0 else TSTEP
+
+    while True:
+        if np.fabs(t) <= 1e-9:
+            break
+        if np.fabs(t) < TSTEP:
+            tt = t
+        x = glorbit(tt, x, geph.acc)
+        t -= tt
+
+        rs = x[0:3]
+        vs = x[3:6]
+
+    if flg_v:
+        return rs, vs, dts
+
+    return rs, dts
+
+
+def geph2clk(time: gtime_t, geph: Geph):
+    """ calculate GLONASS satellite clock offset based on ephemeris """
+    ts = timediff(time, geph.toe)
+    t = ts
+    for i in range(2):
+        t = ts - (-geph.taun+geph.gamn*t)
+    return -geph.taun + geph.gamn*t
+
+
 def eph2pos(t, eph, flg_v=False):
     """ calculate satellite position based on ephemeris """
     sys, prn = sat2prn(eph.sat)
@@ -197,15 +270,18 @@ def satposs(obs, nav, cs=None, orb=None):
         if nav.ephopt == 4:
 
             rs_, dts_, _ = orb.peph2pos(t, sat, nav)
-            if rs_ is None or dts_ is None:
+            if rs_ is None or dts_ is None or np.isnan(dts_[0]):
                 continue
             dt = dts_[0]
 
-            eph = findeph(nav.eph, t, sat)
-            if eph is None:
-                svh[i] = 1
-                continue
-            svh[i] = eph.svh
+            if len(nav.eph) > 0:
+                eph = findeph(nav.eph, t, sat)
+                if eph is None:
+                    svh[i] = 1
+                    continue
+                svh[i] = eph.svh
+            else:
+                svh[i] = 0
 
         else:
 
@@ -214,51 +290,54 @@ def satposs(obs, nav, cs=None, orb=None):
                 if cs.iodssr_c[sCType.ORBIT] == cs.iodssr:
                     if sat not in cs.sat_n:
                         continue
-                    idx = cs.sat_n.index(sat)
-                else:
-                    if cs.iodssr_c[sCType.ORBIT] == cs.iodssr_p:
-                        if sat not in cs.sat_n_p:
-                            continue
-                        idx = cs.sat_n_p.index(sat)
-                    else:
+                elif cs.iodssr_c[sCType.ORBIT] == cs.iodssr_p:
+                    if sat not in cs.sat_n_p:
                         continue
+                else:
+                    continue
+
+                if sat not in cs.lc[0].iode.keys():
+                    continue
 
-                iode = cs.lc[0].iode[idx]
-                dorb = cs.lc[0].dorb[idx, :]  # radial,along-track,cross-track
+                iode = cs.lc[0].iode[sat]
+                dorb = cs.lc[0].dorb[sat]  # radial,along-track,cross-track
+
+                if cs.cssrmode == sc.PVS_PPP:
+                    dorb += cs.lc[0].dvel[sat] * \
+                        (timediff(obs.t, cs.lc[0].t0[sat][sCType.ORBIT]))
 
                 if cs.cssrmode == sc.BDS_PPP:  # consistency check for IOD corr
-                    if cs.lc[0].iodc[idx] == cs.lc[0].iodc_c[idx]:
-                        dclk = cs.lc[0].dclk[idx]
+
+                    if cs.lc[0].iodc[sat] == cs.lc[0].iodc_c[sat]:
+                        dclk = cs.lc[0].dclk[sat]
                     else:
-                        if cs.lc[0].iodc[idx] == cs.lc[0].iodc_c_p[idx]:
-                            dclk = cs.lc[0].dclk_p[idx]
+                        if cs.lc[0].iodc[sat] == cs.lc[0].iodc_c_p[sat]:
+                            dclk = cs.lc[0].dclk_p[sat]
                         else:
                             continue
 
                 else:
 
                     if cs.cssrmode == sc.GAL_HAS_SIS:  # HAS only
-
                         if cs.mask_id != cs.mask_id_clk:  # mask has changed
                             if sat not in cs.sat_n_p:
                                 continue
-                            idx = cs.sat_n_p.index(sat)
-
                     else:
-
                         if cs.iodssr_c[sCType.CLOCK] == cs.iodssr:
                             if sat not in cs.sat_n:
                                 continue
-                            idx = cs.sat_n.index(sat)
                         else:
                             if cs.iodssr_c[sCType.CLOCK] == cs.iodssr_p:
                                 if sat not in cs.sat_n_p:
                                     continue
-                                idx = cs.sat_n_p.index(sat)
                             else:
                                 continue
 
-                    dclk = cs.lc[0].dclk[idx]
+                    dclk = cs.lc[0].dclk[sat]
+
+                    if cs.cssrmode == sc.PVS_PPP:
+                        dclk += cs.lc[0].ddft[sat] * \
+                            (timediff(obs.t, cs.lc[0].t0[sat][sCType.CLOCK]))
 
                 if np.isnan(dclk) or np.isnan(dorb@dorb):
                     continue
@@ -269,27 +348,43 @@ def satposs(obs, nav, cs=None, orb=None):
 
                 mode = 0
 
-            eph = findeph(nav.eph, t, sat, iode, mode=mode)
-            if eph is None:
-                svh[i] = 1
-                continue
+            if sys == uGNSS.GLO:
+                geph = findeph(nav.geph, t, sat, iode, mode=mode)
+                if geph is None:
+                    svh[i] = 1
+                    continue
+
+                svh[i] = geph.svh
+                dt = geph2clk(t, geph)
+
+                if sat not in nav.glo_ch:
+                    nav.glo_ch[sat] = geph.frq
+
+            else:
+                eph = findeph(nav.eph, t, sat, iode, mode=mode)
+                if eph is None:
+                    svh[i] = 1
+                    continue
 
-            svh[i] = eph.svh
-            dt = eph2clk(t, eph)
+                svh[i] = eph.svh
+                dt = eph2clk(t, eph)
 
         t = timeadd(t, -dt)
 
         if nav.ephopt == 4:  # precise ephemeris
 
             rs_, dts_, _ = orb.peph2pos(t, sat, nav)
-            rs[i, :] = rs_[0:3]
-            vs[i, :] = rs_[3:6]
+            rs[i, :] = rs_[0: 3]
+            vs[i, :] = rs_[3: 6]
             dts[i] = dts_[0]
             nsat += 1
 
         else:
 
-            rs[i, :], vs[i, :], dts[i] = eph2pos(t, eph, True)
+            if sys == uGNSS.GLO:
+                rs[i, :], vs[i, :], dts[i] = geph2pos(t, geph, True)
+            else:
+                rs[i, :], vs[i, :], dts[i] = eph2pos(t, eph, True)
 
             # Apply SSR correction
             #
@@ -308,25 +403,33 @@ def satposs(obs, nav, cs=None, orb=None):
                     er = np.cross(ea, ec)
                     A = np.array([er, ea, ec])
 
-                dorb_e = dorb@A
+                if cs.cssrmode == sc.PVS_PPP:
+                    dorb_e = dorb
+                else:
+                    dorb_e = dorb@A
 
                 rs[i, :] -= dorb_e
                 dts[i] += dclk/rCST.CLIGHT
 
-                ers = vnorm(rs[i, :]-nav.x[0:3])
+                if cs.cssrmode == sc.PVS_PPP and sys == uGNSS.GPS:
+                    dts[i] -= eph.tgd
+
+                ers = vnorm(rs[i, :]-nav.x[0: 3])
                 dorb_ = -ers@dorb_e
                 sis = dclk-dorb_
-                if cs.lc[0].t0[1].time % 30 == 0 and \
-                        timediff(cs.lc[0].t0[1], nav.time_p) > 0:
+                if cs.lc[0].t0[sat][sCType.ORBIT].time % 30 == 0 and \
+                        timediff(cs.lc[0].t0[sat][sCType.ORBIT], nav.time_p) > 0:
                     if abs(nav.sis[sat]) > 0:
                         nav.dsis[sat] = sis - nav.sis[sat]
                     nav.sis[sat] = sis
 
                 nav.dorb[sat] = dorb_
                 nav.dclk[sat] = dclk
-                nsat += 1
+
+            nsat += 1
 
     if cs is not None:
-        nav.time_p = cs.lc[0].t0[1]
+        if sat in cs.lc[0].t0 and sCType.ORBIT in cs.lc[0].t0[sat]:
+            nav.time_p = cs.lc[0].t0[sat][sCType.ORBIT]
 
     return rs, vs, dts, svh, nsat
diff --git a/src/cssrlib/gnss.py b/src/cssrlib/gnss.py
index a1ca746..ba4b3e7 100644
--- a/src/cssrlib/gnss.py
+++ b/src/cssrlib/gnss.py
@@ -30,11 +30,14 @@ class rCST():
     OMGE_BDS = 7.2921150E-5
     RE_WGS84 = 6378137.0
     FE_WGS84 = (1.0/298.257223563)
+    RE_GLO = 6378136.0
     AU = 149597870691.0
     D2R = 0.017453292519943295
     R2D = 57.29577951308232
     AS2R = D2R/3600.0
     DAY_SEC = 86400.0
+    WEEK_SEC = 604800.0
+    HALFWEEK_SEC = 302400.0
     CENTURY_SEC = DAY_SEC*36525.0
 
     FREQ_G1 = 1575.42e6      # [Hz] GPS L1
@@ -416,8 +419,8 @@ def str2sig(self, sys, s):
                    (s[2] == '5' and s[2] not in 'IQX'):
                     raise ValueError
             elif sys == uGNSS.GLO:
-                if (s[1] == '1' and s[2] not in 'CP') or \
-                   (s[1] == '2' and s[2] not in 'CP') or \
+                if (s[1] == '1' and s[2] not in 'CPX') or \
+                   (s[1] == '2' and s[2] not in 'CPX') or \
                    (s[1] == '3' and s[2] not in 'IQX') or \
                    (s[1] == '4' and s[2] not in 'ABX') or \
                    (s[1] == '6' and s[2] not in 'ABX'):
@@ -638,9 +641,9 @@ class Eph():
     sattype = 0
     sismai = 0
     code = 0
-    urai = np.zeros(3)
-    sisai = np.zeros(4)
-    isc = np.zeros(6)
+    urai = None
+    sisai = None
+    isc = None
     integ = 0
     # 0:LNAV,INAV,D1/D2, 1:CNAV/CNAV1/FNAV, 2: CNAV2, 3: CNAV3, 4:FDMA, 5:SBAS
     mode = 0
@@ -659,13 +662,14 @@ class Geph():
     age = 0.0
     toe = gtime_t()
     tof = gtime_t()
-    pos = np.zeros(3)
-    vel = np.zeros(3)
-    acc = np.zeros(3)
+    pos = None
+    vel = None
+    acc = None
     taun = 0.0         # SV clock bias [s]
     gamn = 0.0         # relative frq bias
     dtaun = 0.0        # delta between L1 and L2 [s]
     mode = 0
+    status = 0
 
     def __init__(self, sat=0):
         self.sat = sat
@@ -695,6 +699,8 @@ class Nav():
 
     def __init__(self, nf=2):
         self.eph = []
+        self.geph = []
+        self.seph = []
         self.peph = []
         self.ion = np.array([
             [0.1118E-07, -0.7451E-08, -0.5961E-07, 0.1192E-06],
@@ -716,6 +722,24 @@ def __init__(self, nf=2):
         self.ephopt = 2  # ephemeris option 0: BRDC, 1: SBAS, 2: SSR-APC,
         #                  3: SSR-CG, 4: PREC
 
+        # 0:float-ppp,1:continuous,2:instantaneous,3:fix-and-hold
+        self.armode = 0
+        self.thresar = 3.0  # AR acceptance threshold
+        self.elmaskar = np.deg2rad(20.0)  # elevation mask for AR
+
+        # Select tropospheric model
+        #
+        self.trpModel = uTropoModel.SAAST
+
+        # 0: use trop-model, 1: estimate, 2: use cssr correction
+        self.trop_opt = 0
+
+        # 0: use iono-model, 1: estimate, 2: use cssr correction
+        self.iono_opt = 0
+
+        # 0: none, 1: full model, 2: local/regional model
+        self.phw_opt = 1
+
         self.monlevel = 1
         self.cnr_min = 25
         self.cnr_min_gpy = 15
@@ -742,6 +766,9 @@ def __init__(self, nf=2):
         self.tt = 0
         self.t = gtime_t()
 
+        # GLONASS frequency channel table
+        self.glo_ch = {}
+
 
 def epoch2time(ep):
     """ calculate time from epoch """
@@ -927,7 +954,7 @@ def time2epoch(t):
 
 
 def time2doy(t):
-    """ convert time to epoch """
+    """ convert time to day of year (as float value!) """
     ep = time2epoch(t)
     ep[1] = ep[2] = 1.0
     ep[3] = ep[4] = ep[5] = 0.0
@@ -951,6 +978,15 @@ def time2str(t):
         .format(e[0], e[1], e[2], e[3], e[4], int(e[5]))
 
 
+def adjtime(t: gtime_t, tref: gtime_t, dt=rCST.WEEK_SEC):
+    tt = timediff(t, tref)
+    if tt < -dt/2.0:
+        return timeadd(t, dt)
+    if tt > dt/2.0:
+        return timeadd(t, -dt)
+    return t
+
+
 def prn2sat(sys, prn):
     """ convert sys+prn to sat """
     if sys == uGNSS.GPS:
diff --git a/src/cssrlib/peph.py b/src/cssrlib/peph.py
index f041824..8dc6e66 100644
--- a/src/cssrlib/peph.py
+++ b/src/cssrlib/peph.py
@@ -5,9 +5,10 @@
 @author: ruihi
 """
 
-from cssrlib.gnss import Nav, id2sat, sat2id, char2sys
+from cssrlib.gnss import Nav, id2sat, sat2id, char2sys, sat2prn
 from cssrlib.gnss import epoch2time, time2epoch, timeadd, timediff, gtime_t
-from cssrlib.gnss import str2time, gpst2utc, utc2gpst
+from cssrlib.gnss import str2time, time2str, time2doy
+from cssrlib.gnss import gpst2utc, utc2gpst, time2gpst
 from cssrlib.gnss import ecef2enu
 from cssrlib.gnss import rCST, rSigRnx, uGNSS, uTYP, uSIG
 
@@ -28,8 +29,8 @@ def __init__(self, time=None):
             self.time = time
         else:
             self.time = gtime_t()
-        self.pos = np.zeros((uGNSS.MAXSAT, 4))
-        self.vel = np.zeros((uGNSS.MAXSAT, 4))
+        self.pos = np.ones((uGNSS.MAXSAT, 4))*np.nan
+        self.vel = np.ones((uGNSS.MAXSAT, 4))*np.nan
         self.std = np.zeros((uGNSS.MAXSAT, 4))
         self.vst = np.zeros((uGNSS.MAXSAT, 4))
 
@@ -175,7 +176,7 @@ def parse_sp3(self, fname, nav, opt=0):
                                     continue
 
                                 val = float(line[4+j*14:18+j*14])
-                                if val != 0.0 and abs(val-999999.999999) >= 1e-6:
+                                if abs(val-999999.999999) >= 1e-6:
                                     scl = 1e3 if j < 3 else 1e-6
                                     if line[0] == 'P':
                                         v = True
@@ -202,6 +203,92 @@ def parse_sp3(self, fname, nav, opt=0):
 
         return nav
 
+    def write_sp3(self, fname, nav):
+        """
+        Write data to SP3 file
+        """
+
+        with open(fname, "w") as fh:
+
+            # Write header section
+            #
+
+            # Epoch lines
+            #
+            t0 = nav.peph[0].time
+            e = time2epoch(t0)
+            ne = len(nav.peph)
+
+            fh.write("#dP{:04d} {:02d} {:02d} {:02d} {:02d} {:011.8f} {:7d} d+D {:16s}\n"
+                     .format(e[0], e[1], e[2], e[3], e[4], e[5], ne, ' '))
+
+            week, secs = time2gpst(t0)
+            tstep = timediff(nav.peph[1].time, t0)
+            mjd = 44244 + 7*week + int(secs/86400.0)
+            fod = time2doy(t0) % 1
+
+            fh.write("## {:04d} {:15.8f} {:14.8f} {:5n} {:15.13f}\n"
+                     .format(week, secs, tstep, mjd, fod))
+
+            # Satellite list and accuracy indicators
+            #
+            for i in range(int(np.ceil(self.nsat / 17))):
+
+                nsat = "{:4n}".format(self.nsat) if i == 0 else "    "
+                prns = [sat2id(s) for s in self.sat[i*17:i*17+17]]
+                fh.write('+ {}   {:51s}\n'.format(nsat, ''.join(prns)))
+
+            for i in range(int(self.nsat / 17+1)):
+
+                accs = ['  0' for s in self.sat[i*17:i*17+17]]
+                fh.write('++{}   {:51s}\n'.format('    ', ''.join(accs)))
+
+            fh.write(
+                '%c M  cc GPS ccc cccc cccc cccc cccc ccccc ccccc ccccc ccccc\n')
+            fh.write(
+                '%c cc cc ccc ccc cccc cccc cccc cccc ccccc ccccc ccccc ccccc\n')
+
+            fh.write(
+                '%f  1.2500000  1.025000000  0.00000000000  0.000000000000000\n')
+            fh.write(
+                '%f  0.0000000  0.000000000  0.00000000000  0.000000000000000\n')
+            fh.write(
+                '%i    0    0    0    0      0      0      0      0         0\n')
+            fh.write(
+                '%i    0    0    0    0      0      0      0      0         0\n')
+
+            # Comment section
+            #
+            fh.write('/* \n')
+
+            # Write data section
+            #
+            for peph in nav.peph:
+
+                e = time2epoch(peph.time)
+
+                fh.write("*  {:04d} {:02d} {:02d} {:02d} {:02d} {:011.8f}\n"
+                         .format(e[0], e[1], e[2], e[3], e[4], e[5]))
+
+                for sat in self.sat:
+
+                    if np.isnan(peph.pos[sat-1][0:3]).any():
+                        continue
+
+                    clk = 0.999999999999 \
+                        if np.isnan(peph.pos[sat-1][3]) else peph.pos[sat-1][3]
+
+                    fh.write("P{:3s} {:13.6f} {:13.6f} {:13.6f} {:13.6f}\n"
+                             .format(sat2id(sat),
+                                     peph.pos[sat-1][0]*1e-3,
+                                     peph.pos[sat-1][1]*1e-3,
+                                     peph.pos[sat-1][2]*1e-3,
+                                     clk*1e+6))
+
+            # Terminate file
+            #
+            fh.write('EOF\n')
+
     def interppol(self, x, y, n):
         for j in range(1, n):
             for i in range(n-j):
@@ -275,12 +362,12 @@ def pephpos(self, time, sat, nav, vare=False, varc=False):
             dts[0] = c[1]
             if dts[0] != 0.0:
                 std = p_[1].std[sat-1, 3]*rCST.CLIGHT-EXTERR_CLK*t[1]
-        elif c[0] != 0.0 and c[1] != 0.0:
+        elif c[0] != np.nan and c[1] != np.nan:
             dts[0] = (c[1]*t[0]-c[0]*t[1])/(t[0]-t[1])
             i = 0 if t[0] < -t[1] else 1
             std = p_[i].std[sat-1, 3]+EXTERR_CLK*abs(t[i])
         else:
-            dts[0] = 0.0
+            dts[0] = np.nan
 
         if varc:
             varc = std**2
@@ -686,6 +773,7 @@ def antModelTx(nav, e, sigs, sat, time, rs, sig0=None):
         range correction for each specified signal
     """
 
+    sys, _ = sat2prn(sat)
     # Select satellite antenna
     #
     ant = searchpcv(nav.sat_ant, sat, time)
@@ -706,8 +794,10 @@ def antModelTx(nav, e, sigs, sat, time, rs, sig0=None):
     #
     off0 = np.zeros(3)
     if sig0 is not None:
-
-        freq = [s.frequency() for s in sig0]
+        if sys == uGNSS.GLO:
+            freq = [s.frequency(nav.glo_ch[sat]) for s in sig0]
+        else:
+            freq = [s.frequency() for s in sig0]
         fac0 = [1.0 for s in sig0]
 
         if len(freq) == 2:
diff --git a/src/cssrlib/ppp.py b/src/cssrlib/ppp.py
index 47bb600..bb5e027 100644
--- a/src/cssrlib/ppp.py
+++ b/src/cssrlib/ppp.py
@@ -2,9 +2,11 @@
 module for PPP processing
 """
 import cssrlib.gnss as gn
-from cssrlib.peph import gpst2utc
-import numpy as np
+from cssrlib.peph import gpst2utc, time2epoch
+from enum import IntEnum
 from math import sin, cos, atan2, asin
+import numpy as np
+from pysolid.solid import solid_grid
 
 
 def nut_iau1980(t_, f):
@@ -308,6 +310,16 @@ def windupcorr(time, rs, vs, rr, phw, full=False):
     return phw
 
 
+class uTideModel(IntEnum):
+    """
+    Enumeration for Earth tide model selection
+    """
+
+    NONE = -1
+    SIMPLE = 0
+    IERS2010 = 1
+
+
 def tide_pl(eu, rp, GMp, pos):
     """ pole tide correction """
     H3 = 0.293
@@ -372,6 +384,19 @@ def tidedisp(tutc, pos, erpv=None):
     return dr
 
 
+def tidedispIERS2010(tutc, pos, erpv=None):
+    """
+    Wrapper for solid_grid() method of PySolid module to compute Earth tide
+    displacement corrections according to the IERS2010 conventions
+    """
+    e = time2epoch(tutc)
+    disp_e, disp_n, disp_u = solid_grid(e[0], e[1], e[2], e[3], e[4], int(e[5]),
+                                        np.rad2deg(pos[0]), 0, 1,
+                                        np.rad2deg(pos[1]), 0, 1)
+    E = gn.enu2xyz(pos)
+    return E@np.array([disp_e[0, 0], disp_n[0, 0], disp_u[0, 0]])
+
+
 if __name__ == '__main__':
     import matplotlib.pyplot as plt
     from cssrlib.ephemeris import findeph, eph2pos
diff --git a/src/cssrlib/pppigs.py b/src/cssrlib/pppigs.py
deleted file mode 100644
index 7bbcf3c..0000000
--- a/src/cssrlib/pppigs.py
+++ /dev/null
@@ -1,926 +0,0 @@
-"""
-module for standard PPP positioning
-"""
-
-import numpy as np
-
-import cssrlib.gnss as gn
-from cssrlib.ephemeris import satposs
-from cssrlib.gnss import sat2id, sat2prn, uTYP
-from cssrlib.gnss import uTropoModel
-from cssrlib.gnss import time2str
-from cssrlib.ppp import tidedisp, shapiro, windupcorr
-from cssrlib.peph import antModelRx, antModelTx
-from cssrlib.rtk import IB, ddcov, resamb_lambda, valpos, holdamb, initx
-
-
-def IT(na):
-    """ return index of zenith tropospheric delay estimate """
-    return na-gn.uGNSS.MAXSAT-1
-
-
-def II(s, na):
-    """ return index of slant ionospheric delay estimate """
-    return na-gn.uGNSS.MAXSAT+s-1
-
-
-def varerr(nav, el, f):
-    """ variation of measurement """
-    s_el = np.sin(el) if el > np.deg2rad(0.1) else np.sin(np.deg2rad(0.1))
-    fact = nav.eratio[f-nav.nf] if f >= nav.nf else 1
-    a = fact*nav.err[1]
-    b = fact*nav.err[2]
-    return (a**2+(b/s_el)**2)
-
-
-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)
-
-
-def rtkinit(nav, pos0=np.zeros(3), logfile=None):
-    """ initialize variables for PPP """
-
-    # Number of frequencies (actually signals!)
-    #
-    nav.ephopt = 4  # IGS
-
-    # Select tropospheric model
-    #
-    nav.trpModel = uTropoModel.SAAST
-
-    # 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
-
-    # State vector dimensions (including slant iono delay and ambiguities)
-    #
-    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.phw = np.zeros(gn.uGNSS.MAXSAT)
-    nav.el = np.zeros(gn.uGNSS.MAXSAT)
-
-    # Parameters for PPP
-    #
-    # Observation noise parameters
-    #
-    nav.eratio = np.ones(nav.nf)*100  # [-] factor
-    nav.err = [0, 0.000, 0.003]       # [m] sigma
-
-    # Initial sigma for state covariance
-    #
-    nav.sig_p0 = 100.0   # [m]
-    nav.sig_v0 = 1.0     # [m/s]
-    nav.sig_ztd0 = 0.25  # [m]
-    nav.sig_ion0 = 10.0  # [m]
-    nav.sig_n0 = 30.0    # [cyc]
-
-    # 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)]
-    nav.sig_qion = 10.0/np.sqrt(1)        # [m/s/sqrt(s)]
-
-    # Processing options
-    #
-    nav.tidecorr = True
-    nav.useBiases = True
-    nav.useRxPco = True
-
-    nav.thresar = 3.0  # AR acceptance threshold
-    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
-    # 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[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[4:4+gn.uGNSS.MAXSAT] = nav.sig_qion**2
-
-    # Logging level
-    #
-    nav.fout = None
-    if logfile is None:
-        nav.monlevel = 0
-    else:
-        nav.fout = open(logfile, 'w')
-
-
-def sysidx(satlist, sys_ref):
-    """ return index of satellites with sys=sys_ref """
-    idx = []
-    for k, sat in enumerate(satlist):
-        sys, _ = gn.sat2prn(sat)
-        if sys == sys_ref:
-            idx.append(k)
-    return idx
-
-
-def udstate(nav, obs):
-    """ time propagation of states and initialize """
-
-    tt = gn.timediff(obs.t, nav.t)
-
-    ns = len(obs.sat)
-    sys = []
-    sat = obs.sat
-    for sat_i in obs.sat:
-        sys_i, _ = gn.sat2prn(sat_i)
-        sys.append(sys_i)
-
-    # pos,vel,ztd,ion,amb
-    #
-    nx = nav.nx
-    Phi = np.eye(nx)
-    ni = nav.na-gn.uGNSS.MAXSAT
-    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.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
-
-    # Update Kalman filter state elements
-    #
-    for f in range(nav.nf):
-
-        # Reset phase-ambiguity if instantaneous AR
-        # or expire obs outage counter
-        #
-        for i in range(gn.uGNSS.MAXSAT):
-
-            sat_ = i+1
-            sys_i, _ = gn.sat2prn(sat_)
-
-            nav.outc[i, f] += 1
-            reset = (nav.outc[i, f] > nav.maxout or np.any(nav.edt[i, :] > 0))
-            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
-
-                if nav.monlevel > 0:
-                    nav.fout.write("{}  {} - reset ambiguity  {}\n"
-                                   .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:
-                    nav.fout.write("{}  {} - reset ionosphere\n"
-                                   .format(time2str(obs.t), sat2id(sat_)))
-
-        # Ambiguity
-        #
-        bias = np.zeros(ns)
-        ion = np.zeros(ns)
-
-        """
-        offset = 0
-        na = 0
-        """
-        for i in range(ns):
-
-            # Do not initialize invalid observations
-            #
-            if np.any(nav.edt[sat[i]-1, :] > 0):
-                continue
-
-            # 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/lam + \
-                2.0*ion[i]/lam*(sig1.frequency()/sig.frequency())**2
-
-            """
-            amb = nav.x[IB(sat[i], f, nav.na)]
-            if amb != 0.0:
-                offset += bias[i] - amb
-                na += 1
-            """
-        """
-        # 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
-        #
-        for i in range(ns):
-
-            sys_i, _ = sat2prn(sat[i])
-
-            j = IB(sat[i], f, nav.na)
-            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]
-                    nav.fout.write("{}  {} - init  ambiguity  {} {:12.3f}\n"
-                                   .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:
-                    nav.fout.write("{}  {} - init  ionosphere      {:12.3f}\n"
-                                   .format(time2str(obs.t), sat2id(sat[i]), ion[i]))
-
-    return 0
-
-
-def zdres(nav, obs, bsx, rs, vs, dts, 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))
-    el = np.zeros(n)
-    e = np.zeros((n, 3))
-    rr_ = rr.copy()
-
-    # Solid Earth tide corrections
-    #
-    # TODO: add solid earth tide displacements
-    #
-    if nav.tidecorr:
-        pos = gn.ecef2pos(rr_)
-        disp = tidedisp(gn.gpst2utc(obs.t), pos)
-    else:
-        disp = np.zeros(3)
-    rr_ += disp
-
-    # Geodetic position
-    #
-    pos = gn.ecef2pos(rr_)
-
-    # Zenith tropospheric dry and wet delays at user position
-    #
-    trop_hs, trop_wet, _ = gn.tropmodel(obs.t, pos, model=nav.trpModel)
-
-    cpc = np.zeros((n, nf))
-    prc = np.zeros((n, nf))
-
-    for i in range(n):
-
-        sat = obs.sat[i]
-        sys, _ = gn.sat2prn(sat)
-
-        # Skip edited observations
-        #
-        if np.any(nav.edt[sat-1, :] > 0):
-            continue
-
-        # Pseudorange, carrier-phase and C/N0 signals
-        #
-        sigsPR = obs.sig[sys][gn.uTYP.C]
-        sigsCP = obs.sig[sys][gn.uTYP.L]
-
-        # Wavelength
-        #
-        lam = np.array([s.wavelength() for s in sigsCP])
-
-        # Code and phase signal bias, converted from [ns] to [m]
-        #
-        if nav.useBiases:
-            cbias = np.array([bsx.getosb(sat, obs.t, s)*ns2m for s in sigsPR])
-            pbias = np.array([bsx.getosb(sat, obs.t, s)*ns2m for s in sigsCP])
-        else:
-            cbias = np.array([0.0 for s in sigsPR])
-            pbias = np.array([0.0 for s in sigsCP])
-
-        # Check for invalid biases
-        #
-        if np.isnan(cbias).any() or np.isnan(pbias).any():
-            if nav.monlevel > 0:
-                print("skip invalid cbias/pbias for sat={:d}".format(sat))
-            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
-
-        # Shapiro relativistic effect
-        #
-        relatv = shapiro(rs[i, :], rr_)
-
-        # Tropospheric delay mapping functions
-        #
-        mapfh, mapfw = gn.tropmapf(obs.t, pos, el[i], model=nav.trpModel)
-
-        # Tropospheric delay
-        #
-        trop = mapfh*trop_hs + mapfw*trop_wet
-
-        # Phase wind-up effect
-        #
-        nav.phw[sat-1] = windupcorr(obs.t, rs[i, :], vs[i, :], rr_,
-                                    nav.phw[sat-1], full=True)
-
-        # cycle -> m
-        phw = lam*nav.phw[sat-1]
-
-        # Receiver/satellite antenna offset
-        #
-        if nav.useRxPco:
-            antrPR = antModelRx(nav, pos, e[i, :], sigsPR)
-            antrCP = antModelRx(nav, pos, e[i, :], sigsCP)
-        else:
-            antrPR = np.array([0.0 for s in sigsPR])
-            antrCP = np.array([0.0 for s in sigsCP])
-
-        antsPR = antModelTx(nav, e[i, :], sigsPR, sat, obs.t, rs[i, :])
-        antsCP = antModelTx(nav, e[i, :], sigsCP, sat, obs.t, rs[i, :])
-
-        if antrPR is None or antrCP is None or \
-           antsPR is None or antsCP is None:
-            continue
-
-        # Range correction
-        #
-        prc[i, :] = trop + antrPR + antsPR + cbias
-        cpc[i, :] = trop + antrCP + antsCP + pbias + phw
-
-        r += relatv - _c*dts[i]
-
-        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 sdres(nav, obs, 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
-    """
-
-    nf = nav.nf  # number of frequencies (or signals)
-    ns = len(el)  # number of satellites
-    nc = len(obs.sig.keys())  # number of constellations
-
-    nb = np.zeros(2*nc*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])
-
-    # Loop over constellations
-    #
-    for sys in obs.sig.keys():
-
-        # 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
-            #
-            if f < nf:  # carrier
-                sig = obs.sig[sys][uTYP.L][f]
-                mu = -(freq0/sig.frequency())**2
-            else:  # code
-                sig = obs.sig[sys][uTYP.C][f % nf]
-                mu = +(freq0/sig.frequency())**2
-
-            # Select satellites from one constellation only
-            #
-            idx = sysidx(sat, sys)
-
-            # Select reference satellite with highest elevation
-            #
-            if len(idx) > 0:
-                i = idx[np.argmax(el[idx])]
-
-            # Loop over satellites
-            #
-            for j in idx:
-
-                # Skip edited observations
-                #
-                if np.any(nav.edt[sat[j]-1, :] > 0):
-                    continue
-
-                # Skip invalid measurements
-                # NOTE: this additional test is included here, since biases or
-                #       antenna offsets may not be available and this zdres()
-                #       returns zero observation residuals!
-                #
-                if y[i, f] == 0.0 or y[j, f] == 0.0:
-                    continue
-
-                # Skip reference satellite i
-                #
-                if i == j:
-                    continue
-
-                #  Single-difference measurement
-                #
-                v[nv] = y[i, f] - y[j, f]
-
-                # SD line-of-sight vectors
-                #
-                H[nv, 0:3] = -e[i, :] + e[j, :]
-
-                # SD troposphere
-                #
-                _, mapfwi = gn.tropmapf(obs.t, pos, el[i], model=nav.trpModel)
-                _, mapfwj = gn.tropmapf(obs.t, pos, el[j], model=nav.trpModel)
-
-                idx_i = IT(nav.na)
-                H[nv, idx_i] = mapfwi - mapfwj
-                v[nv] -= (mapfwi - mapfwj)*x[idx_i]
-
-                if nav.monlevel > 2:
-                    nav.fout.write("{}         ztd      ({:3d},{:3d}) {:10.3f} {:10.3f} {:10.3f}\n"
-                                   .format(time2str(obs.t),
-                                           idx_i, idx_i, (mapfwi - mapfwj),
-                                           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 > 2:
-                    nav.fout.write("{} {}-{} ion {} ({:3d},{:3d}) {:10.3f} {:10.3f} {:10.3f} {:10.3f} {:10.3f}\n"
-                                   .format(time2str(obs.t),
-                                           sat2id(sat[i]), sat2id(sat[j]),
-                                           sig, 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 = sig.wavelength()
-
-                    H[nv, idx_i] = +lami
-                    H[nv, idx_j] = -lami
-                    v[nv] -= lami*(x[idx_i] - x[idx_j])
-
-                    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 > 2:
-                        nav.fout.write("{} {}-{} amb {} ({:3d},{:3d}) {:10.3f} {:10.3f} {:10.3f} {:10.3f} {:10.3f}\n"
-                                       .format(time2str(obs.t),
-                                               sat2id(sat[i]), sat2id(sat[j]),
-                                               sig, 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
-
-                if nav.monlevel > 1:
-                    nav.fout.write("{} {}-{} res {} ({:3d}) {:10.3f} sig_i {:10.3f} sig_j {:10.3f}\n"
-                                   .format(time2str(obs.t),
-                                           sat2id(sat[i]), sat2id(sat[j]), sig,
-                                           nv, v[nv],
-                                           np.sqrt(Ri[nv]), np.sqrt(Rj[nv])))
-
-                nb[b] += 1  # counter for single-differences per signal
-                nv += 1  # counter for single-difference observations
-
-            b += 1  # counter for signal (pseudrange+carrier-phase)
-
-    v = np.resize(v, nv)
-    H = np.resize(H, (nv, nav.nx))
-    R = ddcov(nb, b, Ri, Rj, nv)
-
-    return v, H, R
-
-
-def kfupdate(x, P, H, v, R):
-    """ Kalman filter measurement update """
-    PHt = P@H.T
-    S = H@PHt+R
-    K = PHt@np.linalg.inv(S)
-    x += K@v
-    P = P - K@H@P
-
-    return x, P, S
-
-
-def qcedit(nav, obs, rs, dts, svh):
-    """ Coarse quality control and editing of observations """
-
-    # Predicted position at next epoch
-    #
-    tt = gn.timediff(obs.t, nav.t)
-    rr_ = nav.x[0:3].copy()
-    if nav.pmode > 0:
-        rr_ += nav.x[3:6]*tt
-
-    # Solid Earth tide corrections
-    #
-    # TODO: add solid earth tide displacements
-    #
-    if nav.tidecorr:
-        pos = gn.ecef2pos(rr_)
-        disp = tidedisp(gn.gpst2utc(obs.t), pos)
-    else:
-        disp = np.zeros(3)
-    rr_ += disp
-
-    # Geodetic position
-    #
-    pos = gn.ecef2pos(rr_)
-
-    # Total number of satellites
-    #
-    ns = gn.uGNSS.MAXSAT
-
-    # Reset previous editing results
-    #
-    nav.edt = np.zeros((ns, nav.nf), dtype=int)
-
-    # Loop over all satellites
-    #
-    sat = []
-    for i in range(ns):
-
-        sat_i = i+1
-        sys_i, _ = sat2prn(sat_i)
-
-        if sat_i not in obs.sat:
-            nav.edt[i, :] = 1
-            continue
-
-        # Check satellite exclusion
-        #
-        if sat_i in nav.excl_sat:
-            nav.edt[i, :] = 1
-            if nav.monlevel > 0:
-                nav.fout.write("{}  {} - edit - satellite excluded\n"
-                               .format(time2str(obs.t), sat2id(sat_i)))
-            continue
-
-        j = np.where(obs.sat == sat_i)[0][0]
-
-        # Check for valid orbit and clock offset
-        #
-        if np.isnan(rs[j, :]).any() or np.isnan(dts[j]):
-            nav.edt[i, :] = 1
-            if nav.monlevel > 0:
-                nav.fout.write("{}  {} - edit - invalid eph\n"
-                               .format(time2str(obs.t), sat2id(sat_i)))
-            continue
-
-        # Check satellite health
-        #
-        if svh[j] > 0:
-            nav.edt[i, :] = 1
-            if nav.monlevel > 0:
-                nav.fout.write("{}  {} - edit - satellite unhealthy\n"
-                               .format(time2str(obs.t), sat2id(sat_i)))
-            continue
-
-        # Check elevation angle
-        #
-        _, e = gn.geodist(rs[j, :], rr_)
-        _, el = gn.satazel(pos, e)
-        if el < nav.elmin:
-            nav.edt[i][:] = 1
-            if nav.monlevel > 0:
-                nav.fout.write("{}  {} - edit - low elevation {:5.1f} deg\n"
-                               .format(time2str(obs.t), sat2id(sat_i),
-                                       np.rad2deg(el)))
-            continue
-
-        # Pseudorange, carrier-phase and C/N0 signals
-        #
-        sigsPR = obs.sig[sys_i][gn.uTYP.C]
-        sigsCP = obs.sig[sys_i][gn.uTYP.L]
-        sigsCN = obs.sig[sys_i][gn.uTYP.S]
-
-        # Loop over signals
-        #
-        for f in range(nav.nf):
-
-            # Cycle  slip check by LLI
-            #
-            if obs.lli[j, f] == 1:
-                nav.edt[i, f] = 1
-                if nav.monlevel > 0:
-                    nav.fout.write("{}  {} - edit {:4s} - LLI\n"
-                                   .format(time2str(obs.t), sat2id(sat_i),
-                                           sigsCP[f].str()))
-                continue
-
-            # Check for measurement consistency
-            #
-            if obs.P[j, f] == 0.0:
-                nav.edt[i, f] = 1
-                if nav.monlevel > 0:
-                    nav.fout.write("{}  {} - edit {:4s} - invalid PR obs\n"
-                                   .format(time2str(obs.t), sat2id(sat_i),
-                                           sigsPR[f].str()))
-                continue
-
-            if obs.L[j, f] == 0.0:
-                nav.edt[i, f] = 1
-                if nav.monlevel > 0:
-                    nav.fout.write("{}  {} - edit {:4s} - invalid CP obs\n"
-                                   .format(time2str(obs.t), sat2id(sat_i),
-                                           sigsCP[f].str()))
-                continue
-
-            # Check C/N0
-            #
-            cnr_min = nav.cnr_min_gpy if sigsCN[f].isGPS_PY() else nav.cnr_min
-            if obs.S[j, f] < cnr_min:
-                nav.edt[i, f] = 1
-                if nav.monlevel > 0:
-                    nav.fout.write("{}  {} - edit {:4s} - low C/N0 {:4.1f} dB-Hz\n"
-                                   .format(time2str(obs.t), sat2id(sat_i),
-                                           sigsCN[f].str(), obs.S[j, f]))
-                continue
-
-        # Store satellite which have passed all tests
-        #
-        if np.any(nav.edt[i, :] > 0):
-            continue
-
-        sat.append(sat_i)
-
-    return np.array(sat, dtype=int)
-
-
-def ppppos(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 for all satellites
-    # in RINEX observations
-    #
-    rs, vs, dts, svh, nsat = satposs(obs, nav, cs=None, orb=orb)
-
-    if nsat < 6:
-        print(" too few satellites < 6: nsat={:d}".format(nsat))
-        return
-
-    # Editing of observations
-    #
-    sat_ed = qcedit(nav, obs, rs, dts, svh)
-
-    # Kalman filter time propagation, initialization of ambiguities and iono
-    #
-    udstate(nav, obs)
-
-    xa = np.zeros(nav.nx)
-    xp = nav.x.copy()
-
-    # Non-differential residuals
-    #
-    yu, eu, elu = zdres(nav, obs, bsx, rs, vs, dts, xp[0:3])
-
-    # Select satellites having passed quality control
-    #
-    # index of valid sats in obs.sat
-    iu = np.where(np.isin(obs.sat, sat_ed))[0]
-    sat = obs.sat[iu]
-    y = yu[iu, :]
-    e = eu[iu, :]
-    el = elu[iu]
-
-    # Store reduced satellite list
-    # NOTE: where are working on a reduced list of observations from here on
-    #
-    nav.sat = sat
-    nav.el[sat-1] = el  # needed in rtk.ddidx()
-    nav.y = y
-    ns = len(sat)
-
-    # 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
-
-    # SD residuals
-    #
-    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 residuals after measurement update
-    #
-    yu, eu, elu = zdres(nav, obs, bsx, rs, vs, dts, xp[0:3])
-    y = yu[iu, :]
-    e = eu[iu, :]
-    ny = y.shape[0]
-    if ny < 6:
-        return -1
-
-    # 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
-        nav.ns = 0
-        for i in range(ns):
-            j = sat[i]-1
-            for f in range(nav.nf):
-                if nav.vsat[j, f] == 0:
-                    continue
-                nav.outc[j, f] = 0
-                if f == 0:
-                    nav.ns += 1
-    else:
-        nav.smode = 0
-
-    nav.smode = 5  # 4: fixed ambiguities, 5: float ambiguities
-
-    if nav.armode > 0:
-        nb, xa = resamb_lambda(nav, sat)
-        if nb > 0:
-            # Use position with fixed ambiguities xa
-            yu, eu, elu = zdres(nav, obs, bsx, rs, vs, dts, xa[0:3])
-            y = yu[iu, :]
-            e = eu[iu, :]
-            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 408aeef..a0c34f7 100644
--- a/src/cssrlib/ppprtk.py
+++ b/src/cssrlib/ppprtk.py
@@ -3,414 +3,22 @@
 """
 
 import numpy as np
-import cssrlib.gnss as gn
-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
+from cssrlib.pppssr import pppos
 
-MAXITR = 10
-ELMIN = 10
-NX = 4
 
+class ppprtkpos(pppos):
+    """ class for PPP-RTK processing """
 
-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 __init__(self, nav, pos0=np.zeros(3), logfile=None):
+        """ initialize variables for PPP-RTK """
 
+        # trop, iono from cssr
+        # phase windup model is local/regional
+        super().__init__(nav=nav, pos0=pos0, logfile=logfile,
+                         trop_opt=2, iono_opt=2, phw_opt=2)
 
-def rtkinit(nav, pos0=np.zeros(3), logfile=None):
-    """ initialize variables for RTK """
-
-    nav.pmode = 1  # 0:static, 1:kinematic
-    nav.monlevel = 1
-
-    nav.na = 3 if nav.pmode == 0 else 6
-    nav.nq = 3 if nav.pmode == 0 else 6
-
-    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.phw = np.zeros(gn.uGNSS.MAXSAT)
-    nav.el = np.zeros(gn.uGNSS.MAXSAT)
-
-    # parameter for PPP-RTK
-    nav.eratio = [50, 50]
-    nav.err = [0, 0.01, 0.005]/np.sqrt(2)
-    nav.sig_p0 = 30.0
-    nav.sig_v0 = 1.0
-    nav.sig_n0 = 30.0
-    nav.sig_qp = 0.01
-    nav.sig_qv = 1.0
-    nav.tidecorr = True
-    nav.armode = 1  # 1:continuous,2:instantaneous,3:fix-and-hold
-    nav.elmaskar = np.deg2rad(20)  # elevation mask for AR
-    nav.x[0:3] = pos0
-
-    dP = np.diag(nav.P)
-    dP.flags['WRITEABLE'] = True
-    dP[0:3] = nav.sig_p0**2
-    nav.q = np.zeros(nav.nq)
-    if nav.pmode >= 1:  # kinematic
-        nav.x[3:6] = 0.0
-        dP[3:6] = nav.sig_v0**2
-        nav.q[0:3] = nav.sig_qp**2
-        nav.q[3:6] = nav.sig_qv**2
-    else:
-        nav.q[0:3] = nav.sig_qp**2
-
-    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]}
-
-    # Logging level
-    #
-    nav.fout = None
-    if logfile is None:
-        nav.monlevel = 0
-    else:
-        nav.fout = open(logfile, 'w')
-
-    if nav.monlevel >= 2:
-        nav.fout = open(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")
-
-
-def udstate(nav, obs, cs):
-    """ time propagation of states and initialize """
-    tt = gn.timediff(obs.t, nav.t)
-
-    ns = len(obs.sat)
-    sys = []
-    sat = obs.sat
-    for sat_i in obs.sat:
-        sys_i, _ = gn.sat2prn(sat_i)
-        sys.append(sys_i)
-
-    # pos,vel
-    na = nav.na
-    Phi = np.eye(nav.na)
-    if nav.na > 3:
-        nav.x[0:3] += nav.x[3:6]*tt
-        Phi[0:3, 3:6] = np.eye(3)*tt
-    nav.P[0:na, 0:na] = Phi@nav.P[0:na, 0:na]@Phi.T
-    dP = np.diag(nav.P)
-    dP.flags['WRITEABLE'] = True
-    dP[0:nav.nq] += nav.q[0:nav.nq]*tt
-
-    # bias
-    for f in range(nav.nf):
-        # reset phase-bias 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 obs.sig.keys():
-                continue
-            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
-        for i in range(ns):
-            if sys[i] not in obs.sig.keys():
-                continue
-            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
-        for i in range(ns):
-            if sat[i] in cs.sat_n:
-                continue
-            initx(nav, 0.0, 0.0, IB(sat[i], f, nav.na))
-        # bias
-        bias = np.zeros(ns)
-        offset = 0
-        na = 0
-        for i in range(ns):
-            if sys[i] not in obs.sig.keys():
-                continue
-
-            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/lam
-            amb = nav.x[IB(sat[i], f, nav.na)]
-            if amb != 0.0:
-                offset += bias[i]-amb
-                na += 1
-        # 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
-        for i in range(ns):
-            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)
-    return 0
-
-
-def zdres(nav, obs, rs, vs, dts, svh, rr, cs):
-    """ non-differential residual """
-    _c = gn.rCST.CLIGHT
-    nf = nav.nf
-    n = len(obs.P)
-    y = np.zeros((n, nf*2))
-    el = np.zeros(n)
-    e = np.zeros((n, 3))
-    rr_ = rr.copy()
-    if nav.tidecorr:
-        pos = gn.ecef2pos(rr_)
-        disp = tidedisp(gn.gpst2utc(obs.t), pos)
-        rr_ += disp
-    pos = gn.ecef2pos(rr_)
-
-    inet = cs.find_grid_index(pos)
-    dlat, dlon = cs.get_dpos(pos)
-    trph, trpw = cs.get_trop(dlat, dlon)
-
-    trop_hs, trop_wet, _ = tropmodel(obs.t, pos)
-    trop_hs0, trop_wet0, _ = tropmodel(obs.t, [pos[0], pos[1], 0])
-    r_hs = trop_hs/trop_hs0
-    r_wet = trop_wet/trop_wet0
-
-    stec = cs.get_stec(dlat, dlon)
-
-    cs.cpc = np.zeros((n, nf))
-    cs.prc = np.zeros((n, nf))
-    cs.osr = np.zeros((n, 4*nf+5))
-
-    for i in range(n):
-        sat = obs.sat[i]
-        sys, _ = gn.sat2prn(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
-        idx_n = np.where(np.array(cs.sat_n) == sat)[0][0]  # global
-        idx_l = np.where(np.array(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:
-            continue
-        # check for measurement consistency
-        flg_m = True
-        for f in range(nav.nf):
-            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
-        if flg_m is False:
-            continue
-
-        r, e[i, :] = gn.geodist(rs[i, :], rr_)
-        _, el[i] = gn.satazel(pos, e[i, :])
-        if el[i] < nav.elmin:
-            continue
-
-        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)
-        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
-        relatv = shapiro(rs[i, :], rr_)
-
-        # 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]
-
-        antrPR = antModelRx(nav, pos, e[i, :], sigsPR)
-        antrCP = antModelRx(nav, pos, e[i, :], sigsCP)
-
-        # range correction
-        prc_c = trop+relatv
-
-        # prc_c += nav.dorb[sat]-nav.dclk[sat]
-        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],
-                        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):
-            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
-
-
-def kfupdate(x, P, H, v, R):
-    """ Kalman filter measurement update """
-    PHt = P@H.T
-    S = H@PHt+R
-    K = PHt@np.linalg.inv(S)
-    x += K@v
-    P = P - K@H@P
-
-    return x, P, S
-
-
-def ppprtkpos(nav, obs, cs):
-    """ PPP-RTK positioning """
-    rs, vs, dts, svh, _ = satposs(obs, nav, cs)
-    # Kalman filter time propagation
-    udstate(nav, obs, cs)
-    # 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)
-
-    iu = np.where(elu >= nav.elmin)[0]
-    sat = obs.sat[iu]
-    y = yu[iu, :]
-    e = eu[iu, :]
-    el = elu[iu]
-    nav.sat = sat
-    nav.el[sat-1] = el
-    nav.y = y
-    ns = len(sat)
-
-    if nav.monlevel > 1:
-        logmon(nav, obs.t, sat, cs, iu)
-
-    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, 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)
-    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, xp, y, e, sat, el)
-    if valpos(nav, v, R):
-        nav.x = xp
-        nav.P = Pp
-        nav.ns = 0
-        for i in range(ns):
-            j = sat[i]-1
-            for f in range(nav.nf):
-                if nav.vsat[j, f] == 0:
-                    continue
-                nav.outc[j, f] = 0
-                if f == 0:
-                    nav.ns += 1
-    else:
-        nav.smode = 0
-
-    nb, xa = resamb_lambda(nav, sat)
-    nav.smode = 5  # float
-    if nb > 0:
-        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, 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
-            nav.smode = 4  # fix
-    nav.t = obs.t
-    return 0
+        self.nav.eratio = np.ones(self.nav.nf)*50  # [-] factor
+        self.nav.err = [0, 0.01, 0.005]/np.sqrt(2)  # [m] sigma
+        self.nav.sig_p0 = 30.0  # [m]
+        self.nav.thresar = 2.0  # AR acceptance threshold
+        self.nav.armode = 1     # AR is enabled
diff --git a/src/cssrlib/pppssr.py b/src/cssrlib/pppssr.py
index b1b1c3d..718f084 100644
--- a/src/cssrlib/pppssr.py
+++ b/src/cssrlib/pppssr.py
@@ -4,1042 +4,1380 @@
 
 import numpy as np
 
-import cssrlib.gnss as gn
 from cssrlib.ephemeris import satposs
-from cssrlib.gnss import sat2id, sat2prn, rSigRnx, uTYP, uGNSS
-from cssrlib.gnss import uTropoModel
-from cssrlib.gnss import time2str
-from cssrlib.ppp import tidedisp, shapiro, windupcorr
+from cssrlib.gnss import sat2id, sat2prn, rSigRnx, uTYP, uGNSS, rCST
+from cssrlib.gnss import uTropoModel, ecef2pos, tropmodel, geodist, satazel
+from cssrlib.gnss import time2str, timediff, gpst2utc, tropmapf
+from cssrlib.ppp import tidedisp, tidedispIERS2010, uTideModel
+from cssrlib.ppp import shapiro, windupcorr
 from cssrlib.peph import antModelRx, antModelTx
-from cssrlib.rtk import IB, ddcov, resamb_lambda, valpos, holdamb, initx
-from cssrlib.cssrlib import sCType, sSigGPS
+from cssrlib.cssrlib import sCType
 from cssrlib.cssrlib import sCSSRTYPE as sc
+from cssrlib.mlambda import mlambda
 
+# format definition for logging
+fmt_ztd = "{}         ztd      ({:3d},{:3d}) {:10.3f} {:10.3f} {:10.3f}\n"
+fmt_ion = "{} {}-{} ion {} ({:3d},{:3d}) {:10.3f} {:10.3f} {:10.3f} " + \
+    "{:10.3f} {:10.3f}\n"
+fmt_res = "{} {}-{} res {} ({:3d}) {:10.3f} sig_i {:10.3f} sig_j {:10.3f}\n"
+fmt_amb = "{} {}-{} amb {} ({:3d},{:3d}) {:10.3f} {:10.3f} {:10.3f} " + \
+    "{:10.3f} {:10.3f} {:10.3f}\n"
 
-def IT(na):
-    """ return index of zenith tropospheric delay estimate """
-    return na-gn.uGNSS.MAXSAT-1
-
-
-def II(s, na):
-    """ return index of slant ionospheric delay estimate """
-    return na-gn.uGNSS.MAXSAT+s-1
-
-
-def varerr(nav, el, f):
-    """ variation of measurement """
-    s_el = np.sin(el) if el > np.deg2rad(0.1) else np.sin(np.deg2rad(0.1))
-    fact = nav.eratio[f-nav.nf] if f >= nav.nf else 1
-    a = fact*nav.err[1]
-    b = fact*nav.err[2]
-    return (a**2+(b/s_el)**2)
-
-
-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)
-
-
-def rtkinit(nav, pos0=np.zeros(3), logfile=None):
-    """ initialize variables for PPP """
-
-    # Number of frequencies (actually signals!)
-    #
-    nav.ephopt = 2  # SSR-APC
-
-    # Select tropospheric model
-    #
-    nav.trpModel = uTropoModel.SAAST
-
-    # 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
-
-    # State vector dimensions (including slant iono delay and ambiguities)
-    #
-    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.phw = np.zeros(gn.uGNSS.MAXSAT)
-    nav.el = np.zeros(gn.uGNSS.MAXSAT)
-
-    # Parameters for PPP
-    #
-    # Observation noise parameters
-    #
-    nav.eratio = np.ones(nav.nf)*100  # [-] factor
-    nav.err = [0, 0.000, 0.003]       # [m] sigma
-
-    # Initial sigma for state covariance
-    #
-    nav.sig_p0 = 100.0   # [m]
-    nav.sig_v0 = 1.0     # [m/s]
-    nav.sig_ztd0 = 0.25  # [m]
-    nav.sig_ion0 = 10.0  # [m]
-    nav.sig_n0 = 30.0    # [cyc]
-
-    # 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)]
-    nav.sig_qion = 10.0/np.sqrt(1)        # [m/s/sqrt(s)]
-
-    nav.tidecorr = True
-    nav.thresar = 3.0  # AR acceptance threshold
-    nav.armode = 0     # 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
-    # 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[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[4:4+gn.uGNSS.MAXSAT] = nav.sig_qion**2
-
-    # Logging level
-    #
-    nav.fout = None
-    if logfile is None:
-        nav.monlevel = 0
-    else:
-        nav.fout = open(logfile, 'w')
-
-
-def sysidx(satlist, sys_ref):
-    """ return index of satellites with sys=sys_ref """
-    idx = []
-    for k, sat in enumerate(satlist):
-        sys, _ = gn.sat2prn(sat)
-        if sys == sys_ref:
-            idx.append(k)
-    return idx
-
-
-def udstate(nav, obs):
-    """ time propagation of states and initialize """
-
-    tt = gn.timediff(obs.t, nav.t)
-
-    ns = len(obs.sat)
-    sys = []
-    sat = obs.sat
-    for sat_i in obs.sat:
-        sys_i, _ = gn.sat2prn(sat_i)
-        sys.append(sys_i)
-
-    # pos,vel,ztd,ion,amb
-    #
-    nx = nav.nx
-    Phi = np.eye(nx)
-    ni = nav.na-gn.uGNSS.MAXSAT
-    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.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
-
-    # Update Kalman filter state elements
-    #
-    for f in range(nav.nf):
-
-        # Reset phase-ambiguity if instantaneous AR
-        # or expire obs outage counter
-        #
-        for i in range(gn.uGNSS.MAXSAT):
-
-            sat_ = i+1
-            sys_i, _ = gn.sat2prn(sat_)
 
-            nav.outc[i, f] += 1
-            reset = (nav.outc[i, f] > nav.maxout or np.any(nav.edt[i, :] > 0))
-            if sys_i not in obs.sig.keys():
-                continue
+class pppos():
+    """ class for PPP processing """
 
-            # 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
+    nav = None
+    VAR_HOLDAMB = 0.001
 
-                if nav.monlevel > 0:
-                    nav.fout.write("{}  {} - reset ambiguity  {}\n"
-                                   .format(time2str(obs.t), sat2id(sat_),
-                                           obs.sig[sys_i][uTYP.L][f]))
+    def __init__(self, nav, pos0=np.zeros(3),
+                 logfile=None, trop_opt=1, iono_opt=1, phw_opt=1):
+        """ initialize variables for PPP """
 
-            # 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)
+        self.nav = nav
 
-                if nav.monlevel > 0:
-                    nav.fout.write("{}  {} - reset ionosphere\n"
-                                   .format(time2str(obs.t), sat2id(sat_)))
+        # Number of frequencies (actually signals!)
+        #
+        self.nav.ephopt = 2  # SSR-APC
 
-        # Ambiguity
+        # Select tropospheric model
         #
-        bias = np.zeros(ns)
-        ion = np.zeros(ns)
+        self.nav.trpModel = uTropoModel.SAAST
 
-        """
-        offset = 0
-        na = 0
-        """
-        for i in range(ns):
+        # 0: use trop-model, 1: estimate, 2: use cssr correction
+        self.nav.trop_opt = trop_opt
 
-            # Do not initialize invalid observations
-            #
-            if np.any(nav.edt[sat[i]-1, :] > 0):
-                continue
+        # 0: use iono-model, 1: estimate, 2: use cssr correction
+        self.nav.iono_opt = iono_opt
 
-            # Get dual-frequency pseudoranges for this constellation
-            #
-            sig1 = obs.sig[sys[i]][uTYP.C][0]
-            sig2 = obs.sig[sys[i]][uTYP.C][1]
+        # 0: none, 1: full model, 2: local/regional model
+        self.nav.phw_opt = phw_opt
 
-            pr1 = obs.P[i, 0]
-            pr2 = obs.P[i, 1]
+        # Position (+ optional velocity), zenith tropo delay and
+        # slant ionospheric delay states
+        #
+        self.nav.ntrop = (1 if self.nav.trop_opt == 1 else 0)
+        self.nav.niono = (uGNSS.MAXSAT if self.nav.iono_opt == 1 else 0)
 
-            # Skip zero observations
-            #
-            if pr1 == 0.0 or pr2 == 0.0:
-                continue
+        self.nav.na = (3 if self.nav.pmode == 0 else 6)
+        self.nav.nq = (3 if self.nav.pmode == 0 else 6)
 
-            # Get iono delay at frequency of first signal
-            #
-            ion[i] = ionoDelay(sig1, sig2, pr1, pr2)
+        self.nav.na += self.nav.ntrop + self.nav.niono
+        self.nav.nq += self.nav.ntrop + self.nav.niono
 
-            # Get pseudorange and carrier-phase observation of signal f
-            #
-            sig = obs.sig[sys[i]][uTYP.L][f]
-            lam = sig.wavelength()
+        # State vector dimensions (including slant iono delay and ambiguities)
+        #
+        self.nav.nx = self.nav.na+uGNSS.MAXSAT*self.nav.nf
 
-            cp = obs.L[i, f]
-            pr = obs.P[i, f]
-            if cp == 0.0 or pr == 0.0 or lam is None:
-                continue
+        self.nav.x = np.zeros(self.nav.nx)
+        self.nav.P = np.zeros((self.nav.nx, self.nav.nx))
 
-            bias[i] = cp - pr/lam + \
-                2.0*ion[i]/lam*(sig1.frequency()/sig.frequency())**2
+        self.nav.xa = np.zeros(self.nav.na)
+        self.nav.Pa = np.zeros((self.nav.na, self.nav.na))
 
-            """
-            amb = nav.x[IB(sat[i], f, nav.na)]
-            if amb != 0.0:
-                offset += bias[i] - amb
-                na += 1
-            """
-        """
-        # Adjust phase-code coherency
+        self.nav.phw = np.zeros(uGNSS.MAXSAT)
+        self.nav.el = np.zeros(uGNSS.MAXSAT)
+
+        # Parameters for PPP
         #
-        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
-        """
+        # Observation noise parameters
+        #
+        self.nav.eratio = np.ones(self.nav.nf)*100  # [-] factor
+        self.nav.err = [0, 0.000, 0.003]       # [m] sigma
 
-        # Initialize ambiguity
+        # Initial sigma for state covariance
         #
-        for i in range(ns):
+        self.nav.sig_p0 = 100.0   # [m]
+        self.nav.sig_v0 = 1.0     # [m/s]
+        self.nav.sig_ztd0 = 0.25  # [m]
+        self.nav.sig_ion0 = 10.0  # [m]
+        self.nav.sig_n0 = 30.0    # [cyc]
 
-            sys_i, _ = sat2prn(sat[i])
+        # Process noise sigma
+        #
+        if self.nav.pmode == 0:
+            self.nav.sig_qp = 100.0/np.sqrt(1)     # [m/sqrt(s)]
+            self.nav.sig_qv = None
+        else:
+            self.nav.sig_qp = 0.01/np.sqrt(1)      # [m/sqrt(s)]
+            self.nav.sig_qv = 1.0/np.sqrt(1)       # [m/s/sqrt(s)]
+        self.nav.sig_qztd = 0.1/np.sqrt(3600)      # [m/sqrt(s)]
+        self.nav.sig_qion = 10.0/np.sqrt(1)        # [m/s/sqrt(s)]
 
-            j = IB(sat[i], f, nav.na)
-            if bias[i] != 0.0 and nav.x[j] == 0.0:
+        # Processing options
+        #
+        self.nav.tidecorr = uTideModel.IERS2010
+        # self.nav.tidecorr = uTideModel.SIMPLE
+        self.nav.thresar = 3.0  # AR acceptance threshold
+        # 0:float-ppp,1:continuous,2:instantaneous,3:fix-and-hold
+        self.nav.armode = 0
+        self.nav.elmaskar = np.deg2rad(20.0)  # elevation mask for AR
+        self.nav.elmin = np.deg2rad(10.0)
+
+        # Initial state vector
+        #
+        self.nav.x[0:3] = pos0
+        if self.nav.pmode >= 1:  # kinematic
+            self.nav.x[3:6] = 0.0  # velocity
 
-                initx(nav, bias[i], nav.sig_n0**2, j)
+        # Diagonal elements of covariance matrix
+        #
+        dP = np.diag(self.nav.P)
+        dP.flags['WRITEABLE'] = True
+
+        dP[0:3] = self.nav.sig_p0**2
+        # Velocity
+        if self.nav.pmode >= 1:  # kinematic
+            dP[3:6] = self.nav.sig_v0**2
+
+        # Tropo delay
+        if self.nav.trop_opt == 1:  # trop is estimated
+            if self.nav.pmode >= 1:  # kinematic
+                dP[6] = self.nav.sig_ztd0**2
+            else:
+                dP[3] = self.nav.sig_ztd0**2
+
+        # Process noise
+        #
+        self.nav.q = np.zeros(self.nav.nq)
+        self.nav.q[0:3] = self.nav.sig_qp**2
+
+        # Velocity
+        if self.nav.pmode >= 1:  # kinematic
+            self.nav.q[3:6] = self.nav.sig_qv**2
+
+        if self.nav.trop_opt == 1:  # trop is estimated
+            # Tropo delay
+            if self.nav.pmode >= 1:  # kinematic
+                self.nav.q[6] = self.nav.sig_qztd**2
+            else:
+                self.nav.q[3] = self.nav.sig_qztd**2
+
+        if self.nav.iono_opt == 1:  # iono is estimated
+            # Iono delay
+            if self.nav.pmode >= 1:  # kinematic
+                self.nav.q[7:7+uGNSS.MAXSAT] = self.nav.sig_qion**2
+            else:
+                self.nav.q[4:4+uGNSS.MAXSAT] = self.nav.sig_qion**2
+
+        # Logging level
+        #
+        self.nav.fout = None
+        if logfile is None:
+            self.nav.monlevel = 0
+        else:
+            self.nav.fout = open(logfile, 'w')
+
+    def valpos(self, v, R, thres=4.0):
+        """ post-fit residual test """
+        nv = len(v)
+        fact = thres**2
+        for i in range(nv):
+            if v[i]**2 <= fact*R[i, i]:
+                continue
+            if self.nav.monlevel > 1:
+                txt = "{:3d} is large: {:8.4f} ({:8.4f})".format(
+                    i, v[i], R[i, i])
+                if self.nav.fout is None:
+                    print(txt)
+                else:
+                    self.nav.fout.write(txt+"\n")
+        return True
+
+    def initx(self, x0, v0, i):
+        """ initialize x and P for index i """
+        self.nav.x[i] = x0
+        for j in range(self.nav.nx):
+            self.nav.P[j, i] = self.nav.P[i, j] = v0 if i == j else 0
+
+    def IB(self, s, f, na=3):
+        """ return index of phase ambguity """
+        idx = na+uGNSS.MAXSAT*f+s-1
+        return idx
+
+    def II(self, s, na):
+        """ return index of slant ionospheric delay estimate """
+        return na-uGNSS.MAXSAT+s-1
+
+    def IT(self, na):
+        """ return index of zenith tropospheric delay estimate """
+        return na-uGNSS.MAXSAT-1
+
+    def varerr(self, nav, el, f):
+        """ variation of measurement """
+        s_el = np.sin(el) if el > np.deg2rad(0.1) else np.sin(np.deg2rad(0.1))
+        fact = nav.eratio[f-nav.nf] if f >= nav.nf else 1
+        a = fact*nav.err[1]
+        b = fact*nav.err[2]
+        return (a**2+(b/s_el)**2)
+
+    def sysidx(self, satlist, sys_ref):
+        """ return index of satellites with sys=sys_ref """
+        idx = []
+        for k, sat in enumerate(satlist):
+            sys, _ = sat2prn(sat)
+            if sys == sys_ref:
+                idx.append(k)
+        return idx
+
+    def udstate(self, obs):
+        """ time propagation of states and initialize """
+
+        tt = timediff(obs.t, self.nav.t)
+
+        ns = len(obs.sat)
+        sys = []
+        sat = obs.sat
+        for sat_i in obs.sat:
+            sys_i, _ = sat2prn(sat_i)
+            sys.append(sys_i)
+
+        # pos,vel,ztd,ion,amb
+        #
+        nx = self.nav.nx
+        Phi = np.eye(nx)
+        # if self.nav.niono > 0:
+        #    ni = self.nav.na-uGNSS.MAXSAT
+        #    Phi[ni:self.nav.na, ni:self.nav.na] = np.zeros(
+        #        (uGNSS.MAXSAT, uGNSS.MAXSAT))
+        if self.nav.pmode > 0:
+            self.nav.x[0:3] += self.nav.x[3:6]*tt
+            Phi[0:3, 3:6] = np.eye(3)*tt
+        self.nav.P[0:nx, 0:nx] = Phi@self.nav.P[0:nx, 0:nx]@Phi.T
+
+        # Process noise
+        #
+        dP = np.diag(self.nav.P)
+        dP.flags['WRITEABLE'] = True
+        dP[0:self.nav.nq] += self.nav.q[0:self.nav.nq]*tt
 
-                if nav.monlevel > 0:
-                    sig = obs.sig[sys_i][uTYP.L][f]
-                    nav.fout.write("{}  {} - init  ambiguity  {} {:12.3f}\n"
-                                   .format(time2str(obs.t), sat2id(sat[i]), sig, bias[i]))
+        # Update Kalman filter state elements
+        #
+        for f in range(self.nav.nf):
 
-            j = II(sat[i], nav.na)
-            if ion[i] != 0 and nav.x[j] == 0.0:
+            # Reset phase-ambiguity if instantaneous AR
+            # or expire obs outage counter
+            #
+            for i in range(uGNSS.MAXSAT):
 
-                initx(nav, ion[i], nav.sig_ion0**2, j)
+                sat_ = i+1
+                sys_i, _ = sat2prn(sat_)
 
-                if nav.monlevel > 0:
-                    nav.fout.write("{}  {} - init  ionosphere      {:12.3f}\n"
-                                   .format(time2str(obs.t), sat2id(sat[i]), ion[i]))
+                self.nav.outc[i, f] += 1
+                reset = (self.nav.outc[i, f] >
+                         self.nav.maxout or np.any(self.nav.edt[i, :] > 0))
+                if sys_i not in obs.sig.keys():
+                    continue
 
-    return 0
+                # Reset ambiguity estimate
+                #
+                j = self.IB(sat_, f, self.nav.na)
+                if reset and self.nav.x[j] != 0.0:
+                    self.initx(0.0, 0.0, j)
+                    self.nav.outc[i, f] = 0
+
+                    if self.nav.monlevel > 0:
+                        self.nav.fout.write(
+                            "{}  {} - reset ambiguity  {}\n"
+                            .format(time2str(obs.t), sat2id(sat_),
+                                    obs.sig[sys_i][uTYP.L][f]))
+
+                if self.nav.niono > 0:
+                    # Reset slant ionospheric delay estimate
+                    #
+                    j = self.II(sat_, self.nav.na)
+                    if reset and self.nav.x[j] != 0.0:
+                        self.initx(0.0, 0.0, j)
+
+                        if self.nav.monlevel > 0:
+                            self.nav.fout.write("{}  {} - reset ionosphere\n"
+                                                .format(time2str(obs.t),
+                                                        sat2id(sat_)))
+
+            # Ambiguity
+            #
+            bias = np.zeros(ns)
+            ion = np.zeros(ns)
+            f1 = 0
 
+            """
+            offset = 0
+            na = 0
+            """
+            for i in range(ns):
 
-def find_corr_idx(cs, nf, ctype, sigref, sat):
-    nsig = 0
-    kidx = [-1]*nf
-    idx_n = []
+                # Do not initialize invalid observations
+                #
+                if np.any(self.nav.edt[sat[i]-1, :] > 0):
+                    continue
 
-    sys, _ = sat2prn(sat)
+                if self.nav.niono > 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
+
+                    if sys[i] == uGNSS.GLO:
+                        if sat[i] not in self.nav.glo_ch:
+                            print("glonass channed not found: {:d}"
+                                  .format(sat[i]))
+                            continue
+                        f1 = sig1.frequency(self.nav.glo_ch[sat[i]])
+                        f2 = sig2.frequency(self.nav.glo_ch[sat[i]])
+                    else:
+                        f1 = sig1.frequency()
+                        f2 = sig2.frequency()
+
+                    # Get iono delay at frequency of first signal
+                    #
+                    ion[i] = (pr1-pr2)/(1.0-(f1/f2)**2)
+
+                # Get pseudorange and carrier-phase observation of signal f
+                #
+                sig = obs.sig[sys[i]][uTYP.L][f]
 
-    if cs.iodssr_c[ctype] == cs.iodssr:
-        idx_n_ = np.where(np.array(cs.sat_n) == sat)[0]
-    else:  # work-around for mask change case
-        if cs.iodssr_c[ctype] == cs.iodssr_p:
-            idx_n_ = np.where(np.array(cs.sat_n_p) == sat)[0]
-        else:
-            return nsig, idx_n, kidx
-
-    if len(idx_n_) == 0:
-        return nsig, idx_n, kidx
-    idx_n = idx_n_[0]
-
-    if cs.iodssr_c[ctype] == cs.iodssr:
-        sig_n = cs.sig_n[idx_n]
-    elif cs.iodssr_c[ctype] == cs.iodssr_p:
-        sig_n = cs.sig_n_p[idx_n]
-    else:
-        return nsig, idx_n, kidx
-
-    utype = uTYP.C if ctype == sCType.CBIAS else uTYP.L
-
-    for k, sig in enumerate(sig_n):
-        if sig < 0:
-            continue
-        for f in range(nf):
-            if cs.cssrmode == sc.GAL_HAS_SIS and \
-               sys == uGNSS.GPS and sig == sSigGPS.L2P:
-                sig = sSigGPS.L2W  # work-around
-            sig_ = cs.ssig2rsig(sys, utype, sig)
-            if sig_ == sigref[f] or sig_ == sigref[f].toAtt('X'):
-                kidx[f] = k
-                nsig += 1
-
-    return nsig, idx_n, kidx
-
-
-def zdres(nav, obs, cs, bsx, rs, vs, dts, 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))
-    el = np.zeros(n)
-    e = np.zeros((n, 3))
-    rr_ = rr.copy()
-
-    # Solid Earth tide corrections
-    #
-    # TODO: add solid earth tide displacements
-    #
-    if nav.tidecorr:
-        pos = gn.ecef2pos(rr_)
-        disp = tidedisp(gn.gpst2utc(obs.t), pos)
-    else:
-        disp = np.zeros(3)
-    rr_ += disp
-
-    # Geodetic position
-    #
-    pos = gn.ecef2pos(rr_)
-
-    # Zenith tropospheric dry and wet delays at user position
-    #
-    trop_hs, trop_wet, _ = gn.tropmodel(obs.t, pos, model=nav.trpModel)
-
-    cpc = np.zeros((n, nf))
-    prc = np.zeros((n, nf))
-
-    for i in range(n):
-
-        sat = obs.sat[i]
-        sys, _ = gn.sat2prn(sat)
-
-        # Skip edited observations
-        #
-        if np.any(nav.edt[sat-1, :] > 0):
-            continue
+                if sys[i] == uGNSS.GLO:
+                    fi = sig.frequency(self.nav.glo_ch[sat[i]])
+                else:
+                    fi = sig.frequency()
 
-        # Pseudorange, carrier-phase and C/N0 signals
-        #
-        sigsPR = obs.sig[sys][gn.uTYP.C]
-        sigsCP = obs.sig[sys][gn.uTYP.L]
+                lam = rCST.CLIGHT/fi
 
-        # Wavelength
-        #
-        lam = np.array([s.wavelength() for s in sigsCP])
+                cp = obs.L[i, f]
+                pr = obs.P[i, f]
+                if cp == 0.0 or pr == 0.0 or lam is None:
+                    continue
 
-        cbias = np.zeros(nav.nf)
-        pbias = np.zeros(nav.nf)
+                bias[i] = cp - pr/lam + 2.0*ion[i]/lam*(f1/fi)**2
 
-        if nav.ephopt == 4:
+                """
+                amb = nav.x[IB(sat[i], f, nav.na)]
+                if amb != 0.0:
+                    offset += bias[i] - amb
+                    na += 1
+                """
+            """
+            # Adjust phase-code coherency
+            #
+            if na > 0:
+                db = offset/na
+                for i in range(uGNSS.MAXSAT):
+                    if nav.x[IB(i+1, f, nav.na)] != 0.0:
+                        nav.x[IB(i+1, f, nav.na)] += db
+            """
 
-            # Code and phase signal bias, converted from [ns] to [m]
+            # Initialize ambiguity
             #
-            cbias = np.array([bsx.getosb(sat, obs.t, s)*ns2m for s in sigsPR])
-            pbias = np.array([bsx.getosb(sat, obs.t, s)*ns2m for s in sigsCP])
+            for i in range(ns):
 
-        else:  # from CSSR
+                sys_i, _ = sat2prn(sat[i])
 
-            if cs.lc[0].cstat & (1 << sCType.CBIAS) == (1 << sCType.CBIAS):
-                nsig, idx_n, kidx = find_corr_idx(cs, nav.nf, sCType.CBIAS,
-                                                  sigsPR, sat)
+                j = self.IB(sat[i], f, self.nav.na)
+                if bias[i] != 0.0 and self.nav.x[j] == 0.0:
 
-                if cs.cssrmode == sc.BDS_PPP and sys == uGNSS.GPS:
-                    # BDS PPP not including cbias for GPS
-                    cbias = np.zeros(nav.nf)
-                else:
-                    cbias = np.ones(nav.nf)*np.nan
-                if nsig >= nav.nf:
-                    cbias = cs.lc[0].cbias[idx_n][kidx]
-                elif nav.monlevel > 1:
-                    print("skip cbias for sat={:d}".format(sat))
-
-                # - IS-QZSS-MDC-001 sec 5.5.3.3
-                # - HAS SIS ICD sec 7.4, 7.5
-                # - HAS IDD ICD sec 3.3.4
-                if cs.cssrmode in [sc.GAL_HAS_IDD, sc.GAL_HAS_SIS, sc.QZS_MADOCA]:
-                    cbias = -cbias
+                    self.initx(bias[i], self.nav.sig_n0**2, j)
 
-            if cs.lc[0].cstat & (1 << sCType.PBIAS) == (1 << sCType.PBIAS):
-                nsig, idx_n, kidx = find_corr_idx(cs, nav.nf, sCType.PBIAS,
-                                                  sigsCP, sat)
-
-                pbias = np.ones(nav.nf)*np.nan
-                if nsig >= nav.nf:
-                    pbias = cs.lc[0].pbias[idx_n][kidx]
-                    # for Gal HAS (cycle -> m)
-                    if cs.cssrmode == sc.GAL_HAS_SIS:
-                        pbias *= lam
-                elif nav.monlevel > 1:
-                    print("skip pbias for sat={:d}".format(sat))
-
-                # - IS-QZSS-MDC-001 sec 5.5.3.3
-                # - HAS SIS ICD sec 7.4, 7.5
-                if cs.cssrmode in [sc.GAL_HAS_SIS, sc.QZS_MADOCA]:
-                    pbias = -pbias
+                    if self.nav.monlevel > 0:
+                        sig = obs.sig[sys_i][uTYP.L][f]
+                        self.nav.fout.write(
+                            "{}  {} - init  ambiguity  {} {:12.3f}\n"
+                            .format(time2str(obs.t), sat2id(sat[i]),
+                                    sig, bias[i]))
 
-        # Check for invalid biases
-        #
-        if np.isnan(cbias).any() or np.isnan(pbias).any():
-            if nav.monlevel > 3:
-                print("skip invalid cbias/pbias for sat={:d}".format(sat))
-            continue
+                if self.nav.niono > 0:
+                    j = self.II(sat[i], self.nav.na)
+                    if ion[i] != 0 and self.nav.x[j] == 0.0:
 
-        # 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
+                        self.initx(ion[i], self.nav.sig_ion0**2, j)
 
-        # Shapiro relativistic effect
-        #
-        relatv = shapiro(rs[i, :], rr_)
+                        if self.nav.monlevel > 0:
+                            self.nav.fout.write(
+                                "{}  {} - init  ionosphere      {:12.3f}\n"
+                                .format(time2str(obs.t), sat2id(sat[i]),
+                                        ion[i]))
+
+        return 0
+
+    def find_bias(self, cs, sigref, sat, inet=0):
+        """ find satellite signal bias from correction """
+        nf = len(sigref)
+        v = np.zeros(nf)
+
+        if nf == 0:
+            return v
 
-        # Tropospheric delay mapping functions
+        ctype = sigref[0].typ
+        if ctype == uTYP.C:
+            if cs.lc[inet].cbias is None or \
+                    sat not in cs.lc[inet].cbias.keys():
+                return v
+            sigc = cs.lc[inet].cbias[sat]
+        else:
+            if cs.lc[inet].pbias is None or \
+                    sat not in cs.lc[inet].pbias.keys():
+                return v
+            sigc = cs.lc[inet].pbias[sat]
+
+        # work-around for Galileo HAS: L2P -> L2W
+        if cs.cssrmode in [sc.GAL_HAS_SIS, sc.GAL_HAS_IDD]:
+            if ctype == uTYP.C and rSigRnx('GC2P') in sigc.keys():
+                sigc[rSigRnx('GC2W')] = sigc[rSigRnx('GC2P')]
+            if ctype == uTYP.L and rSigRnx('GL2P') in sigc.keys():
+                sigc[rSigRnx('GL2W')] = sigc[rSigRnx('GL2P')]
+
+        for k, sig in enumerate(sigref):
+            if sig in sigc.keys():
+                v[k] = sigc[sig]
+            elif sig.toAtt('X') in sigc.keys():
+                v[k] = sigc[sig.toAtt('X')]
+        return v
+
+    def zdres(self, obs, cs, bsx, rs, vs, dts, rr, rtype=1):
+        """ non-differential residual """
+
+        _c = rCST.CLIGHT
+        ns2m = _c*1e-9
+
+        nf = self.nav.nf
+        n = len(obs.P)
+        y = np.zeros((n, nf*2))
+        el = np.zeros(n)
+        e = np.zeros((n, 3))
+        rr_ = rr.copy()
+
+        # Solid Earth tide corrections
         #
-        mapfh, mapfw = gn.tropmapf(obs.t, pos, el[i], model=nav.trpModel)
+        if self.nav.tidecorr == uTideModel.SIMPLE:
+            pos = ecef2pos(rr_)
+            disp = tidedisp(gpst2utc(obs.t), pos)
+        elif self.nav.tidecorr == uTideModel.IERS2010:
+            pos = ecef2pos(rr_)
+            disp = tidedispIERS2010(gpst2utc(obs.t), pos)
+        else:
+            disp = np.zeros(3)
+        rr_ += disp
 
-        # Tropospheric delay
+        # Geodetic position
         #
-        trop = mapfh*trop_hs + mapfw*trop_wet
+        pos = ecef2pos(rr_)
 
-        # Phase wind-up effect
+        # Zenith tropospheric dry and wet delays at user position
         #
-        nav.phw[sat-1] = windupcorr(obs.t, rs[i, :], vs[i, :], rr_,
-                                    nav.phw[sat-1], full=True)
+        trop_hs, trop_wet, _ = tropmodel(obs.t, pos,
+                                         model=self.nav.trpModel)
 
-        # cycle -> m
-        phw = lam*nav.phw[sat-1]
+        if self.nav.trop_opt == 2 or self.nav.iono_opt == 2:  # from cssr
+            inet = cs.find_grid_index(pos)
+            dlat, dlon = cs.get_dpos(pos)
+        else:
+            inet = -1
 
-        # Select APC reference signals
-        #
-        sig0 = None
-        if cs is not None:
+        if self.nav.trop_opt == 2:  # trop from cssr
+            trph, trpw = cs.get_trop(dlat, dlon)
+            trop_hs0, trop_wet0, _ = tropmodel(obs.t, [pos[0], pos[1], 0],
+                                               model=self.nav.trpModel)
+            r_hs = trop_hs/trop_hs0
+            r_wet = trop_wet/trop_wet0
 
-            if (cs.cssrmode == sc.GAL_HAS_SIS or
-                cs.cssrmode == sc.GAL_HAS_IDD or
-                    cs.cssrmode == sc.QZS_MADOCA):
+        if self.nav.iono_opt == 2:  # iono from cssr
+            stec = cs.get_stec(dlat, dlon)
 
-                if sys == uGNSS.GPS:
-                    sig0 = (rSigRnx("GC1W"), rSigRnx("GC2W"))
-                elif sys == uGNSS.GAL:
-                    sig0 = (rSigRnx("EC1C"), rSigRnx("EC7Q"))
-                elif sys == uGNSS.QZS:
-                    sig0 = (rSigRnx("JC1C"), rSigRnx("JC2S"))
+        cpc = np.zeros((n, nf))
+        prc = np.zeros((n, nf))
 
-            elif cs.cssrmode == sc.BDS_PPP:
+        for i in range(n):
 
-                if sys == uGNSS.GPS:
-                    sig0 = (rSigRnx("GC1W"), rSigRnx("GC2W"))
-                elif sys == uGNSS.BDS:
-                    sig0 = (rSigRnx("CC6I"),)
+            sat = obs.sat[i]
+            sys, _ = sat2prn(sat)
 
-        # Receiver/satellite antenna offset
-        #
-        antrPR = antModelRx(nav, pos, e[i, :], sigsPR)
-        antrCP = antModelRx(nav, pos, e[i, :], sigsCP)
+            # Skip edited observations
+            #
+            if np.any(self.nav.edt[sat-1, :] > 0):
+                continue
 
-        if nav.ephopt == 4:
+            if inet > 0 and sat not in cs.lc[inet].sat_n:
+                continue
 
-            antsPR = antModelTx(nav, e[i, :], sigsPR, sat, obs.t, rs[i, :])
-            antsCP = antModelTx(nav, e[i, :], sigsCP, sat, obs.t, rs[i, :])
+            # Pseudorange, carrier-phase and C/N0 signals
+            #
+            sigsPR = obs.sig[sys][uTYP.C]
+            sigsCP = obs.sig[sys][uTYP.L]
 
-        elif cs is not None and (cs.cssrmode == sc.GAL_HAS_SIS or
-                                 cs.cssrmode == sc.GAL_HAS_IDD or
-                                 cs.cssrmode == sc.QZS_MADOCA or
-                                 cs.cssrmode == sc.BDS_PPP):
+            # Wavelength
+            #
+            if sys == uGNSS.GLO:
+                lam = np.array([s.wavelength(self.nav.glo_ch[sat])
+                                for s in sigsCP])
+                frq = np.array([s.frequency(self.nav.glo_ch[sat])
+                               for s in sigsCP])
+            else:
+                lam = np.array([s.wavelength() for s in sigsCP])
+                frq = np.array([s.frequency() for s in sigsCP])
 
-            antsPR = antModelTx(nav, e[i, :], sigsPR,
-                                sat, obs.t, rs[i, :], sig0)
-            antsCP = antModelTx(nav, e[i, :], sigsCP,
-                                sat, obs.t, rs[i, :], sig0)
+            cbias = np.zeros(self.nav.nf)
+            pbias = np.zeros(self.nav.nf)
 
-        else:
+            if self.nav.ephopt == 4:
 
-            antsPR = [0.0 for _ in sigsPR]
-            antsCP = [0.0 for _ in sigsCP]
+                # Code and phase signal bias, converted from [ns] to [m]
+                # note: IGS uses sign convention different with RTCM
+                cbias = np.array(
+                    [-bsx.getosb(sat, obs.t, s)*ns2m for s in sigsPR])
+                pbias = np.array(
+                    [-bsx.getosb(sat, obs.t, s)*ns2m for s in sigsCP])
 
-        # Check for invalid values
-        #
-        if antrPR is None or antrCP is None or \
-           antsPR is None or antsCP is None:
-            continue
+            elif cs is not None:  # from CSSR
 
-        # Range correction
-        #
-        prc[i, :] = trop + antrPR + antsPR + cbias
-        cpc[i, :] = trop + antrCP + antsCP + pbias + phw
-
-        r += relatv - _c*dts[i]
-
-        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 sdres(nav, obs, 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
-    """
-
-    nf = nav.nf  # number of frequencies (or signals)
-    ns = len(el)  # number of satellites
-    nc = len(obs.sig.keys())  # number of constellations
-
-    nb = np.zeros(2*nc*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])
-
-    # Loop over constellations
-    #
-    for sys in obs.sig.keys():
-
-        # Slant ionospheric delay reference frequency
-        #
-        freq0 = obs.sig[sys][uTYP.L][0].frequency()
+                if cs.lc[0].cstat & (1 << sCType.CBIAS) == (1 << sCType.CBIAS):
+                    cbias = self.find_bias(cs, sigsPR, sat)
 
-        # Loop over twice the number of frequencies
-        #   first for all carrier-phase observations
-        #   second all pseudorange observations
-        #
-        for f in range(0, nf*2):
+                if inet > 0 and cs.lc[inet].cstat & (1 << sCType.CBIAS) == \
+                        (1 << sCType.CBIAS):
+                    cbias += self.find_bias(cs, sigsPR, sat, inet)
 
-            # Select carrier-phase frequency and iono frequency ratio
-            #
-            if f < nf:  # carrier
-                sig = obs.sig[sys][uTYP.L][f]
-                mu = -(freq0/sig.frequency())**2
-            else:  # code
-                sig = obs.sig[sys][uTYP.C][f % nf]
-                mu = +(freq0/sig.frequency())**2
-
-            # Select satellites from one constellation only
-            #
-            idx = sysidx(sat, sys)
+                if cs.lc[0].cstat & (1 << sCType.PBIAS) == (1 << sCType.PBIAS):
+                    pbias = self.find_bias(cs, sigsCP, sat)
+
+                if inet > 0 and cs.lc[inet].cstat & (1 << sCType.PBIAS) == \
+                        (1 << sCType.PBIAS):
+                    pbias += self.find_bias(cs, sigsCP, sat, inet)
 
-            # Select reference satellite with highest elevation
+                # note: some services use sign convention different with RTCM
+                if cs.cssrmode in [sc.QZS_CLAS, sc.BDS_PPP, sc.PVS_PPP]:
+                    pbias = -pbias
+                    cbias = -cbias
+
+            # Check for invalid biases
             #
-            if len(idx) > 0:
-                i = idx[np.argmax(el[idx])]
+            if np.isnan(cbias).any() or np.isnan(pbias).any():
+                if self.nav.monlevel > 3:
+                    print("skip invalid cbias/pbias for sat={:d}".format(sat))
+                continue
 
-            # Loop over satellites
+            # Geometric distance corrected for Earth rotation
+            # during flight time
             #
-            for j in idx:
+            r, e[i, :] = geodist(rs[i, :], rr_)
+            _, el[i] = satazel(pos, e[i, :])
+            if el[i] < self.nav.elmin:
+                continue
 
-                # Skip edited observations
-                #
-                if np.any(nav.edt[sat[j]-1, :] > 0):
-                    continue
+            # Shapiro relativistic effect
+            #
+            relatv = shapiro(rs[i, :], rr_)
 
-                # Skip invalid measurements
-                # NOTE: this additional test is included here, since biases or
-                #       antenna offsets may not be available and this zdres()
-                #       returns zero observation residuals!
-                #
-                if y[i, f] == 0.0 or y[j, f] == 0.0:
-                    continue
+            # Tropospheric delay mapping functions
+            #
+            mapfh, mapfw = tropmapf(obs.t, pos, el[i],
+                                    model=self.nav.trpModel)
 
-                # Skip reference satellite i
-                #
-                if i == j:
-                    continue
+            # Tropospheric delay
+            #
+            if self.nav.iono_opt == 2:  # from cssr
+                trop = mapfh*trph*r_hs+mapfw*trpw*r_wet
+            else:
+                trop = mapfh*trop_hs + mapfw*trop_wet
+
+            # Ionospheric delay
+            if self.nav.iono_opt == 2:  # from cssr
+                idx_l = cs.lc[inet].sat_n.index(sat)
+                iono = np.array([40.3e16/(f*f)*stec[idx_l] for f in frq])
+            else:
+                iono = np.zeros(nf)
+
+            # Phase wind-up effect
+            #
+            if self.nav.phw_opt > 0:
+                phw_mode = (False if self.nav.phw_opt == 2 else True)
+                self.nav.phw[sat-1] = windupcorr(obs.t, rs[i, :], vs[i, :],
+                                                 rr_, self.nav.phw[sat-1],
+                                                 full=phw_mode)
+
+                # cycle -> m
+                phw = lam*self.nav.phw[sat-1]
+            else:
+                phw = np.zeros(nf)
+
+            # Select APC reference signals
+            #
+            sig0 = None
+            if cs is not None:
 
-                #  Single-difference measurement
-                #
-                v[nv] = y[i, f] - y[j, f]
+                if cs.cssrmode in (sc.GAL_HAS_SIS, sc.GAL_HAS_IDD,
+                                   sc.QZS_MADOCA):
 
-                # SD line-of-sight vectors
-                #
-                H[nv, 0:3] = -e[i, :] + e[j, :]
+                    if sys == uGNSS.GPS:
+                        sig0 = (rSigRnx("GC1W"), rSigRnx("GC2W"))
+                    elif sys == uGNSS.GAL:
+                        sig0 = (rSigRnx("EC1C"), rSigRnx("EC7Q"))
+                    elif sys == uGNSS.QZS:
+                        sig0 = (rSigRnx("JC1C"), rSigRnx("JC2S"))
+                    elif sys == uGNSS.GLO:
+                        sig0 = (rSigRnx("RC1C"), rSigRnx("RC2C"))
 
-                # SD troposphere
-                #
-                _, mapfwi = gn.tropmapf(obs.t, pos, el[i], model=nav.trpModel)
-                _, mapfwj = gn.tropmapf(obs.t, pos, el[j], model=nav.trpModel)
-
-                idx_i = IT(nav.na)
-                H[nv, idx_i] = mapfwi - mapfwj
-                v[nv] -= (mapfwi - mapfwj)*x[idx_i]
-
-                if nav.monlevel > 2:
-                    nav.fout.write("{}         ztd      ({:3d},{:3d}) {:10.3f} {:10.3f} {:10.3f}\n"
-                                   .format(time2str(obs.t),
-                                           idx_i, idx_i, (mapfwi - mapfwj),
-                                           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 > 2:
-                    nav.fout.write("{} {}-{} ion {} ({:3d},{:3d}) {:10.3f} {:10.3f} {:10.3f} {:10.3f} {:10.3f}\n"
-                                   .format(time2str(obs.t),
-                                           sat2id(sat[i]), sat2id(sat[j]),
-                                           sig, 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
+                elif cs.cssrmode == sc.BDS_PPP:
 
-                    idx_i = IB(sat[i], f, nav.na)
-                    idx_j = IB(sat[j], f, nav.na)
+                    if sys == uGNSS.GPS:
+                        sig0 = (rSigRnx("GC1W"), rSigRnx("GC2W"))
+                    elif sys == uGNSS.BDS:
+                        sig0 = (rSigRnx("CC6I"),)
 
-                    lami = sig.wavelength()
+            # Receiver/satellite antenna offset
+            #
+            antrPR = antModelRx(self.nav, pos, e[i, :], sigsPR, rtype)
+            antrCP = antModelRx(self.nav, pos, e[i, :], sigsCP, rtype)
 
-                    H[nv, idx_i] = +lami
-                    H[nv, idx_j] = -lami
-                    v[nv] -= lami*(x[idx_i] - x[idx_j])
+            if self.nav.ephopt == 4:
 
-                    Ri[nv] = varerr(nav, el[i], f)  # measurement variance
-                    Rj[nv] = varerr(nav, el[j], f)  # measurement variance
+                antsPR = antModelTx(
+                    self.nav, e[i, :], sigsPR, sat, obs.t, rs[i, :])
+                antsCP = antModelTx(
+                    self.nav, e[i, :], sigsCP, sat, obs.t, rs[i, :])
 
-                    nav.vsat[sat[i]-1, f] = 1
-                    nav.vsat[sat[j]-1, f] = 1
+            elif cs is not None and (cs.cssrmode == sc.GAL_HAS_SIS or
+                                     cs.cssrmode == sc.GAL_HAS_IDD or
+                                     cs.cssrmode == sc.QZS_MADOCA or
+                                     cs.cssrmode == sc.BDS_PPP):
 
-                    if nav.monlevel > 2:
-                        nav.fout.write("{} {}-{} amb {} ({:3d},{:3d}) {:10.3f} {:10.3f} {:10.3f} {:10.3f} {:10.3f}\n"
-                                       .format(time2str(obs.t),
-                                               sat2id(sat[i]), sat2id(sat[j]),
-                                               sig, 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])))
+                antsPR = antModelTx(self.nav, e[i, :], sigsPR,
+                                    sat, obs.t, rs[i, :], sig0)
+                antsCP = antModelTx(self.nav, e[i, :], sigsCP,
+                                    sat, obs.t, rs[i, :], sig0)
 
-                else:  # pseudorange
+            else:
 
-                    Ri[nv] = varerr(nav, el[i], f)  # measurement variance
-                    Rj[nv] = varerr(nav, el[j], f)  # measurement variance
+                antsPR = [0.0 for _ in sigsPR]
+                antsCP = [0.0 for _ in sigsCP]
 
-                if nav.monlevel > 1:
-                    nav.fout.write("{} {}-{} res {} ({:3d}) {:10.3f} sig_i {:10.3f} sig_j {:10.3f}\n"
-                                   .format(time2str(obs.t),
-                                           sat2id(sat[i]), sat2id(sat[j]), sig,
-                                           nv, v[nv],
-                                           np.sqrt(Ri[nv]), np.sqrt(Rj[nv])))
+            # Check for invalid values
+            #
+            if antrPR is None or antrCP is None or \
+               antsPR is None or antsCP is None:
+                continue
 
-                nb[b] += 1  # counter for single-differences per signal
-                nv += 1  # counter for single-difference observations
+            # Range correction
+            #
+            prc[i, :] = trop + antrPR + antsPR + iono - cbias
+            cpc[i, :] = trop + antrCP + antsCP - iono - pbias + phw
 
-            b += 1  # counter for signal (pseudrange+carrier-phase)
+            r += relatv - _c*dts[i]
 
-    v = np.resize(v, nv)
-    H = np.resize(H, (nv, nav.nx))
-    R = ddcov(nb, b, Ri, Rj, nv)
+            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 v, H, R
+        return y, e, el
 
+    def sdres(self, obs, x, y, e, sat, el):
+        """
+        SD phase/code residuals
+
+        Parameters
+        ----------
+
+        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
+        """
 
-def kfupdate(x, P, H, v, R):
-    """ Kalman filter measurement update """
-    PHt = P@H.T
-    S = H@PHt+R
-    K = PHt@np.linalg.inv(S)
-    x += K@v
-    P = P - K@H@P
+        nf = self.nav.nf  # number of frequencies (or signals)
+        ns = len(el)  # number of satellites
+        nc = len(obs.sig.keys())  # number of constellations
 
-    return x, P, S
+        mode = 1 if len(y) == ns else 0  # 0:DD,1:SD
 
+        nb = np.zeros(2*nc*nf, dtype=int)
 
-def qcedit(nav, obs, rs, dts, svh):
-    """ Coarse quality control and editing of observations """
+        Ri = np.zeros(ns*nf*2)
+        Rj = np.zeros(ns*nf*2)
 
-    # Predicted position at next epoch
-    #
-    tt = gn.timediff(obs.t, nav.t)
-    rr_ = nav.x[0:3].copy()
-    if nav.pmode > 0:
-        rr_ += nav.x[3:6]*tt
+        nv = 0
+        b = 0
 
-    # Solid Earth tide corrections
-    #
-    # TODO: add solid earth tide displacements
-    #
-    if nav.tidecorr:
-        pos = gn.ecef2pos(rr_)
-        disp = tidedisp(gn.gpst2utc(obs.t), pos)
-    else:
-        disp = np.zeros(3)
-    rr_ += disp
+        H = np.zeros((ns*nf*2, self.nav.nx))
+        v = np.zeros(ns*nf*2)
 
-    # Geodetic position
-    #
-    pos = gn.ecef2pos(rr_)
+        # Geodetic position
+        #
+        pos = ecef2pos(x[0:3])
 
-    # Total number of satellites
-    #
-    ns = gn.uGNSS.MAXSAT
+        # Loop over constellations
+        #
+        for sys in obs.sig.keys():
 
-    # Reset previous editing results
-    #
-    nav.edt = np.zeros((ns, nav.nf), dtype=int)
+            # 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 satellites from one constellation only
+                #
+                idx = self.sysidx(sat, sys)
 
-    # Loop over all satellites
-    #
-    sat = []
-    for i in range(ns):
+                if len(idx) == 0:
+                    continue
 
-        sat_i = i+1
-        sys_i, _ = sat2prn(sat_i)
+                # Select reference satellite with highest elevation
+                #
+                i = idx[np.argmax(el[idx])]
 
-        if sat_i not in obs.sat:
-            nav.edt[i, :] = 1
-            continue
+                # Loop over satellites
+                #
+                for j in idx:
+
+                    # Slant ionospheric delay reference frequency
+                    #
+                    if sys == uGNSS.GLO:
+                        freq0 = obs.sig[sys][uTYP.L][0].frequency(0)
+                    else:
+                        freq0 = obs.sig[sys][uTYP.L][0].frequency()
+
+                    # Select carrier-phase frequency and iono frequency ratio
+                    #
+                    if f < nf:  # carrier
+                        sig = obs.sig[sys][uTYP.L][f]
+                        if sys == uGNSS.GLO:
+                            freq = sig.frequency(self.nav.glo_ch[sat[j]])
+                        else:
+                            freq = sig.frequency()
+                        mu = -(freq0/freq)**2
+                    else:  # code
+                        sig = obs.sig[sys][uTYP.C][f % nf]
+                        if sys == uGNSS.GLO:
+                            freq = sig.frequency(self.nav.glo_ch[sat[j]])
+                        else:
+                            freq = sig.frequency()
+                        mu = +(freq0/freq)**2
+
+                    # Skip edited observations
+                    #
+                    if np.any(self.nav.edt[sat[j]-1, :] > 0):
+                        continue
+
+                    # Skip invalid measurements
+                    # NOTE: this additional test is included here,
+                    #       since biases or antenna offsets may not be
+                    #       available and this zdres()
+                    #       returns zero observation residuals!
+                    #
+                    if y[i, f] == 0.0 or y[j, f] == 0.0:
+                        continue
+
+                    # Skip reference satellite i
+                    #
+                    if i == j:
+                        continue
+
+                    if mode == 0:  # DD
+                        v[nv] = (y[i, f]-y[i+ns, f])-(y[j, f]-y[j+ns, f])
+                    else:
+                        #  Single-difference measurement
+                        #
+                        v[nv] = y[i, f] - y[j, f]
+
+                    # SD line-of-sight vectors
+                    #
+                    H[nv, 0:3] = -e[i, :] + e[j, :]
+
+                    if self.nav.ntrop > 0:  # tropo is estimated
+
+                        # SD troposphere
+                        #
+                        _, mapfwi = tropmapf(
+                            obs.t, pos, el[i], model=self.nav.trpModel)
+                        _, mapfwj = tropmapf(
+                            obs.t, pos, el[j], model=self.nav.trpModel)
+
+                        idx_i = self.IT(self.nav.na)
+                        H[nv, idx_i] = mapfwi - mapfwj
+                        v[nv] -= (mapfwi - mapfwj)*x[idx_i]
+
+                        if self.nav.monlevel > 2:
+                            self.nav.fout.write(
+                                fmt_ztd
+                                .format(time2str(obs.t), idx_i, idx_i,
+                                        (mapfwi - mapfwj),
+                                        x[self.IT(self.nav.na)],
+                                        np.sqrt(self.nav.P[
+                                            self.IT(self.nav.na),
+                                            self.IT(self.nav.na)])))
+
+                    if self.nav.niono > 0:  # iono is estimated
+                        # SD ionosphere
+                        #
+                        idx_i = self.II(sat[i], self.nav.na)
+                        idx_j = self.II(sat[j], self.nav.na)
+                        H[nv, idx_i] = +mu
+                        H[nv, idx_j] = -mu
+                        v[nv] -= mu*(x[idx_i] - x[idx_j])
+
+                        if self.nav.monlevel > 2:
+                            self.nav.fout.write(
+                                fmt_ion
+                                .format(time2str(obs.t),
+                                        sat2id(sat[i]), sat2id(sat[j]),
+                                        sig, idx_i, idx_j, mu,
+                                        x[idx_i], x[idx_j],
+                                        np.sqrt(self.nav.P[idx_i, idx_i]),
+                                        np.sqrt(self.nav.P[idx_j, idx_j])))
+
+                    # SD ambiguity
+                    #
+                    if f < nf:  # carrier-phase
+
+                        idx_i = self.IB(sat[i], f, self.nav.na)
+                        idx_j = self.IB(sat[j], f, self.nav.na)
+
+                        if sys == uGNSS.GLO:
+                            lami = sig.wavelength(self.nav.glo_ch[sat[i]])
+                            lamj = sig.wavelength(self.nav.glo_ch[sat[j]])
+                        else:
+                            lami = sig.wavelength()
+                            lamj = lami
+
+                        H[nv, idx_i] = +lami
+                        H[nv, idx_j] = -lamj
+                        v[nv] -= lami*(x[idx_i] - x[idx_j])
+
+                        # measurement variance
+                        Ri[nv] = self.varerr(self.nav, el[i], f)
+                        # measurement variance
+                        Rj[nv] = self.varerr(self.nav, el[j], f)
+
+                        self.nav.vsat[sat[i]-1, f] = 1
+                        self.nav.vsat[sat[j]-1, f] = 1
+
+                        if self.nav.monlevel > 2:
+                            self.nav.fout.write(
+                                fmt_amb
+                                .format(time2str(obs.t),
+                                        sat2id(sat[i]), sat2id(sat[j]),
+                                        sig, idx_i, idx_j, lami, lamj,
+                                        x[idx_i], x[idx_j],
+                                        np.sqrt(self.nav.P[idx_i, idx_i]),
+                                        np.sqrt(self.nav.P[idx_j, idx_j])))
+
+                    else:  # pseudorange
+
+                        # measurement variance
+                        Ri[nv] = self.varerr(self.nav, el[i], f)
+                        # measurement variance
+                        Rj[nv] = self.varerr(self.nav, el[j], f)
+
+                    if self.nav.monlevel > 1:
+                        self.nav.fout.write(
+                            fmt_res
+                            .format(time2str(obs.t),
+                                    sat2id(sat[i]), sat2id(sat[j]), sig,
+                                    nv, v[nv],
+                                    np.sqrt(Ri[nv]), np.sqrt(Rj[nv])))
+
+                    nb[b] += 1  # counter for single-differences per signal
+                    nv += 1  # counter for single-difference observations
+
+                b += 1  # counter for signal (pseudrange+carrier-phase)
+
+        v = np.resize(v, nv)
+        H = np.resize(H, (nv, self.nav.nx))
+        R = self.ddcov(nb, b, Ri, Rj, nv)
+
+        return v, H, R
+
+    def ddcov(self, 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
+
+    def kfupdate(self, x, P, H, v, R):
+        """
+        Kalman filter measurement update.
+
+        Parameters:
+        x (ndarray): State estimate vector
+        P (ndarray): State covariance matrix
+        H (ndarray): Observation model matrix
+        v (ndarray): Innovation vector
+                     (residual between measurement and prediction)
+        R (ndarray): Measurement noise covariance
+
+        Returns:
+        x (ndarray): Updated state estimate vector
+        P (ndarray): Updated state covariance matrix
+        S (ndarray): Innovation covariance matrix
+        """
 
-        # Check satellite exclusion
+        PHt = P@H.T
+        S = H@PHt+R
+        K = PHt@np.linalg.inv(S)
+        x += K@v
+        # P = P - K@H@P
+        IKH = np.eye(P.shape[0])-K@H
+        P = IKH@P@IKH.T + K@R@K.T  # Joseph stabilized version
+
+        return x, P, S
+
+    def restamb(self, bias, nb):
+        """ restore SD ambiguity """
+        nv = 0
+        xa = self.nav.x.copy()
+        xa[0:self.nav.na] = self.nav.xa[0:self.nav.na]
+
+        for m in range(uGNSS.GNSSMAX):
+            for f in range(self.nav.nf):
+                n = 0
+                index = []
+                for i in range(uGNSS.MAXSAT):
+                    sys, _ = sat2prn(i+1)
+                    if sys != m or self.nav.fix[i, f] != 2:
+                        continue
+                    index.append(self.IB(i+1, f, self.nav.na))
+                    n += 1
+                if n < 2:
+                    continue
+                xa[index[0]] = self.nav.x[index[0]]
+                for i in range(1, n):
+                    xa[index[i]] = xa[index[0]]-bias[nv]
+                    nv += 1
+        return xa
+
+    def ddidx(self, nav, sat):
+        """ index for SD to DD transformation matrix D """
+        nb = 0
+        n = uGNSS.MAXSAT
+        na = nav.na
+        ix = np.zeros((n, 2), dtype=int)
+        nav.fix = np.zeros((n, nav.nf), dtype=int)
+        for m in range(uGNSS.GNSSMAX):
+            k = na
+            for f in range(nav.nf):
+                for i in range(k, k+n):
+                    sat_i = i-k+1
+                    sys, _ = sat2prn(sat_i)
+                    if (sys != m):
+                        continue
+                    if sat_i not in sat or nav.x[i] == 0.0 \
+                       or nav.vsat[sat_i-1, f] == 0:
+                        continue
+                    if nav.el[sat_i-1] >= nav.elmaskar:
+                        nav.fix[sat_i-1, f] = 2
+                        break
+                    else:
+                        nav.fix[sat_i-1, f] = 1
+                for j in range(k, k+n):
+                    sat_j = j-k+1
+                    sys, _ = sat2prn(sat_j)
+                    if (sys != m):
+                        continue
+                    if i == j or sat_j not in sat or nav.x[j] == 0.0 \
+                       or nav.vsat[sat_j-1, f] == 0:
+                        continue
+                    if nav.el[sat_j-1] >= nav.elmaskar:
+                        ix[nb, :] = [i, j]
+                        nb += 1
+                        nav.fix[sat_j-1, f] = 2
+                k += n
+        ix = np.resize(ix, (nb, 2))
+        return ix
+
+    def resamb_lambda(self, sat):
+        """ resolve integer ambiguity using LAMBDA method """
+        nx = self.nav.nx
+        na = self.nav.na
+        xa = np.zeros(na)
+        ix = self.ddidx(self.nav, sat)
+        nb = len(ix)
+        if nb <= 0:
+            print("no valid DD")
+            return -1, -1
+
+        # y=D*xc, Qb=D*Qc*D', Qab=Qac*D'
+        y = self.nav.x[ix[:, 0]]-self.nav.x[ix[:, 1]]
+        DP = self.nav.P[ix[:, 0], na:nx]-self.nav.P[ix[:, 1], na:nx]
+        Qb = DP[:, ix[:, 0]-na]-DP[:, ix[:, 1]-na]
+        Qab = self.nav.P[0:na, ix[:, 0]]-self.nav.P[0:na, ix[:, 1]]
+
+        # MLAMBDA ILS
+        b, s = mlambda(y, Qb)
+        if s[0] <= 0.0 or s[1]/s[0] >= self.nav.thresar:
+            self.nav.xa = self.nav.x[0:na].copy()
+            self.nav.Pa = self.nav.P[0:na, 0:na].copy()
+            bias = b[:, 0]
+            y -= b[:, 0]
+            K = Qab@np.linalg.inv(Qb)
+            self.nav.xa -= K@y
+            self.nav.Pa -= K@Qab.T
+
+            # restore SD ambiguity
+            xa = self.restamb(bias, nb)
+        else:
+            nb = 0
+
+        return nb, xa
+
+    def holdamb(self, xa):
+        """ hold integer ambiguity """
+        nb = self.nav.nx-self.nav.na
+        v = np.zeros(nb)
+        H = np.zeros((nb, self.nav.nx))
+        nv = 0
+        for m in range(uGNSS.GNSSMAX):
+            for f in range(self.nav.nf):
+                n = 0
+                index = []
+                for i in range(uGNSS.MAXSAT):
+                    sys, _ = sat2prn(i+1)
+                    if sys != m or self.nav.fix[i, f] != 2:
+                        continue
+                    index.append(self.IB(i+1, f, self.nav.na))
+                    n += 1
+                    self.nav.fix[i, f] = 3  # hold
+                # constraint to fixed ambiguity
+                for i in range(1, n):
+                    v[nv] = (xa[index[0]]-xa[index[i]]) - \
+                        (self.nav.x[index[0]]-self.nav.x[index[i]])
+                    H[nv, index[0]] = 1.0
+                    H[nv, index[i]] = -1.0
+                    nv += 1
+        if nv > 0:
+            R = np.eye(nv)*self.VAR_HOLDAMB
+            # update states with constraints
+            self.nav.x, self.nav.P, _ = self.kfupdate(
+                self.nav.x, self.nav.P, H[0:nv, :], v[0:nv], R)
+        return 0
+
+    def qcedit(self, obs, rs, dts, svh, rr=None):
+        """ Coarse quality control and editing of observations """
+
+        # Predicted position at next epoch
         #
-        if sat_i in nav.excl_sat:
-            nav.edt[i, :] = 1
-            if nav.monlevel > 0:
-                nav.fout.write("{}  {} - edit - satellite excluded\n"
-                               .format(time2str(obs.t), sat2id(sat_i)))
-            continue
-
-        j = np.where(obs.sat == sat_i)[0][0]
+        tt = timediff(obs.t, self.nav.t)
+        if rr is None:
+            rr_ = self.nav.x[0:3].copy()
+            if self.nav.pmode > 0:
+                rr_ += self.nav.x[3:6]*tt
+        else:
+            rr_ = rr
 
-        # Check for valid orbit and clock offset
+        # Solid Earth tide corrections
         #
-        if np.isnan(rs[j, :]).any() or np.isnan(dts[j]):
-            nav.edt[i, :] = 1
-            if nav.monlevel > 0:
-                nav.fout.write("{}  {} - edit - invalid eph\n"
-                               .format(time2str(obs.t), sat2id(sat_i)))
-            continue
-
-        # Check satellite health
+        if self.nav.tidecorr == uTideModel.SIMPLE:
+            pos = ecef2pos(rr_)
+            disp = tidedisp(gpst2utc(obs.t), pos)
+        elif self.nav.tidecorr == uTideModel.IERS2010:
+            pos = ecef2pos(rr_)
+            disp = tidedispIERS2010(gpst2utc(obs.t), pos)
+        else:
+            disp = np.zeros(3)
+        rr_ += disp
+
+        # Geodetic position
         #
-        if svh[j] > 0:
-            nav.edt[i, :] = 1
-            if nav.monlevel > 0:
-                nav.fout.write("{}  {} - edit - satellite unhealthy\n"
-                               .format(time2str(obs.t), sat2id(sat_i)))
-            continue
-
-        # Check elevation angle
+        pos = ecef2pos(rr_)
+
+        # Total number of satellites
         #
-        _, e = gn.geodist(rs[j, :], rr_)
-        _, el = gn.satazel(pos, e)
-        if el < nav.elmin:
-            nav.edt[i][:] = 1
-            if nav.monlevel > 0:
-                nav.fout.write("{}  {} - edit - low elevation {:5.1f} deg\n"
-                               .format(time2str(obs.t), sat2id(sat_i),
-                                       np.rad2deg(el)))
-            continue
-
-        # Pseudorange, carrier-phase and C/N0 signals
+        ns = uGNSS.MAXSAT
+
+        # Reset previous editing results
         #
-        sigsPR = obs.sig[sys_i][gn.uTYP.C]
-        sigsCP = obs.sig[sys_i][gn.uTYP.L]
-        sigsCN = obs.sig[sys_i][gn.uTYP.S]
+        self.nav.edt = np.zeros((ns, self.nav.nf), dtype=int)
 
-        # Loop over signals
+        # Loop over all satellites
         #
-        for f in range(nav.nf):
+        sat = []
+        for i in range(ns):
+
+            sat_i = i+1
+            sys_i, _ = sat2prn(sat_i)
+
+            if sat_i not in obs.sat:
+                self.nav.edt[i, :] = 1
+                continue
 
-            # Cycle  slip check by LLI
+            # Check satellite exclusion
             #
-            if obs.lli[j, f] == 1:
-                nav.edt[i, f] = 1
-                if nav.monlevel > 0:
-                    nav.fout.write("{}  {} - edit {:4s} - LLI\n"
-                                   .format(time2str(obs.t), sat2id(sat_i),
-                                           sigsCP[f].str()))
+            if sat_i in self.nav.excl_sat:
+                self.nav.edt[i, :] = 1
+                if self.nav.monlevel > 0:
+                    self.nav.fout.write("{}  {} - edit - satellite excluded\n"
+                                        .format(time2str(obs.t),
+                                                sat2id(sat_i)))
                 continue
 
-            # Check for measurement consistency
+            j = np.where(obs.sat == sat_i)[0][0]
+
+            # Check for valid orbit and clock offset
             #
-            if obs.P[j, f] == 0.0:
-                nav.edt[i, f] = 1
-                if nav.monlevel > 0:
-                    nav.fout.write("{}  {} - edit {:4s} - invalid PR obs\n"
-                                   .format(time2str(obs.t), sat2id(sat_i),
-                                           sigsPR[f].str()))
+            if np.isnan(rs[j, :]).any() or np.isnan(dts[j]):
+                self.nav.edt[i, :] = 1
+                if self.nav.monlevel > 0:
+                    self.nav.fout.write("{}  {} - edit - invalid eph\n"
+                                        .format(time2str(obs.t),
+                                                sat2id(sat_i)))
                 continue
 
-            if obs.L[j, f] == 0.0:
-                nav.edt[i, f] = 1
-                if nav.monlevel > 0:
-                    nav.fout.write("{}  {} - edit {:4s} - invalid CP obs\n"
-                                   .format(time2str(obs.t), sat2id(sat_i),
-                                           sigsCP[f].str()))
+            # Check satellite health
+            #
+            if svh[j] > 0:
+                self.nav.edt[i, :] = 1
+                if self.nav.monlevel > 0:
+                    self.nav.fout.write("{}  {} - edit - satellite unhealthy\n"
+                                        .format(time2str(obs.t),
+                                                sat2id(sat_i)))
                 continue
 
-            # Check C/N0
+            # Check elevation angle
             #
-            cnr_min = nav.cnr_min_gpy if sigsCN[f].isGPS_PY() else nav.cnr_min
-            if obs.S[j, f] < cnr_min:
-                nav.edt[i, f] = 1
-                if nav.monlevel > 0:
-                    nav.fout.write("{}  {} - edit {:4s} - low C/N0 {:4.1f} dB-Hz\n"
-                                   .format(time2str(obs.t), sat2id(sat_i),
-                                           sigsCN[f].str(), obs.S[j, f]))
+            _, e = geodist(rs[j, :], rr_)
+            _, el = satazel(pos, e)
+            if el < self.nav.elmin:
+                self.nav.edt[i][:] = 1
+                if self.nav.monlevel > 0:
+                    self.nav.fout.write(
+                        "{}  {} - edit - low elevation {:5.1f} deg\n"
+                        .format(time2str(obs.t), sat2id(sat_i),
+                                np.rad2deg(el)))
                 continue
 
-        # Store satellite which have passed all tests
-        #
-        if np.any(nav.edt[i, :] > 0):
-            continue
-
-        sat.append(sat_i)
-
-    return np.array(sat, dtype=int)
-
-
-def ppppos(nav, obs, cs=None, orb=None, bsx=None):
-    """
-    PPP positioning
-    """
-
-    # Skip empty epochs
-    #
-    if len(obs.sat) == 0:
-        return
-
-    # GNSS satellite positions, velocities and clock offsets for all satellites
-    # in RINEX observations
-    #
-    rs, vs, dts, svh, nsat = satposs(obs, nav, cs=cs, orb=orb)
-
-    if nsat < 6:
-        print(" too few satellites < 6: nsat={:d}".format(nsat))
-        return
-
-    # Editing of observations
-    #
-    sat_ed = qcedit(nav, obs, rs, dts, svh)
-
-    # Kalman filter time propagation, initialization of ambiguities and iono
-    #
-    udstate(nav, obs)
-
-    xa = np.zeros(nav.nx)
-    xp = nav.x.copy()
-
-    # Non-differential residuals
-    #
-    yu, eu, elu = zdres(nav, obs, cs, bsx, rs, vs, dts, xp[0:3])
-
-    # Select satellites having passed quality control
-    #
-    # index of valid sats in obs.sat
-    iu = np.where(np.isin(obs.sat, sat_ed))[0]
-    sat = obs.sat[iu]
-    y = yu[iu, :]
-    e = eu[iu, :]
-    el = elu[iu]
-
-    # Store reduced satellite list
-    # NOTE: where are working on a reduced list of observations from here on
-    #
-    nav.sat = sat
-    nav.el[sat-1] = el  # needed in rtk.ddidx()
-    nav.y = y
-    ns = len(sat)
-
-    # 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
-
-    # SD residuals
-    #
-    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 residuals after measurement update
-    #
-    yu, eu, elu = zdres(nav, obs, cs, bsx, rs, vs, dts, xp[0:3])
-    y = yu[iu, :]
-    e = eu[iu, :]
-    ny = y.shape[0]
-    if ny < 6:
-        return -1
-
-    # 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
-        nav.ns = 0
-        for i in range(ns):
-            j = sat[i]-1
-            for f in range(nav.nf):
-                if nav.vsat[j, f] == 0:
+            # Pseudorange, carrier-phase and C/N0 signals
+            #
+            sigsPR = obs.sig[sys_i][uTYP.C]
+            sigsCP = obs.sig[sys_i][uTYP.L]
+            sigsCN = obs.sig[sys_i][uTYP.S]
+
+            # Loop over signals
+            #
+            for f in range(self.nav.nf):
+
+                # Cycle  slip check by LLI
+                #
+                if obs.lli[j, f] == 1:
+                    self.nav.edt[i, f] = 1
+                    if self.nav.monlevel > 0:
+                        self.nav.fout.write("{}  {} - edit {:4s} - LLI\n"
+                                            .format(time2str(obs.t),
+                                                    sat2id(sat_i),
+                                                    sigsCP[f].str()))
+                    continue
+
+                # Check for measurement consistency
+                #
+                if obs.P[j, f] == 0.0:
+                    self.nav.edt[i, f] = 1
+                    if self.nav.monlevel > 0:
+                        self.nav.fout.write(
+                            "{}  {} - edit {:4s} - invalid PR obs\n"
+                            .format(time2str(obs.t),
+                                    sat2id(sat_i),
+                                    sigsPR[f].str()))
+                    continue
+
+                if obs.L[j, f] == 0.0:
+                    self.nav.edt[i, f] = 1
+                    if self.nav.monlevel > 0:
+                        self.nav.fout.write(
+                            "{}  {} - edit {:4s} - invalid CP obs\n"
+                            .format(time2str(obs.t),
+                                    sat2id(sat_i),
+                                    sigsCP[f].str()))
+                    continue
+
+                # Check C/N0
+                #
+                cnr_min = self.nav.cnr_min_gpy \
+                    if sigsCN[f].isGPS_PY() else self.nav.cnr_min
+                if obs.S[j, f] < cnr_min:
+                    self.nav.edt[i, f] = 1
+                    if self.nav.monlevel > 0:
+                        self.nav.fout.write(
+                            "{}  {} - edit {:4s} - low C/N0 {:4.1f} dB-Hz\n"
+                            .format(time2str(obs.t),
+                                    sat2id(sat_i),
+                                    sigsCN[f].str(),
+                                    obs.S[j, f]))
                     continue
-                nav.outc[j, f] = 0
-                if f == 0:
-                    nav.ns += 1
-    else:
-        nav.smode = 0
-
-    nav.smode = 5  # 4: fixed ambiguities, 5: float ambiguities
-
-    if nav.armode > 0:
-        nb, xa = resamb_lambda(nav, sat)
-        if nb > 0:
-            # Use position with fixed ambiguities xa
-            yu, eu, elu = zdres(nav, obs, cs, bsx, rs, vs, dts, xa[0:3])
-            y = yu[iu, :]
-            e = eu[iu, :]
-            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
+
+            # Store satellite which have passed all tests
+            #
+            if np.any(self.nav.edt[i, :] > 0):
+                continue
+
+            sat.append(sat_i)
+
+        return np.array(sat, dtype=int)
+
+    def base_process(self, obs, obsb, rs, dts, svh):
+        """ processing for base station in RTK
+            (implemented in rtkpos) """
+        return None, None, None, None
+
+    def process(self, obs, cs=None, orb=None, bsx=None, obsb=None):
+        """
+        PPP/PPP-RTK/RTK positioning
+        """
+
+        # Skip empty epochs
+        #
+        if len(obs.sat) == 0:
+            return
+
+        # GNSS satellite positions, velocities and clock offsets
+        # for all satellite in RINEX observations
+        #
+        rs, vs, dts, svh, nsat = satposs(obs, self.nav, cs=cs, orb=orb)
+
+        if nsat < 6:
+            print(" too few satellites < 6: nsat={:d}".format(nsat))
+            return
+
+        # Editing of observations
+        #
+        sat_ed = self.qcedit(obs, rs, dts, svh)
+
+        if obsb is None:  # PPP/PPP-RTK
+            # Select satellites having passed quality control
+            #
+            # index of valid sats in obs.sat
+            iu = np.where(np.isin(obs.sat, sat_ed))[0]
+            ns = len(iu)
+            y = np.zeros((ns, self.nav.nf*2))
+            e = np.zeros((ns, 3))
+
+            obs_ = obs
+        else:  # RTK
+            y, e, iu, obs_ = self.base_process(obs, obsb, rs, dts, svh)
+            ns = len(iu)
+
+        if ns < 6:
+            print(" too few satellites < 6: ns={:d}".format(ns))
+            return
+
+        # Kalman filter time propagation, initialization of ambiguities
+        # and iono
+        #
+        self.udstate(obs_)
+
+        xa = np.zeros(self.nav.nx)
+        xp = self.nav.x.copy()
+
+        # Non-differential residuals
+        #
+        yu, eu, elu = self.zdres(obs, cs, bsx, rs, vs, dts, xp[0:3])
+
+        # Select satellites having passed quality control
+        #
+        # index of valid sats in obs.sat
+        sat = obs.sat[iu]
+        y[:ns, :] = yu[iu, :]
+        e[:ns, :] = eu[iu, :]
+        el = elu[iu]
+
+        # Store reduced satellite list
+        # NOTE: where are working on a reduced list of observations
+        # from here on
+        #
+        self.nav.sat = sat
+        self.nav.el[sat-1] = el  # needed in rtk.ddidx()
+        self.nav.y = y
+        ns = len(sat)
+
+        # Check if observations of at least 6 satellites are left over
+        # after editing
+        #
+        ny = y.shape[0]
+        if ny < 6:
+            self.nav.P[np.diag_indices(3)] = 1.0
+            self.nav.smode = 5
+            return -1
+
+        # SD residuals
+        #
+        v, H, R = self.sdres(obs, xp, y, e, sat, el)
+        Pp = self.nav.P.copy()
+
+        # Kalman filter measurement update
+        #
+        xp, Pp, _ = self.kfupdate(xp, Pp, H, v, R)
+
+        # Non-differential residuals after measurement update
+        #
+        yu, eu, elu = self.zdres(obs, cs, bsx, rs, vs, dts, xp[0:3])
+        y = yu[iu, :]
+        e = eu[iu, :]
+        ny = y.shape[0]
+        if ny < 6:
+            return -1
+
+        # Residuals for float solution
+        #
+        v, H, R = self.sdres(obs, xp, y, e, sat, el)
+        if self.valpos(v, R):
+            self.nav.x = xp
+            self.nav.P = Pp
+            self.nav.ns = 0
+            for i in range(ns):
+                j = sat[i]-1
+                for f in range(self.nav.nf):
+                    if self.nav.vsat[j, f] == 0:
+                        continue
+                    self.nav.outc[j, f] = 0
+                    if f == 0:
+                        self.nav.ns += 1
+        else:
+            self.nav.smode = 0
+
+        self.nav.smode = 5  # 4: fixed ambiguities, 5: float ambiguities
+
+        if self.nav.armode > 0:
+            nb, xa = self.resamb_lambda(sat)
+            if nb > 0:
+                # Use position with fixed ambiguities xa
+                yu, eu, elu = self.zdres(obs, cs, bsx, rs, vs, dts, xa[0:3])
+                y = yu[iu, :]
+                e = eu[iu, :]
+                v, H, R = self.sdres(obs, xa, y, e, sat, el)
+                # R <= Q=H'PH+R  chisq<max_inno[3] (0.5)
+                if self.valpos(v, R):
+                    if self.nav.armode == 3:     # fix and hold
+                        self.holdamb(xa)    # hold fixed ambiguity
+                    self.nav.smode = 4           # fix
+
+        # Store epoch for solution
+        #
+        self.nav.t = obs.t
+
+        return 0
diff --git a/src/cssrlib/rinex.py b/src/cssrlib/rinex.py
index c91167d..31f631d 100644
--- a/src/cssrlib/rinex.py
+++ b/src/cssrlib/rinex.py
@@ -7,7 +7,8 @@
 from cssrlib.gnss import bdt2gpst, time2bdt
 from cssrlib.gnss import gpst2time, bdt2time, epoch2time, timediff, gtime_t
 from cssrlib.gnss import prn2sat, char2sys, timeget, utc2gpst, time2epoch
-from cssrlib.gnss import Eph, Obs, sat2id, sat2prn, gpst2bdt, time2gpst
+from cssrlib.gnss import Eph, Geph, Obs, sat2id, sat2prn, gpst2bdt, time2gpst
+from cssrlib.gnss import timeadd, id2sat
 
 
 class pclk_t:
@@ -43,6 +44,7 @@ def __init__(self):
         # 0:LNAV,INAV,D1/D2, 1:CNAV/CNAV1/FNAV, 2: CNAV2, 3: CNAV3,
         # 4:FDMA, 5:SBAS
         self.mode_nav = 0
+        self.glo_ch = {}
 
     def setSignals(self, sigList):
         """ define the signal list for each constellation """
@@ -114,6 +116,14 @@ def flt(self, u, c=-1):
             return 0.0
         return float(u.replace("D", "E"))
 
+    def adjday(self, t: gtime_t, t0: gtime_t):
+        tt = timediff(t, t0)
+        if tt < -43200.0:
+            return timeadd(t, 86400.0)
+        if tt > 43200.0:
+            return timeadd(t, -86400.0)
+        return t
+
     def decode_time(self, s, ofst=0, slen=2):
         year = int(s[ofst+0:ofst+4])
         month = int(s[ofst+5:ofst+7])
@@ -124,7 +134,7 @@ def decode_time(self, s, ofst=0, slen=2):
         t = epoch2time([year, month, day, hour, minute, sec])
         return t
 
-    def decode_nav(self, navfile, nav):
+    def decode_nav(self, navfile, nav, append=False):
         """
         Decode RINEX Navigation message from file
 
@@ -132,7 +142,10 @@ def decode_nav(self, navfile, nav):
 
         """
 
-        nav.eph = []
+        if not append:
+            nav.eph = []
+            nav.geph = []
+
         with open(navfile, 'rt') as fnav:
             for line in fnav:
                 if line[60:73] == 'END OF HEADER':
@@ -235,18 +248,86 @@ def decode_nav(self, navfile, nav):
                         elif m == 'CNV3':
                             self.mode_nav = 3
                         elif m == 'FDMA':
-                            self.mode_nav = 4
+                            self.mode_nav = 0
                         elif m == 'SBAS':
-                            self.mode_nav = 5
+                            self.mode_nav = 0
                         line = fnav.readline()
 
+                elif self.ver >= 3.0:  # RINEX 3.0.x
+                    self.mode_nav = 0
+
                 # Process ephemeris information
                 #
                 sys = char2sys(line[0])
 
                 # Skip undesired constellations
                 #
-                if sys not in (uGNSS.GPS, uGNSS.GAL, uGNSS.QZS, uGNSS.BDS):
+                if sys == uGNSS.GLO:
+                    prn = int(line[1:3])
+                    sat = prn2sat(sys, prn)
+                    geph = Geph(sat)
+                    geph.pos = np.zeros(3)
+                    geph.vel = np.zeros(3)
+                    geph.acc = np.zeros(3)
+
+                    geph.mode = self.mode_nav
+                    toc = self.decode_time(line, 4)
+                    week, tocs = time2gpst(toc)
+                    toc = gpst2time(week,
+                                    np.floor((tocs+450.0)/900.0)*900.0)
+                    dow = int(tocs//86400.0)
+
+                    geph.taun = -self.flt(line, 1)
+                    geph.gamn = self.flt(line, 2)
+                    t0 = self.flt(line, 3)
+
+                    tod = t0 % 86400.0
+                    tof = gpst2time(week, tod + dow*86400.0)
+                    tof = self.adjday(tof, toc)
+
+                    geph.toe = utc2gpst(toc)
+                    geph.tof = utc2gpst(tof)
+
+                    # iode = Tb(7bit)
+                    geph.iode = int(((tocs+10800.0) % 86400)/900.0+0.5)
+
+                    line = fnav.readline()  # line #1
+                    geph.pos[0] = self.flt(line, 0)*1e3
+                    geph.vel[0] = self.flt(line, 1)*1e3
+                    geph.acc[0] = self.flt(line, 2)*1e3
+                    geph.svh = int(self.flt(line, 3))
+
+                    line = fnav.readline()  # line #2
+                    geph.pos[1] = self.flt(line, 0)*1e3
+                    geph.vel[1] = self.flt(line, 1)*1e3
+                    geph.acc[1] = self.flt(line, 2)*1e3
+                    geph.frq = int(self.flt(line, 3))
+
+                    if geph.frq > 128:
+                        geph.frq -= 256
+
+                    line = fnav.readline()  # line #3
+                    geph.pos[2] = self.flt(line, 0)*1e3
+                    geph.vel[2] = self.flt(line, 1)*1e3
+                    geph.acc[2] = self.flt(line, 2)*1e3
+                    geph.age = int(self.flt(line, 3))
+
+                    # Use GLONASS line #4 only from RINEX v3.05 onwards
+                    #
+                    if self.ver >= 3.05:
+
+                        line = fnav.readline()  # line #4
+
+                        # b7-8: M, b6: P4, b5: P3, b4: P2, b2-3: P1, b0-1: P
+                        geph.status = int(self.flt(line, 0))
+                        geph.dtaun = -self.flt(line, 1)
+                        geph.urai = int(self.flt(line, 2))
+                        # svh = int(self.flt(line, 3))
+
+                    nav.geph.append(geph)
+                    continue
+
+                elif sys not in (uGNSS.GPS, uGNSS.GAL, uGNSS.QZS, uGNSS.BDS):
                     continue
 
                 prn = int(line[1:3])
@@ -255,6 +336,10 @@ def decode_nav(self, navfile, nav):
                 sat = prn2sat(sys, prn)
                 eph = Eph(sat)
 
+                eph.urai = np.zeros(3, dtype=int)
+                eph.sisai = np.zeros(4, dtype=int)
+                eph.isc = np.zeros(6)
+
                 eph.mode = self.mode_nav
 
                 eph.toc = self.decode_time(line, 4)
@@ -264,10 +349,15 @@ def decode_nav(self, navfile, nav):
 
                 line = fnav.readline()  # line #1
 
-                if self.mode_nav > 0:
-                    eph.Adot = self.flt(line, 0)
-                else:
+                if sys == uGNSS.GAL:
                     eph.iode = int(self.flt(line, 0))
+                    eph.iodc = eph.iode
+                else:
+                    if self.mode_nav > 0:
+                        eph.Adot = self.flt(line, 0)
+                    else:
+                        eph.iode = int(self.flt(line, 0))
+
                 eph.crs = self.flt(line, 1)
                 eph.deln = self.flt(line, 2)
                 eph.M0 = self.flt(line, 3)
@@ -293,7 +383,15 @@ def decode_nav(self, navfile, nav):
 
                 line = fnav.readline()  # line #5
                 eph.idot = self.flt(line, 0)
-                if self.mode_nav > 0:
+
+                if sys == uGNSS.GAL or self.mode_nav == 0:
+                    eph.code = int(self.flt(line, 1))  # source for GAL
+                    eph.week = int(self.flt(line, 2))
+
+                    if sys == uGNSS.GAL and self.ver < 4.0:
+                        eph.mode = 1 if eph.code & 0x2 else 0
+
+                else:
                     eph.delnd = self.flt(line, 1)
                     if sys == uGNSS.BDS:
                         eph.sattype = int(self.flt(line, 2))
@@ -301,9 +399,6 @@ def decode_nav(self, navfile, nav):
                     else:
                         eph.urai[0] = int(self.flt(line, 2))
                         eph.urai[1] = int(self.flt(line, 3))
-                else:
-                    eph.code = int(self.flt(line, 1))  # source for GAL
-                    eph.week = int(self.flt(line, 2))
 
                 line = fnav.readline()  # line #6
                 if sys == uGNSS.BDS and self.mode_nav > 0:
@@ -322,13 +417,18 @@ def decode_nav(self, navfile, nav):
                             eph.urai[2] = int(self.flt(line, 3))
                     elif sys == uGNSS.GAL:
                         tgd_b = float(self.flt(line, 3))
-                        if (eph.code >> 9) & 1:
+                        if (eph.code >> 9) & 1:  # E5b,E1
+                            eph.tgd_b = eph.tgd
                             eph.tgd = tgd_b
-                        else:
-                            eph.iodc = int(self.flt(line, 3))
+                        else:  # E5a,E1
+                            eph.tgd_b = tgd_b
                     elif sys == uGNSS.BDS:
                         eph.tgd_b = float(self.flt(line, 3))  # tgd2 B2/B3
 
+                    if sys == uGNSS.QZS:
+                        eph.code = eph.svh & 0x11  # L1C/A:0x01 or L1C/B:0x10
+                        eph.svh = eph.svh & 0xEE   # mask L1C/A, L1C/B health
+
                 if self.mode_nav < 3:
                     line = fnav.readline()  # line #7
                     if sys == uGNSS.BDS:
@@ -344,6 +444,9 @@ def decode_nav(self, navfile, nav):
                             eph.tgd = float(self.flt(line, 2))    # tgd_B1Cp
                             eph.tgd_b = float(self.flt(line, 3))  # tgd_B2ap
 
+                    elif sys == uGNSS.GAL:
+                        tot = int(self.flt(line, 0))
+
                     else:
                         if self.mode_nav > 0 and sys != uGNSS.GAL:
                             eph.isc[0] = self.flt(line, 0)
@@ -488,6 +591,15 @@ 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()])
+            elif 'GLONASS SLOT / FRQ #' in line[60:]:
+                nsat = int(line[0:3])
+                for i in range(nsat):
+                    if i > 0 and i % 8 == 0:
+                        line = self.fobs.readline()
+                    j = i % 8
+                    sat = id2sat(line[4+7*j:7+7*j])
+                    ch = int(line[8+7*j: 10+7*j])
+                    self.glo_ch[sat] = ch
 
         return 0
 
diff --git a/src/cssrlib/rtcm.py b/src/cssrlib/rtcm.py
index ff2cc7b..7a02f4d 100644
--- a/src/cssrlib/rtcm.py
+++ b/src/cssrlib/rtcm.py
@@ -9,10 +9,11 @@
 import numpy as np
 import struct as st
 import bitstruct as bs
-from cssrlib.cssrlib import cssr, sCSSR, sCSSRTYPE, prn2sat, sCType
+from cssrlib.cssrlib import cssr, sCSSR, prn2sat, sCType
 from cssrlib.gnss import uGNSS, sat2id, gpst2time, timediff, time2str, sat2prn
 from cssrlib.gnss import uTYP, uSIG, rSigRnx, bdt2time, bdt2gpst, glo2time
 from cssrlib.gnss import time2bdt, gpst2bdt, rCST, time2gpst, utc2gpst, timeadd
+from cssrlib.gnss import gtime_t
 from crccheck.crc import Crc24LteA
 from enum import IntEnum
 from cssrlib.gnss import Eph, Obs, Geph, Seph
@@ -54,12 +55,22 @@ def __init__(self, foutname=None):
             uGNSS.QZS: 1246, uGNSS.SBS: 1252, uGNSS.BDS: 1258
         }
 
+        self.eph_t = {
+            uGNSS.GPS: 1019, uGNSS.GLO: 1020, uGNSS.BDS: 1042,
+            uGNSS.QZS: 1044, uGNSS.GAL: 1046
+        }
+
     def is_msmtype(self, msgtype):
         for sys_ in self.msm_t.keys():
             if msgtype >= self.msm_t[sys_] and msgtype <= self.msm_t[sys_]+6:
                 return True
         return False
 
+    def adjustweek(self, week: int, tref: gtime_t):
+        week_, _ = time2gpst(tref)
+        week_ref = (week_//1024)*1024
+        return (week % 1024) + week_ref
+
     def msmtype(self, msgtype):
         sys = uGNSS.NONE
         msm = 0
@@ -80,6 +91,14 @@ def ssrtype(self, msgtype):
                 break
         return sys, ssr
 
+    def svid2sat(self, sys, svid):
+        prn = svid
+        if sys == uGNSS.QZS:
+            prn += 192
+        elif sys == uGNSS.SBS:
+            prn += 119
+        return prn2sat(sys, prn)
+
     def ssig2rsig(self, sys: uGNSS, utyp: uTYP, ssig):
         gps_tbl = {
             0: uSIG.L1C,
@@ -373,9 +392,10 @@ def decode_sat(self, msg, i, sys=uGNSS.NONE):
             else:
                 blen = 6
 
-        prn = bs.unpack_from('u'+str(blen), msg, i)[0]
+        svid = bs.unpack_from('u'+str(blen), msg, i)[0]
         i += blen
-        sat = prn2sat(sys, prn)
+
+        sat = self.svid2sat(sys, svid)
         return i, sat
 
     def decode_orb_sat(self, msg, i, k, sys=uGNSS.NONE, inet=0):
@@ -461,12 +481,14 @@ def decode_cssr_orb(self, msg, i, inet=0):
             self.lc[0].iode[sat_] = self.iode_n[k]
             self.lc[0].dorb[sat_] = self.dorb_n[k, :]
 
+            self.set_t0(0, sat_, sCType.ORBIT, self.time)
+
         self.nsat_n += nsat
         self.sys_n += [sys]*nsat
         self.sat_n += sat
 
+        self.iodssr_c[sCType.ORBIT] = v['iodssr']
         self.lc[0].cstat |= (1 << sCType.ORBIT)
-        self.lc[0].t0[sCType.ORBIT] = self.time
 
         return i
 
@@ -486,8 +508,10 @@ def decode_cssr_clk(self, msg, i, inet=0):
             i = self.decode_clk_sat(msg, i, k)
             self.lc[0].dclk[sat_] = self.dclk_n[k]
 
+            self.set_t0(0, sat_, sCType.CLOCK, self.time)
+
+        self.iodssr_c[sCType.CLOCK] = v['iodssr']
         self.lc[0].cstat |= (1 << sCType.CLOCK)
-        self.lc[0].t0[sCType.CLOCK] = self.time
 
         return i
 
@@ -505,6 +529,8 @@ def decode_cssr_cbias(self, msg, i, inet=0):
 
         for k in range(nsat):
             i, sat_ = self.decode_sat(msg, i, sys)
+            self.set_t0(0, sat_, sCType.CBIAS, self.time)
+
             nsig = bs.unpack_from('u5', msg, i)[0]
             i += 5
             if sat_ not in self.sat_b:
@@ -515,13 +541,12 @@ def decode_cssr_cbias(self, msg, i, inet=0):
                 sig, cb = bs.unpack_from('u5s14', msg, i)
                 i += 19
 
-                self.lc[0].cbias[sat_][sig] = self.sval(cb, 14, 0.01)
-                # if self.cssrmode == sCSSRTYPE.GAL_HAS_IDD:
-                # work-around for HAS IDD
-                #    self.lc[0].cbias[sat_][sig] *= -1.0
+                rsig = self.ssig2rsig(sys, uTYP.C, sig)
+                self.lc[0].cbias[sat_][rsig] = self.sval(cb, 14, 0.01)
 
+        self.iodssr_c[sCType.CBIAS] = v['iodssr']
         self.lc[0].cstat |= (1 << sCType.CBIAS)
-        self.lc[0].t0[sCType.CBIAS] = self.time
+
         return i
 
     def decode_cssr_pbias(self, msg, i, inet=0):
@@ -533,12 +558,14 @@ def decode_cssr_pbias(self, msg, i, inet=0):
         # if self.iodssr != v['iodssr']:
         #    return -1
 
-        if timediff(self.time, self.lc[inet].t0[sCType.CBIAS]) > 0:
+        if timediff(self.time, self.lc[inet].t0[sCType.PBIAS]) > 0:
             self.sat_b = []
-            self.lc[0].cbias = {}
+            self.lc[0].pbias = {}
 
         for k in range(nsat):
             i, sat_ = self.decode_sat(msg, i, sys)
+            self.set_t0(0, sat_, sCType.PBIAS, self.time)
+
             nsig = bs.unpack_from('u5', msg, i)[0]
             i += 5
             yaw, dyaw = bs.unpack_from('u9s8', msg, i)
@@ -548,19 +575,18 @@ def decode_cssr_pbias(self, msg, i, inet=0):
             i += 17
             if sat_ not in self.sat_b:
                 self.sat_b.append(sat_)
-                self.lc[0].cbias[sat_] = {}
+                self.lc[0].pbias[sat_] = {}
 
             for j in range(nsig):
                 sig, si, wl, ci, pb = bs.unpack_from('u5u1u2u4s20', msg, i)
                 i += 32
 
-                self.lc[0].pbias[sat_][sig] = self.sval(pb, 20, 1e-4)
-                # if self.cssrmode == sCSSRTYPE.GAL_HAS_IDD:
-                # work-around for HAS IDD
-                #    self.lc[0].pbias[sat_][sig] *= -1.0
+                rsig = self.ssig2rsig(sys, uTYP.L, sig)
+                self.lc[0].pbias[sat_][rsig] = self.sval(pb, 20, 1e-4)
 
+        self.iodssr_c[sCType.PBIAS] = v['iodssr']
         self.lc[0].cstat |= (1 << sCType.PBIAS)
-        self.lc[0].t0[sCType.PBIAS] = self.time
+
         return i
 
     def decode_cssr_comb(self, msg, i, inet=0):
@@ -584,6 +610,9 @@ def decode_cssr_comb(self, msg, i, inet=0):
         sat = []
         for k in range(nsat):
             i, sat_ = self.decode_sat(msg, i, sys)
+            self.set_t0(0, sat_, sCType.ORBIT, self.time)
+            self.set_t0(0, sat_, sCType.CLOCK, self.time)
+
             i = self.decode_orb_sat(msg, i, k, sys)
             i = self.decode_clk_sat(msg, i, k)
             sat.append(sat_)
@@ -595,10 +624,11 @@ def decode_cssr_comb(self, msg, i, inet=0):
         self.sys_n += [sys]*nsat
         self.sat_n += sat
 
+        self.iodssr_c[sCType.ORBIT] = v['iodssr']
+        self.iodssr_c[sCType.CLOCK] = v['iodssr']
+
         self.lc[0].cstat |= (1 << sCType.ORBIT)
-        self.lc[0].t0[sCType.ORBIT] = self.time
         self.lc[0].cstat |= (1 << sCType.CLOCK)
-        self.lc[0].t0[sCType.CLOCK] = self.time
 
         return i
 
@@ -611,13 +641,14 @@ def decode_cssr_ura(self, msg, i):
             self.ura = np.zeros(self.nsat_n)
 
         for k in range(nsat):
-            i, sat = self.decode_sat(msg, i, sys)
-            s = self.sat_n.index(sat)
+            i, sat_ = self.decode_sat(msg, i, sys)
+            s = self.sat_n.index(sat_)
             cls_, val = bs.unpack_from_dict('u3u3', msg, i)
             self.ura[s] = self.quality_idx(cls_, val)
+            self.set_t0(0, sat_, sCType.URA, self.time)
 
+        self.iodssr_c[sCType.URA] = v['iodssr']
         self.lc[0].cstat |= (1 << sCType.URA)
-        self.lc[0].t0[sCType.URA] = self.time
 
         return i
 
@@ -634,9 +665,10 @@ def decode_cssr_hclk(self, msg, i, inet=0):
             i, sat_ = self.decode_sat(msg, i, sys)
             i = self.decode_hclk_sat(msg, i, k)
             self.lc[0].dclk[sat_] = self.dclk_n[k]
+            self.set_t0(0, sat_, sCType.HCLOCK, self.time)
 
+        self.iodssr_c[sCType.HCLOCK] = v['iodssr']
         self.lc[0].cstat |= (1 << sCType.HCLOCK)
-        self.lc[0].t0[sCType.HCLOCK] = self.time
 
         return i
 
@@ -667,8 +699,10 @@ def decode_vtec(self, msg, i, inet=0):
                     i += 16
                     s[l_, j] = self.sval(s_, 16, 5e-3)
 
+        self.iodssr_c[sCType.VTEC] = v['iodssr']
         self.lc[0].cstat |= (1 << sCType.VTEC)
         self.lc[0].t0[sCType.VTEC] = self.time
+        self.set_t0(ctype=sCType.VTEC, t=self.time)
 
         return i
 
@@ -748,10 +782,9 @@ def decode_nrtk_residual(self, msg, i=0):
 
         for k in range(self.nsat):
             sz = 6 if sys != uGNSS.QZS else 4
-            prn = bs.unpack_from('u'+str(sz), msg, i)[0]
-            if sys == uGNSS.QZS:
-                prn += 192
-            sat = prn2sat(sys, prn)
+            svid = bs.unpack_from('u'+str(sz), msg, i)[0]
+            sat = self.svid2sat(sys, svid)
+
             i += sz
             s0c, s0d, s0h, sic, sid = bs.unpack_from('u8u9u6u10u10', msg, i)
             i += 43
@@ -791,21 +824,17 @@ def out_log(self, obs=None, eph=None, geph=None, seph=None):
 
         if self.subtype == sCSSR.CLOCK:
             self.fh.write(" {:s}\t{:s}\n".format("SatID", "dclk [m]"))
-            for k, sat_ in enumerate(self.sat_n):
+            for k, sat_ in enumerate(self.lc[0].dclk.keys()):
                 sys_, _ = sat2prn(sat_)
-                if sys != uGNSS.NONE and sys_ != sys:
-                    continue
                 self.fh.write(" {:s}\t{:5.3f}\n".format(sat2id(sat_),
                                                         self.lc[0].dclk[sat_]))
 
-        elif self.subtype == sCSSR.ORBIT:
+        if self.subtype == sCSSR.ORBIT:
             self.fh.write(" {:s}\t{:s}\t{:s}\t{:s}\t{:s}\n"
                           .format("SatID", "IODE", "Radial[m]",
                                   "Along[m]", "Cross[m]"))
-            for k, sat_ in enumerate(self.sat_n):
+            for k, sat_ in enumerate(self.lc[0].dorb.keys()):
                 sys_, _ = sat2prn(sat_)
-                if sys != uGNSS.NONE and sys_ != sys:
-                    continue
                 self.fh.write(" {:s}\t{:3d}\t{:5.3f}\t{:5.3f}\t{:5.3f}\n".
                               format(sat2id(sat_),
                                      self.lc[0].iode[sat_],
@@ -813,14 +842,12 @@ def out_log(self, obs=None, eph=None, geph=None, seph=None):
                                      self.lc[0].dorb[sat_][1],
                                      self.lc[0].dorb[sat_][2]))
 
-        elif self.subtype == sCSSR.COMBINED:
+        if self.subtype == sCSSR.COMBINED:
             self.fh.write(" {:s}\t{:s}\t{:s}\t{:s}\t{:s}\t{:s}\n"
                           .format("SatID", "IODE", "Radial[m]",
                                   "Along[m]", "Cross[m]", "dclk[m]"))
-            for k, sat_ in enumerate(self.sat_n):
+            for k, sat_ in enumerate(self.lc[0].dorb.keys()):
                 sys_, _ = sat2prn(sat_)
-                if sys != uGNSS.NONE and sys_ != sys:
-                    continue
                 self.fh.write(
                     " {:s}\t{:3d}\t{:6.3f}\t{:6.3f}\t{:6.3f}\t{:6.3f}\n".
                     format(sat2id(sat_),
@@ -830,22 +857,33 @@ def out_log(self, obs=None, eph=None, geph=None, seph=None):
                            self.lc[0].dorb[sat_][2],
                            self.lc[0].dclk[sat_]))
 
-        elif self.subtype == sCSSR.CBIAS:
+        if self.subtype == sCSSR.CBIAS or self.subtype == sCSSR.BIAS:
             self.fh.write(" {:s}\t{:s}\t{:s}\t{:s}\n"
-                          .format("SatID", "SigID", "Bias[m]", "..."))
-            for k, sat_ in enumerate(self.sat_b):
+                          .format("SatID", "SigID", "CBias[m]", "..."))
+            for k, sat_ in enumerate(self.lc[0].cbias.keys()):
                 sys_, _ = sat2prn(sat_)
-                if sys != uGNSS.NONE and sys_ != sys:
-                    continue
                 self.fh.write(" {:s}\t".format(sat2id(sat_)))
                 for j, sig in enumerate(self.lc[0].cbias[sat_].keys()):
-                    sig_ = self.ssig2rsig(sys, uTYP.C, sig)
+                    sig_ = self.ssig2rsig(sys_, uTYP.C, sig)
                     self.fh.write(
                         "{:s}\t{:5.2f}\t".format(sig_.str(),
                                                  self.lc[0].cbias[sat_][sig]))
                 self.fh.write("\n")
 
-        elif self.subtype == sRTCM.ANT_DESC:
+        if self.subtype == sCSSR.PBIAS or self.subtype == sCSSR.BIAS:
+            self.fh.write(" {:s}\t{:s}\t{:s}\t{:s}\n"
+                          .format("SatID", "SigID", "PBias[m]", "..."))
+            for k, sat_ in enumerate(self.lc[0].pbias.keys()):
+                sys_, _ = sat2prn(sat_)
+                self.fh.write(" {:s}\t".format(sat2id(sat_)))
+                for j, sig in enumerate(self.lc[0].pbias[sat_].keys()):
+                    sig_ = self.ssig2rsig(sys_, uTYP.L, sig)
+                    self.fh.write(
+                        "{:s}\t{:5.2f}\t".format(sig_.str(),
+                                                 self.lc[0].pbias[sat_][sig]))
+                self.fh.write("\n")
+
+        if self.subtype == sRTCM.ANT_DESC:
             self.fh.write(" {:20s}{:6d}\n".format("StationID:", self.refid))
             self.fh.write(" {:20s}{:s}\n".format("Antenna Descriptor]:",
                                                  self.ant_desc))
@@ -862,7 +900,7 @@ def out_log(self, obs=None, eph=None, geph=None, seph=None):
                 self.fh.write(" {:20s}{:s}\n".format("Receiver Serial:",
                                                      self.rcv_serial))
 
-        elif self.subtype == sRTCM.NRTK_RES:
+        if self.subtype == sRTCM.NRTK_RES:
             sys, nrtk = self.nrtktype(self.msgtype)
             self.fh.write(" {:20s}{:6d} ({:s})\n".format("NRTK Residual:",
                                                          self.msgtype,
@@ -875,7 +913,8 @@ def out_log(self, obs=None, eph=None, geph=None, seph=None):
             self.fh.write(" {:20s}{:5d}\n".format("Number of satellites:",
                                                   self.nsat))
 
-            self.fh.write(" PRN  s0c[mm] s0d[ppm] s0h[ppm]  sic[mm] sid[ppm]\n")
+            self.fh.write(
+                " PRN  s0c[mm] s0d[ppm] s0h[ppm]  sic[mm] sid[ppm]\n")
             sat_ = self.nrtk_r.keys()
             for sat in sat_:
                 sys_, _ = sat2prn(sat)
@@ -886,7 +925,7 @@ def out_log(self, obs=None, eph=None, geph=None, seph=None):
                 self.fh.write("  {:7.1f} {:8.2f} {:8.1f} {:8.1f} {:8.2f}\n"
                               .format(v[0]*1e3, v[1], v[2], v[3]*1e3, v[4]))
 
-        elif self.subtype == sRTCM.ANT_POS:
+        if self.subtype == sRTCM.ANT_POS:
             self.fh.write(" {:20s}{:6d}\n".format("StationID:", self.refid))
             if self.msgtype in (1005, 1006):
                 self.fh.write(" {:20s}{:6d}\n".format("Station Indicator:",
@@ -896,7 +935,7 @@ def out_log(self, obs=None, eph=None, geph=None, seph=None):
                                   self.pos_arp[0], self.pos_arp[1],
                                   self.pos_arp[2]))
 
-        elif self.subtype == sRTCM.GLO_BIAS:
+        if self.subtype == sRTCM.GLO_BIAS:
             self.fh.write(" {:20s}{:6d}\n".format("StationID:", self.refid))
             self.fh.write(" {:20s} {:8.2f}\n"
                           .format("GLONASS BIAS L1 C/A [m]:",
@@ -911,7 +950,7 @@ def out_log(self, obs=None, eph=None, geph=None, seph=None):
                           .format("GLONASS BIAS L2 P   [m]:",
                                   self.glo_bias[3]))
 
-        elif self.subtype == sRTCM.MSM:
+        if self.subtype == sRTCM.MSM:
             sys, msm = self.msmtype(self.msgtype)
             self.fh.write(" {:20s}{:6d} ({:s})\n".format("MSM:", msm,
                                                          self.sys2str(sys)))
@@ -1197,7 +1236,7 @@ def decode_glo_bias(self, msg, i):
 
     def decode_gps_eph(self, msg, i):
         """ GPS Satellite Ephemeris Message """
-        prn, week, sva, code, idot = bs.unpack_from('u6u10s22s16s8', msg, i)
+        prn, week, sva, code, idot = bs.unpack_from('u6u10u4u2s14', msg, i)
         i += 6+10+4+2+14
         iode, toc, af2, af1, af0 = bs.unpack_from('u8u16s8s16s22', msg, i)
         i += 8+16+8+16+22
@@ -1218,7 +1257,7 @@ def decode_gps_eph(self, msg, i):
         eph.idot = idot*rCST.P2_43*rCST.SC2RAD
         eph.iode = iode
         toc *= 16.0
-        eph.af2 = af0*rCST.P2_55
+        eph.af2 = af2*rCST.P2_55
         eph.af1 = af1*rCST.P2_43
         eph.af0 = af0*rCST.P2_31
         eph.iodc = iodc
@@ -1229,7 +1268,7 @@ def decode_gps_eph(self, msg, i):
         eph.e = e*rCST.P2_33
         eph.cus = cus*rCST.P2_29
         sqrtA = Asq*rCST.P2_19
-        eph.toes = toe*60.0
+        eph.toes = toe*16.0
         eph.cic = cic*rCST.P2_29
         eph.OMG0 = Omg0*rCST.P2_31*rCST.SC2RAD
         eph.cis = cis*rCST.P2_29
@@ -1242,6 +1281,8 @@ def decode_gps_eph(self, msg, i):
         eph.flag = flag
         eph.fit = fit
 
+        eph.week = self.adjustweek(eph.week, self.time)
+
         eph.toe = gpst2time(eph.week, eph.toes)
         eph.toc = gpst2time(eph.week, toc)
         eph.ttr = self.time
@@ -1271,6 +1312,10 @@ def decode_glo_eph(self, msg, i):
         i += 102
 
         geph = Geph()
+        geph.pos = np.zeros(3)
+        geph.vel = np.zeros(3)
+        geph.acc = np.zeros(3)
+
         geph.frq = frq-7
         tk_s *= 30.0
         geph.vel[0] = xd*rCST.P2_20*1e3
@@ -1308,6 +1353,9 @@ def decode_glo_eph(self, msg, i):
         geph.toe = utc2gpst(gpst2time(week, tow+toe))
         geph.mode = 0
 
+        # b7-8: M, b6: P4, b5: P3, b4: P2, b2-3: P1, b0-1: P
+        geph.status = M << 7 | P4 << 6 | P3 << 5 | P2 << 4 | P1 << 2 | P
+
         return i, geph
 
     def decode_gal_eph(self, msg, i):
@@ -1339,7 +1387,7 @@ def decode_gal_eph(self, msg, i):
         eph.idot = idot*rCST.P2_43*rCST.SC2RAD
         toc *= 60.0
 
-        eph.af2 = af0*rCST.P2_59
+        eph.af2 = af2*rCST.P2_59
         eph.af1 = af1*rCST.P2_46
         eph.af0 = af0*rCST.P2_34
         eph.crs = crs*rCST.P2_5
@@ -1361,6 +1409,8 @@ def decode_gal_eph(self, msg, i):
         if self.msgtype == 1046:
             eph.tgd_b = tgd2*rCST.P2_32
 
+        eph.week = self.adjustweek(eph.week, self.time)
+
         eph.toe = gpst2time(eph.week, eph.toes)
         eph.toc = gpst2time(eph.week, toc)
         eph.ttr = self.time
@@ -1377,7 +1427,7 @@ def decode_gal_eph(self, msg, i):
 
     def decode_qzs_eph(self, msg, i):
         """ QZS Satellite Ephemeris Message """
-        prn, toc, af2, af1, af0 = bs.unpack_from('u4u16s8s16s22', msg, i)
+        svid, toc, af2, af1, af0 = bs.unpack_from('u4u16s8s16s22', msg, i)
         i += 4+16+8+16+22
         iode, crs, deln, M0, cuc = bs.unpack_from('u8s16s16s32s16', msg, i)
         i += 8+16+16+32+16
@@ -1391,9 +1441,9 @@ def decode_qzs_eph(self, msg, i):
         i += 6+8+10+1
 
         eph = Eph()
-        eph.sat = prn2sat(uGNSS.QZS, prn+192)
+        eph.sat = prn2sat(uGNSS.QZS, svid+192)
         toc *= 16.0
-        eph.af2 = af0*rCST.P2_55
+        eph.af2 = af2*rCST.P2_55
         eph.af1 = af1*rCST.P2_43
         eph.af0 = af0*rCST.P2_31
         eph.iode = iode
@@ -1422,6 +1472,8 @@ def decode_qzs_eph(self, msg, i):
         eph.iodc = iodc
         eph.fit = fit
 
+        eph.week = self.adjustweek(eph.week, self.time)
+
         eph.toe = gpst2time(eph.week, eph.toes)
         eph.toc = gpst2time(eph.week, toc)
         eph.ttr = self.time
@@ -1527,6 +1579,8 @@ def decode_irn_eph(self, msg, i):
         eph.OMGd = OMGd*rCST.P2_41*rCST.SC2RAD
         eph.i0 = i0*rCST.P2_31*rCST.SC2RAD
 
+        eph.week = self.adjustweek(eph.week, self.time)
+
         eph.toe = gpst2time(eph.week, eph.toes)
         eph.toc = gpst2time(eph.week, toc)
         eph.ttr = self.time
@@ -1538,7 +1592,7 @@ def decode_irn_eph(self, msg, i):
 
     def decode_sbs_eph(self, msg, i):
         """ SBAS Satellite Ephemeris Message """
-        prn, iodn, toc, ura = bs.unpack_from('u6u8u13u4', msg, i)
+        svid, iodn, toc, ura = bs.unpack_from('u6u8u13u4', msg, i)
         i += 6+8+13+4
         x, y, z = bs.unpack_from('s30s30s25', msg, i)
         i += 30+30+25
@@ -1550,7 +1604,7 @@ def decode_sbs_eph(self, msg, i):
         i += 12+8
 
         seph = Seph()
-        seph.sat = prn2sat(uGNSS.SBS, prn+119)
+        seph.sat = prn2sat(uGNSS.SBS, svid+119)
         seph.iode = iodn
         toc *= 16.0
         seph.ura = ura
diff --git a/src/cssrlib/rtk.py b/src/cssrlib/rtk.py
index df34bc3..ad0f70b 100644
--- a/src/cssrlib/rtk.py
+++ b/src/cssrlib/rtk.py
@@ -3,544 +3,54 @@
 
 """
 
+from cssrlib.pppssr import pppos
 import numpy as np
-
-import cssrlib.gnss as gn
+from copy import deepcopy
 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
-
-VAR_HOLDAMB = 0.001
-
-
-def rtkinit(nav, pos0=np.zeros(3), logfile=None):
-    """ initalize RTK-GNSS parameters """
-    nav.pmode = 1  # 0:static, 1:kinematic
-    nav.monlevel = 1
-
-    nav.na = 3 if nav.pmode == 0 else 6
-    nav.nq = 3 if nav.pmode == 0 else 6
-
-    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.el = np.zeros(gn.uGNSS.MAXSAT)
-
-    # parameter for RTK
-    nav.eratio = [50, 50]
-    nav.err = [0, 0.01, 0.005]/np.sqrt(2)
-    nav.sig_p0 = 30.0
-    nav.sig_v0 = 1.0
-    nav.sig_n0 = 30.0
-    nav.sig_qp = 0.01
-    nav.sig_qv = 1.0
-
-    nav.armode = 1  # 1:continuous,2:instantaneous,3:fix-and-hold
-    nav.elmaskar = np.deg2rad(20)  # elevation mask for AR
-    nav.x[0:3] = pos0
-    nav.x[3:6] = 0.0
-
-    dP = np.diag(nav.P)
-    dP.flags['WRITEABLE'] = True
-    dP[0:3] = nav.sig_p0**2
-    nav.q = np.zeros(nav.nq)
-    if nav.pmode >= 1:  # kinematic
-        dP[3:6] = nav.sig_v0**2
-        nav.q[0:3] = nav.sig_qp**2
-        nav.q[3:6] = nav.sig_qv**2
-    else:
-        nav.q[0:3] = nav.sig_qp**2
-
-    # Logging level
-    #
-    nav.fout = None
-    if logfile is None:
-        nav.monlevel = 0
-    else:
-        nav.fout = open(logfile, 'w')
-
-
-def zdres(nav, obs, rs, dts, svh, rr, rtype=1):
-    """ non-differencial residual """
-    _c = gn.rCST.CLIGHT
-    nf = nav.nf
-    n = len(obs.P)
-    y = np.zeros((n, nf*2))
-    el = np.zeros(n)
-    e = np.zeros((n, 3))
-    rr_ = rr.copy()
-    if nav.tidecorr:
-        pos = gn.ecef2pos(rr_)
-        disp = tidedisp(gn.gpst2utc(obs.t), pos)
-        rr_ += disp
-    pos = gn.ecef2pos(rr_)
-    for i in range(n):
-        sys, _ = gn.sat2prn(obs.sat[i])
-        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, :])
-        if el[i] < nav.elmin:
-            continue
-        r += -_c*dts[i]
-        zhd, _, _ = gn.tropmodel(obs.t, pos, np.deg2rad(90.0), 0.0)
-        mapfh, _ = gn.tropmapf(obs.t, pos, el[i])
-        r += mapfh*zhd
-
-        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):
-            if obs.L[i, f] == 0.0:
-                y[i, f] = 0.0
-            else:
-                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-dantPR[f]
-    return y, e, el
-
-
-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
-
-
-def sysidx(satlist, sys_ref):
-    """ return index of satellites with sys=sys_ref """
-    idx = []
-    for k, sat in enumerate(satlist):
-        sys, _ = gn.sat2prn(sat)
-        if sys == sys_ref:
-            idx.append(k)
-    return idx
-
-
-def IB(s, f, na=3):
-    """ return index of phase ambguity """
-    idx = na+gn.uGNSS.MAXSAT*f+s-1
-    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 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(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)
-
-    for sys in obs.sig.keys():
-        for f in range(0, nf*2):
-            # reference satellite
-            idx = sysidx(sat, sys)
-            if len(idx) > 0:
-                i = idx[np.argmax(el[idx])]
-            for j in idx:
-                if i == j:
-                    continue
-                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]
-                H[nv, 0:3] = -e[i, :]+e[j, :]
-                if f < nf:  # carrier
-                    idx_i = IB(sat[i], f, nav.na)
-                    idx_j = IB(sat[j], f, nav.na)
-                    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
-                    Ri[nv] = varerr(nav, el[i], f)
-                    Rj[nv] = varerr(nav, el[j], f)
-                    nav.vsat[sat[i]-1, f] = 1
-                    nav.vsat[sat[j]-1, f] = 1
-                else:
-                    Ri[nv] = varerr(nav, el[i], f)
-                    Rj[nv] = varerr(nav, el[j], f)
-                nb[b] += 1
-                nv += 1
-            b += 1
-    v = np.resize(v, nv)
-    H = np.resize(H, (nv, nav.nx))
-    R = ddcov(nb, b, Ri, Rj, nv)
-
-    return v, H, R
-
-
-def valpos(nav, v, R, thres=4.0):
-    """ post-fit residual test """
-    nv = len(v)
-    fact = thres**2
-    for i in range(nv):
-        if v[i]**2 <= fact*R[i, i]:
-            continue
-        if nav.monlevel > 1:
-            txt = "{:3d} is large: {:8.4f} ({:8.4f})".format(i, v[i], R[i, i])
-            if nav.fout is None:
-                print(txt)
-            else:
-                nav.fout.write(txt+"\n")
-    return True
-
-
-def ddidx(nav, sat):
-    """ index for SD to DD transformation matrix D """
-    nb = 0
-    n = gn.uGNSS.MAXSAT
-    na = nav.na
-    ix = np.zeros((n, 2), dtype=int)
-    nav.fix = np.zeros((n, nav.nf), dtype=int)
-    for m in range(gn.uGNSS.GNSSMAX):
-        k = na
-        for f in range(nav.nf):
-            for i in range(k, k+n):
-                sat_i = i-k+1
-                sys, _ = gn.sat2prn(sat_i)
-                if (sys != m):
-                    continue
-                if sat_i not in sat or nav.x[i] == 0.0 \
-                   or nav.vsat[sat_i-1, f] == 0:
-                    continue
-                if nav.el[sat_i-1] >= nav.elmaskar:
-                    nav.fix[sat_i-1, f] = 2
-                    break
-                else:
-                    nav.fix[sat_i-1, f] = 1
-            for j in range(k, k+n):
-                sat_j = j-k+1
-                sys, _ = gn.sat2prn(sat_j)
-                if (sys != m):
-                    continue
-                if i == j or sat_j not in sat or nav.x[j] == 0.0 \
-                   or nav.vsat[sat_j-1, f] == 0:
-                    continue
-                if nav.el[sat_j-1] >= nav.elmaskar:
-                    ix[nb, :] = [i, j]
-                    nb += 1
-                    nav.fix[sat_j-1, f] = 2
-            k += n
-    ix = np.resize(ix, (nb, 2))
-    return ix
-
-
-def restamb(nav, bias, nb):
-    """ restore SD ambiguity """
-    nv = 0
-    xa = nav.x.copy()
-    xa[0:nav.na] = nav.xa[0:nav.na]
-
-    for m in range(gn.uGNSS.GNSSMAX):
-        for f in range(nav.nf):
-            n = 0
-            index = []
-            for i in range(gn.uGNSS.MAXSAT):
-                sys, _ = gn.sat2prn(i+1)
-                if sys != m or nav.fix[i, f] != 2:
-                    continue
-                index.append(IB(i+1, f, nav.na))
-                n += 1
-            if n < 2:
-                continue
-            xa[index[0]] = nav.x[index[0]]
-            for i in range(1, n):
-                xa[index[i]] = xa[index[0]]-bias[nv]
-                nv += 1
-    return xa
-
-
-def resamb_lambda(nav, sat):
-    """ resolve integer ambiguity using LAMBDA method """
-    nx = nav.nx
-    na = nav.na
-    xa = np.zeros(na)
-    ix = ddidx(nav, sat)
-    nb = len(ix)
-    if nb <= 0:
-        print("no valid DD")
-        return -1, -1
-
-    # y=D*xc, Qb=D*Qc*D', Qab=Qac*D'
-    y = nav.x[ix[:, 0]]-nav.x[ix[:, 1]]
-    DP = nav.P[ix[:, 0], na:nx]-nav.P[ix[:, 1], na:nx]
-    Qb = DP[:, ix[:, 0]-na]-DP[:, ix[:, 1]-na]
-    Qab = nav.P[0:na, ix[:, 0]]-nav.P[0:na, ix[:, 1]]
-
-    # MLAMBDA ILS
-    b, s = mlambda(y, Qb)
-    if s[0] <= 0.0 or s[1]/s[0] >= nav.thresar:
-        nav.xa = nav.x[0:na].copy()
-        nav.Pa = nav.P[0:na, 0:na].copy()
-        bias = b[:, 0]
-        y -= b[:, 0]
-        K = Qab@np.linalg.inv(Qb)
-        nav.xa -= K@y
-        nav.Pa -= K@Qab.T
-
-        # restore SD ambiguity
-        xa = restamb(nav, bias, nb)
-    else:
-        nb = 0
-
-    return nb, xa
-
-
-def initx(nav, x0, v0, i):
-    """ initialize x and P for index i """
-    nav.x[i] = x0
-    for j in range(nav.nx):
-        nav.P[j, i] = nav.P[i, j] = v0 if i == j else 0
-
-
-def kfupdate(x, P, H, v, R):
-    """ Kalman filter measurement update """
-    PHt = P@H.T
-    S = H@PHt+R
-    K = PHt@np.linalg.inv(S)
-    x += K@v
-    P = P - K@H@P
-
-    return x, P, S
-
-
-def udstate(nav, obs, obsb, iu, ir):
-    """ states propagation for kalman filter """
-    tt = gn.timediff(obs.t, nav.t)
-    ns = len(iu)
-    sys = []
-    sat = obs.sat[iu]
-    for sat_i in obs.sat[iu]:
-        sys_i, _ = gn.sat2prn(sat_i)
-        sys.append(sys_i)
-
-    # pos,vel
-    na = nav.na
-    Phi = np.eye(nav.nx)
-    if nav.na > 3:
-        nav.x[0:3] += tt*nav.x[3:6]
-        Phi[0:3, 3:6] = np.eye(3)*tt
-    nav.P = Phi@nav.P@Phi.T
-    dP = np.diag(nav.P)
-    dP.flags['WRITEABLE'] = True
-    dP[0:nav.nq] += nav.q[0:nav.nq]*tt
-
-    # bias
-    for f in range(nav.nf):
-        # reset phase-bias if instantaneous AR or
-        # expire obs outage counter
-        for i in range(gn.uGNSS.MAXSAT):
-            nav.outc[i, f] += 1
-            reset = (nav.outc[i, f] > nav.maxout)
-            sys_i, _ = gn.sat2prn(i+1)
-            if sys_i not in obs.sig.keys():
-                continue
-            j = IB(i+1, 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
-        for i in range(ns):
-            if sys[i] not in obs.sig.keys():
-                continue
-            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))
-        # bias
-        bias = np.zeros(ns)
-        offset = 0
-        na = 0
-        for i in range(ns):
-            if sys[i] not in obs.sig.keys():
-                continue
-            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/lam
-            amb = nav.x[IB(sat[i], f, nav.na)]
-            if amb != 0.0:
-                offset += bias[i]-amb
-                na += 1
-        # 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
-        for i in range(ns):
-            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)
-    return 0
-
-
-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, _ 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:
-            continue
-        idx_u.append(k)
-
-    # exclude satellite with missing observation and cycle slip for base
-    idx_r = []
-    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 \
-           elb[k] < nav.elmin:
-            continue
-        idx_r.append(k)
-
-    idx = np.intersect1d(obs.sat[idx_u], obsb.sat[idx_r], return_indices=True)
-    k = len(idx[0])
-    iu = np.array(idx_u)[idx[1]]
-    ir = np.array(idx_r)[idx[2]]
-    return k, iu, ir
-
-
-def holdamb(nav, xa):
-    """ hold integer ambiguity """
-    nb = nav.nx-nav.na
-    v = np.zeros(nb)
-    H = np.zeros((nb, nav.nx))
-    nv = 0
-    for m in range(gn.uGNSS.GNSSMAX):
-        for f in range(nav.nf):
-            n = 0
-            index = []
-            for i in range(gn.uGNSS.MAXSAT):
-                sys, _ = gn.sat2prn(i+1)
-                if sys != m or nav.fix[i, f] != 2:
-                    continue
-                index.append(IB(i+1, f, nav.na))
-                n += 1
-                nav.fix[i, f] = 3  # hold
-            # constraint to fixed ambiguity
-            for i in range(1, n):
-                v[nv] = (xa[index[0]]-xa[index[i]]) - \
-                    (nav.x[index[0]]-nav.x[index[i]])
-                H[nv, index[0]] = 1.0
-                H[nv, index[i]] = -1.0
-                nv += 1
-    if nv > 0:
-        R = np.eye(nv)*VAR_HOLDAMB
-        # update states with constraints
-        nav.x, nav.P, _ = kfupdate(nav.x, nav.P, H[0:nv, :], v[0:nv], R)
-    return 0
-
 
-def relpos(nav, obs, obsb):
-    """ relative positioning for RTK-GNSS """
-    nf = nav.nf
-    if gn.timediff(obs.t, obsb.t) != 0:
-        return -1
 
-    rs, _, dts, svh, _ = satposs(obs, nav)
-    rsb, _, dtsb, svhb, _ = satposs(obsb, nav)
+class rtkpos(pppos):
+    """ class for RTK processing """
 
-    # non-differencial residual for base
-    yr, er, elr = zdres(nav, obsb, rsb, dtsb, svhb, nav.rb, 0)
-    ns, iu, ir = selsat(nav, obs, obsb, elr)
-    y = np.zeros((ns*2, nf*2))
-    e = np.zeros((ns*2, 3))
+    def __init__(self, nav, pos0=np.zeros(3), logfile=None):
+        """ initialize variables for PPP-RTK """
 
-    if ns < 4:
-        return -1
+        # trop, iono from cssr
+        # phase windup model is local/regional
+        super().__init__(nav=nav, pos0=pos0, logfile=logfile,
+                         trop_opt=0, iono_opt=0, phw_opt=0)
 
-    y[ns:, :] = yr[ir, :]
-    e[ns:, :] = er[ir, :]
+        self.nav.eratio = np.ones(self.nav.nf)*50  # [-] factor
+        self.nav.err = [0, 0.01, 0.005]/np.sqrt(2)  # [m] sigma
+        self.nav.sig_p0 = 30.0  # [m]
+        self.nav.thresar = 2.0  # AR acceptance threshold
+        self.nav.armode = 1     # AR is enabled
 
-    # Kalman filter time propagation
-    udstate(nav, obs, obsb, iu, ir)
+    def base_process(self, obs, obsb, rs, dts, svh):
+        """ processing for base station in RTK """
+        rsb, vsb, dtsb, svhb, _ = satposs(obsb, self.nav)
+        yr, er, elr = self.zdres(
+            obsb, None, None, rsb, vsb, dtsb, self.nav.rb, 0)
 
-    #xa = np.zeros(nav.nx)
-    xp = nav.x
+        # Editing observations (base/rover)
+        sat_ed_r = self.qcedit(obsb, rsb, dtsb, svhb, rr=self.nav.rb)
+        sat_ed_u = self.qcedit(obs, rs, dts, svh)
 
-    # non-differencial residual for rover
-    yu, eu, el = zdres(nav, obs, rs, dts, svh, xp[0:3])
+        # define common satellite between base and rover
+        sat_ed = np.intersect1d(sat_ed_u, sat_ed_r, True)
+        ir = np.intersect1d(obsb.sat, sat_ed, True, True)[1]
+        iu = np.intersect1d(obs.sat, sat_ed, True, True)[1]
+        ns = len(iu)
 
-    y[:ns, :] = yu[iu, :]
-    e[:ns, :] = eu[iu, :]
-    el = el[iu]
-    sat = obs.sat[iu]
-    nav.el[sat-1] = el
-    # DD residual
-    v, H, R = ddres(nav, obs, xp, y, e, sat, el)
-    Pp = nav.P
+        y = np.zeros((ns*2, self.nav.nf*2))
+        e = np.zeros((ns*2, 3))
 
-    # Kalman filter measurement update
-    xp, Pp, _ = kfupdate(xp, Pp, H, v, R)
+        y[ns:, :] = yr[ir, :]
+        e[ns:, :] = er[ir, :]
 
-    # non-differencial residual for rover after measurement update
-    yu, eu, _ = zdres(nav, obs, rs, dts, svh, xp[0:3])
-    y[:ns, :] = yu[iu, :]
-    e[:ns, :] = eu[iu, :]
-    # residual for float solution
-    v, H, R = ddres(nav, obs, xp, y, e, sat, el)
-    if valpos(nav, v, R):
-        nav.x = xp
-        nav.P = Pp
-    else:
-        nav.smode = 0
+        obs_ = deepcopy(obs)
+        obs_.sat = obs.sat[iu]
+        obs_.L = obs.L[iu, :]-obsb.L[ir, :]
+        obs_.P = obs.P[iu, :]-obsb.P[ir, :]
 
-    nb, xa = resamb_lambda(nav, sat)
-    nav.smode = 5  # float
-    if nb > 0:
-        yu, eu, _ = zdres(nav, obs, rs, dts, svh, xa[0:3])
-        y[:ns, :] = yu[iu, :]
-        e[:ns, :] = eu[iu, :]
-        v, H, R = ddres(nav, obs, xa, y, e, sat, el)
-        if valpos(nav, v, R):
-            if nav.armode == 3:
-                holdamb(nav, xa)
-            nav.smode = 4  # fix
-    nav.t = obs.t
-    return 0
+        return y, e, iu, obs_


Test patch:
diff --git a/src/cssrlib/test/test_eph.py b/src/cssrlib/test/test_eph.py
index bf21628..40fbacc 100644
--- a/src/cssrlib/test/test_eph.py
+++ b/src/cssrlib/test/test_eph.py
@@ -2,7 +2,7 @@
 import cartopy.crs as ccrs
 import numpy as np
 from cssrlib.rinex import rnxdec
-from cssrlib.gnss import Nav, epoch2time, prn2sat, uGNSS, sat2prn,\
+from cssrlib.gnss import Nav, epoch2time, prn2sat, uGNSS, sat2prn, \
     timeadd, ecef2pos
 from cssrlib.ephemeris import findeph, eph2pos
 
diff --git a/src/cssrlib/test/test_gnss.py b/src/cssrlib/test/test_gnss.py
index 5b66a72..6e54b43 100644
--- a/src/cssrlib/test/test_gnss.py
+++ b/src/cssrlib/test/test_gnss.py
@@ -53,5 +53,8 @@
 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())
+try:
+    sig = rSigRnx(uGNSS.GPS, "D1X", "what??")
+    print(sys2char(sig.sys), sig.str())
+except ValueError:
+    print("error in value.")
diff --git a/src/cssrlib/test/test_peph.py b/src/cssrlib/test/test_peph.py
index e62a886..911d11d 100644
--- a/src/cssrlib/test/test_peph.py
+++ b/src/cssrlib/test/test_peph.py
@@ -2,29 +2,30 @@
 Test script for peph module
 """
 import numpy as np
-from os.path import expanduser
 
 from cssrlib.peph import atxdec, biasdec, peph
 from cssrlib.peph import searchpcv, antModelRx, antModelTx
 from cssrlib.rinex import rnxdec
 from cssrlib.gnss import Nav
-from cssrlib.gnss import epoch2time, time2epoch, timeadd
+from cssrlib.gnss import epoch2time, time2epoch, time2doy, timeadd
 from cssrlib.gnss import sat2id, id2sat, sys2char
 from cssrlib.gnss import rSigRnx
 from cssrlib.gnss import pos2ecef, enu2xyz
 
+ep = [2021, 9, 22, 12, 0, 0]
+time = epoch2time(ep)
+doy = int(time2doy(time))
 
-bdir = expanduser('../../../../cssrlib-data/data/')
+bdir = '../../../../cssrlib-data/data/'
+orbfile = bdir+"COD0MGXFIN_{:4d}{:03d}0000_01D_05M_ORB.SP3".format(ep[0], doy)
+clkfile = bdir+"COD0MGXFIN_{:4d}{:03d}0000_01D_30S_CLK.CLK".format(ep[0], doy)
+dcbfile = bdir+"COD0MGXFIN_{:4d}{:03d}0000_01D_01D_OSB.BIA".format(ep[0], doy)
 atxfile = bdir+"igs14.atx"
-orbfile = bdir+"COD0OPSRAP_20210780000_01D_15M_ORB.SP3"
-clkfile = bdir+"COD0OPSRAP_20210780000_01D_30S_CLK.CLK"
-dcbfile = bdir+"COD0OPSRAP_20210780000_01D_01D_OSB.BIA"
 
-time = epoch2time([2021, 3, 19, 12, 0, 0])
 sat = id2sat("G01")
 sig = rSigRnx("GC1C")
 
-if False:
+if True:
 
     print("Test SP3 and Clock-RINEX module")
     print()
@@ -34,16 +35,20 @@
     sp = peph()
 
     nav = sp.parse_sp3(orbfile, nav)
+    sp.write_sp3('test_peph.sp3', nav)
+
     nav = rnx.decode_clk(clkfile, nav)
 
     n = 10
     rs = np.zeros((1, 6))
     dts = np.zeros((1, 2))
     for k in range(n):
+
         t = timeadd(time, 30*k)
+        ep = time2epoch(t)
+
         rs[0, :], dts[0, :], var = sp.peph2pos(t, sat, nav)
 
-        ep = time2epoch(t)
         print("{:4d}-{:02d}-{:02d} {:02d}:{:02d}:{:02d}  {:s}  "
               "{:15.4f} {:15.4f} {:15.4f} {:15.6f}"
               .format(ep[0], ep[1], ep[2], ep[3], ep[4], ep[5], sat2id(sat),
@@ -51,7 +56,7 @@
 
     print()
 
-if True:
+if False:
 
     print("Test ANTEX module")
     print()
diff --git a/src/cssrlib/test/test_rtk.py b/src/cssrlib/test/test_rtk.py
index a4ba3fa..7393a17 100644
--- a/src/cssrlib/test/test_rtk.py
+++ b/src/cssrlib/test/test_rtk.py
@@ -8,7 +8,7 @@
 
 import cssrlib.rinex as rn
 import cssrlib.gnss as gn
-from cssrlib.rtk import rtkinit, relpos
+from cssrlib.rtk import rtkpos
 from cssrlib.gnss import rSigRnx
 from cssrlib.peph import atxdec, searchpcv
 
@@ -25,13 +25,13 @@
 #
 sigs = [rSigRnx("GC1C"), rSigRnx("GC2W"),
         rSigRnx("EC1C"), rSigRnx("EC5Q"),
-        rSigRnx("GL1C"), rSigRnx("GL2W"),
-        rSigRnx("EL1C"), rSigRnx("EL5Q")]
+        rSigRnx("GL1C"), rSigRnx("GL2W"), rSigRnx("GS1C"), rSigRnx("GS2W"),
+        rSigRnx("EL1C"), rSigRnx("EL5Q"), rSigRnx("ES1C"), rSigRnx("ES5Q")]
 
 sigsb = [rSigRnx("GC1C"), rSigRnx("GC2W"),
          rSigRnx("EC1X"), rSigRnx("EC5X"),
-         rSigRnx("GL1C"), rSigRnx("GL2W"),
-         rSigRnx("EL1X"), rSigRnx("EL5X")]
+         rSigRnx("GL1C"), rSigRnx("GL2W"), rSigRnx("GS1C"), rSigRnx("GS2W"),
+         rSigRnx("EL1X"), rSigRnx("EL5X"), rSigRnx("ES1X"), rSigRnx("ES5X")]
 
 # rover
 #
@@ -79,13 +79,14 @@
 enu = np.zeros((nep, 3))
 smode = np.zeros(nep, dtype=int)
 
-rtkinit(nav, dec.pos)
+rtk = rtkpos(nav, dec.pos, 'test_rtk.log')
 rr = dec.pos
+
 for ne in range(nep):
     obs, obsb = rn.sync_obs(dec, decb)
     if ne == 0:
         t0 = nav.t = obs.t
-    relpos(nav, obs, obsb)
+    rtk.process(obs, obsb=obsb)
     t[ne] = gn.timediff(nav.t, t0)
     sol = nav.xa[0:3] if nav.smode == 4 else nav.x[0:3]
     enu[ne, :] = gn.ecef2enu(pos_ref, sol-xyz_ref)