- optimized by eliminating frequent update of sun position, antTx position. (suggested by mfinneman)

Author: hirokawa

src/cssrlib/gnss.py

@@ -9,7 +9,7 @@
 from datetime import datetime, timezone
 import bitstruct as bs
 import yaml
-
+from numba import njit
 
 gpst0 = [1980, 1, 6, 0, 0, 0]  # GPS system time reference
 gst0 = [1999, 8, 22, 0, 0, 0]  # Galileo system time reference
@@ -951,6 +951,9 @@ def __init__(self, nf=2):
 
         # number of satellite (observed, calculated, corrected)
         self.nsat = [0, 0, 0]
+        
+        self.t_rsun = None # reference time of rsun
+        self.rsun = None # placeholder of sun position
 
 
 def epoch2time(ep):
@@ -1208,7 +1211,7 @@ def prn2sat(sys, prn):
         sat = 0
     return sat
 
-
+@njit
 def sat2prn(sat):
     """ convert sat to sys+prn """
     if sat > uGNSS.MAXSAT:

src/cssrlib/peph.py

@@ -795,7 +795,7 @@ def antModelTx(nav, e, sigs, sat, time, rs, sig0=None):
 
     # Rotation matrix from satellite antenna frame to ECEF frame [ex, ey, ez]
     #
-    A = orb2ecef(time, rs)
+    A = orb2ecef(time, rs, nav.rsun)
     ez = A[2, :]
 
     # Zenith angle and zenith angle grid
@@ -1124,14 +1124,14 @@ def utc2gmst(t: gtime_t, ut1_utc):
     return np.fmod(gmst, 86400.0)*np.pi/43200.0  # 0 <= gmst <= 2*PI
 
 
-def orb2ecef(time, rs):
+def orb2ecef(time, rs, rsun=None):
     """
     Rotation matrix from satellite antenna frame to ECEF frame assuming
     standard yaw attitude law
     """
 
-    erpv = np.zeros(5)
-    rsun, _, _ = sunmoonpos(gpst2utc(time), erpv, True)
+    if rsun is None:
+        rsun, _, _ = sunmoonpos(gpst2utc(time), np.zeros(5), True)
     r = -rs
     ez = r/np.linalg.norm(r)
     r = rsun-rs

src/cssrlib/ppp.py

@@ -274,13 +274,12 @@ def shapiro(rsat, rrcv):
     return corr
 
 
-def windupcorr(time, rs, vs, rr, phw, full=False):
+def windupcorr(time, rs, vs, rr, rsun, phw, full=False):
     """ calculate windup correction """
     ek = vnorm(rr-rs)
     if full:
         # Satellite antenna frame unit vectors assuming standard yaw attitude law
         #
-        rsun, _, _ = sunmoonpos(gpst2utc(time))
         r = -rs
         ezs = r/np.linalg.norm(r)
         r = rsun-rs
@@ -443,7 +442,8 @@ def tidedispIERS2010(tutc, pos, erpv=None):
             tn = timeadd(tgps_, k*300)
             eph = findeph(nav.eph, tn, sat)
             rs_, vs_, dts = eph2pos(tn, eph, True)
-            phw_ = windupcorr(tn, rs_, vs_, rr_, phw_)
+            rsun_, _, _ = sunmoonpos(gpst2utc(tn))
+            phw_ = windupcorr(tn, rs_, vs_, rr_, rsun_, phw_)
             t[k] = timediff(tn, tgps_)
             ph_[k] = phw_
             d[k] = shapiro(rs_, rr_)

src/cssrlib/pppssr.py

@@ -9,7 +9,7 @@
 from cssrlib.gnss import ecef2pos, tropmodel, geodist, satazel, uTideModel
 from cssrlib.gnss import time2str, timediff, gpst2utc, tropmapf
 from cssrlib.gnss import trop_model_tbl, iono_model_tbl, default_config
-from cssrlib.ppp import tidedisp, tidedispIERS2010, shapiro, windupcorr
+from cssrlib.ppp import tidedisp, tidedispIERS2010, shapiro, windupcorr, sunmoonpos
 from cssrlib.peph import antModelRx, antModelTx
 from cssrlib.cssrlib import sCType
 from cssrlib.cssrlib import sCSSRTYPE as sc
@@ -656,7 +656,8 @@ def zdres(self, obs, cs, bsx, rs, vs, dts, rr, rtype=1):
             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],
+                                                 rr_, self.nav.rsun, 
+                                                 self.nav.phw[sat-1],
                                                  full=phw_mode)
 
                 # cycle -> m
@@ -738,11 +739,11 @@ def zdres(self, obs, cs, bsx, rs, vs, dts, rr, rtype=1):
                                                     sc.RTCM3_SSR,
                                                     sc.BDS_PPP,
                                                     sc.PVS_PPP):
-
-                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)
+                
+                ants = antModelTx(self.nav, e[i, :], (sigsPR + sigsCP),
+                                  sat, obs.t, rs[i, :], sig0)                
+                antsPR = ants[0:len(sigsPR)]
+                antsCP = ants[len(sigsPR):]
 
             else:
 
@@ -1390,6 +1391,13 @@ def process(self, obs, cs=None, orb=None, bsx=None, obsb=None):
 
         self.nav.nsat[0] = len(obs.sat)
 
+        # sun position is updated every 1sec
+        if (self.nav.rsun is None) or \
+            (np.abs(timediff(self.nav.t_rsun,obs.t))>=1.0):
+            self.nav.rsun, _, _ = sunmoonpos(gpst2utc(obs.t))
+            self.nav.t_rsun = obs.t
+                
+
         # GNSS satellite positions, velocities and clock offsets
         # for all satellite in RINEX observations
         #