Merge pull request #46 from AndreHauschild/devel

Added workflow for automated testing

.github/workflows/test.yml

@@ -0,0 +1,39 @@
+name: Execute Python Scripts
+
+on:
+  push:
+    branches:
+      - main
+      - devel
+
+jobs:
+  execute-scripts:
+    name: Execute Python Scripts
+    runs-on: ubuntu-latest
+
+    steps:
+      - name: Checkout cssrlib repository
+        uses: actions/checkout@v2
+
+      - name: Set up Python
+        uses: actions/setup-python@v2
+        with:
+          python-version: '3.10'
+
+      - name: Install dependencies from repository
+        run: pip install -r requirements.txt
+
+      - name: Install cssrlib from cloned repository
+        run: pip install .
+
+      - name: Execute Python scripts
+        run: |
+          cd ./src/cssrlib/test
+          for file in $(find . -name "*.py"); do
+            echo "# Execute $file"
+            echo "#"
+            python "$file" > /dev/null
+            echo "# ...done!"
+            echo "#"
+            echo
+          done

README.md

@@ -23,17 +23,10 @@ Click this button for a quick demo in Google Colab
 
 Prerequisites
 =============
-Additional python packages are required as prerequisites and can be installed via the following commands
+Additional python packages are required as prerequisites and can be installed via the following command
 
 ```
-pip install bitstruct galois crccheck pysolid
-pip install notebook numpy matplotlib
-```
-
-Optionally, on linux, users can install the `cartopy` package and the `PySolid` package
-
-```
-pip install cartopy pysolid
+pip install -r requirements.txt
 ```
 
 If the installation of `cartopy` fails, try installing `libgeos++-dev` first.

requirements.txt

@@ -0,0 +1,8 @@
+bitstruct
+galois
+crccheck
+notebook
+numpy
+matplotlib
+cartopy
+pysolid

src/cssrlib/cssrlib.py

@@ -210,8 +210,8 @@ def __init__(self):
         self.inet = -1
         self.inet_ref = -1
         self.ng = -1
-        self.pbias = None
-        self.cbias = None
+        self.pbias = {}
+        self.cbias = {}
         self.iode = None
         self.dorb = None
         self.dclk = None

src/cssrlib/ephemeris.py

@@ -122,9 +122,26 @@ def geph2clk(time: gtime_t, geph: Geph):
     return -geph.taun + geph.gamn*t
 
 
-def eph2pos(t: gtime_t, eph: Eph, flg_v=False):
-    """ calculate satellite position based on ephemeris """
-    sys, prn = sat2prn(eph.sat)
+def geph2rel(rs, vs):
+    return - 2.0*(rs@vs)/(rCST.CLIGHT**2)
+
+
+def eccentricAnomaly(M, e):
+    """
+    Compute eccentric anomaly based on mean anomaly and eccentricity
+    """
+    E = M
+    for _ in range(10):
+        Eold = E
+        sE = np.sin(E)
+        E = M+e*sE
+        if abs(Eold-E) < 1e-12:
+            break
+
+    return E, sE
+
+
+def sys2MuOmega(sys):
     if sys == uGNSS.GAL:
         mu = rCST.MU_GAL
         omge = rCST.OMGE_GAL
@@ -134,6 +151,13 @@ def eph2pos(t: gtime_t, eph: Eph, flg_v=False):
     else:  # GPS,QZS
         mu = rCST.MU_GPS
         omge = rCST.OMGE
+    return mu, omge
+
+
+def eph2pos(t: gtime_t, eph: Eph, flg_v=False):
+    """ calculate satellite position based on ephemeris """
+    sys, prn = sat2prn(eph.sat)
+    mu, omge = sys2MuOmega(sys)
     dt = dtadjust(t, eph.toe)
     n0 = np.sqrt(mu/eph.A**3)
     dna = eph.deln
@@ -143,13 +167,7 @@ def eph2pos(t: gtime_t, eph: Eph, flg_v=False):
         Ak += dt*eph.Adot
     n = n0+dna
     M = eph.M0+n*dt
-    E = M
-    for _ in range(10):
-        Eold = E
-        sE = np.sin(E)
-        E = M+eph.e*sE
-        if abs(Eold-E) < 1e-12:
-            break
+    E, sE = eccentricAnomaly(M, eph.e)
     cE = np.cos(E)
     dtc = dtadjust(t, eph.toc)
     dtrel = -2.0*np.sqrt(mu*eph.A)*eph.e*sE/rCST.CLIGHT**2
@@ -223,13 +241,29 @@ def eph2clk(time, eph):
     return dts
 
 
+def eph2rel(time, eph):
+    sys, _ = sat2prn(eph.sat)
+    mu, _ = sys2MuOmega(sys)
+    dt = dtadjust(time, eph.toe)
+    n0 = np.sqrt(mu/eph.A**3)
+    dna = eph.deln
+    Ak = eph.A
+    if eph.mode > 0:
+        dna += 0.5*dt*eph.delnd
+        Ak += dt*eph.Adot
+    n = n0+dna
+    M = eph.M0+n*dt
+    _, sE = eccentricAnomaly(M, eph.e)
+    mu, _ = sys2MuOmega(sys)
+    return -2.0*np.sqrt(mu*eph.A)*eph.e*sE/rCST.CLIGHT**2
+
+
 def satpos(sat, t, nav, cs=None, orb=None):
     """
     Calculate pos/vel/clk for single satellite
 
     The satellite position, velocity and clock offset are computed at epoch.
-    The satellite clock is already corrected for the relativistic effects. The
-    satellite health indicator is extracted from the broadcast navigation
+    The satellite health indicator is extracted from the broadcast navigation
     message.
 
     Parameters
@@ -272,15 +306,33 @@ def satpos(sat, t, nav, cs=None, orb=None):
         rs_, dts_, _ = orb.peph2pos(t, sat, nav)
         if rs_ is None or dts_ is None or np.isnan(dts_[0]):
             return rs, vs, dts, svh
-        dt = dts_[0]
 
-        if len(nav.eph) > 0:
+        # Health indicator from BRDC
+        #
+
+        if sys == uGNSS.GLO and len(nav.geph) > 0:
+
+            geph = findeph(nav.geph, t, sat)
+            if geph is None:
+                svh[i] = 1
+                return rs, vs, dts, svh
+
+            svh[i] = geph.svh
+
+            if sat not in nav.glo_ch:
+                nav.glo_ch[sat] = geph.frq
+
+        elif len(nav.eph) > 0:
+
             eph = findeph(nav.eph, t, sat)
             if eph is None:
                 svh[i] = 1
                 return rs, vs, dts, svh
+
             svh[i] = eph.svh
+
         else:
+
             svh[i] = 0
 
     else:
@@ -355,41 +407,41 @@ def satpos(sat, t, nav, cs=None, orb=None):
             mode = 0
 
         if sys == uGNSS.GLO:
+
             geph = findeph(nav.geph, t, sat, iode, mode=mode)
             if geph is None:
                 svh[i] = 1
                 return rs, vs, dts, svh
 
             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
                 return rs, vs, dts, svh
 
             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]
-        dts[i] = dts_[0]
+        dts[i] = dts_[0] - orb.pephrel(rs_)  # Remove relativistic correction!
 
     else:
 
         if sys == uGNSS.GLO:
             rs[i, :], vs[i, :], dts[i] = geph2pos(t, geph, True)
+            dts[i] -= geph2rel(rs[i, :], vs[i, :])
         else:
             rs[i, :], vs[i, :], dts[i] = eph2pos(t, eph, True)
+            dts[i] -= eph2rel(t, eph)
 
         # Apply SSR correction
         #

src/cssrlib/test/test_peph.py

@@ -16,11 +16,11 @@
 time = epoch2time(ep)
 doy = int(time2doy(time))
 
-bdir = '../../../../cssrlib-data/data/'
+bdir = '../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"
+atxfile = bdir+"test.atx"
 
 sat = id2sat("G01")
 sig = rSigRnx("GC1C")
@@ -56,7 +56,7 @@
 
     print()
 
-if False:
+if True:
 
     print("Test ANTEX module")
     print()
@@ -175,7 +175,7 @@
 
     print()
 
-if False:
+if True:
 
     print("Test Bias-SINEX module")
     print()
@@ -184,20 +184,19 @@
     bd.parse(dcbfile)
 
     sat = id2sat("G03")
-    sig = rSigRnx("GC1W")
+    sig = rSigRnx("GC1C")
 
-    bias, std, = bd.getosb(sat, time, sig)
-    assert bias == 7.6934
-    assert std == 0.0
+    bias = bd.getosb(sat, time, sig)
 
-    print("{:s} {:s} {:8.5f} {:6.4f}"
-          .format(sat2id(sat), sig.str(), bias, std))
+    print("{:s} {:s} {:8.5f}"
+          .format(sat2id(sat), sig.str(), bias))
+    assert bias == -0.7324
 
     sig = rSigRnx("GL1W")
-    bias, std, = bd.getosb(sat, time, sig)
-    assert bias == 0.00038
-    assert std == 0.0
+    bias = bd.getosb(sat, time, sig)
+
+    print("{:s} {:s} {:8.5f}"
+          .format(sat2id(sat), sig.str(), bias))
+    assert bias == 0.67184
 
-    print("{:s} {:s} {:8.5f} {:6.4f}"
-          .format(sat2id(sat), sig.str(), bias, std))
     print()