[Python] GPS satellite position calculation from RINEX navigation file produces unexpected...

For a college project I have to calculate the position of GPS satellites for a period of time with navigation data. There may be or 4 or 6 hours in a given navigation file - it depends on number of epochs for a specific GPS. The more epoch the more time. I must calculate the GPS position over 6 hours with a navigation file.

I have to compute GPS satellite positions for a period of 6 hours using a RINEX 3.02 navigation file.
I wrote code to parse the navigation data and compute satellite positions in both ECI and ECEF frames.

After days of working to read and save rinex navigation files I finally managed to calculate satellite position, and plotted positions (x, y, z in both eci and ecef).

However, when I plot the positions over time, the satellite trajectory looks incorrect or chaotic, like this:



What I expect is a smooth orbital path like this:



Some key details:

• 
  I compute tk = t - toe, where toe is the time of ephemeris and t is seconds of week for the epoch of interest.
  
  
• 
  I attempted using the nearest ephemeris for each time epoch.
  
  
• 
  When I fix the navigation data to only one epoch and calculate for all times using that, the shape looks better (still not perfect).
  

My main issue is that if I compute the satellite positions using only one ephemeris epoch (within about ±2 hours of its reference time), the results are correct. This makes sense because the navigation parameters are valid for roughly that time span.

However, for a 6‑hour time span or more, I need to use multiple ephemeris epochs. My approach has been: for each computation time (in my code I store this as te — time of emission), I first find the nearest ephemeris epoch to te and then use that epoch’s navigation parameters.

But when I do this, the trajectory becomes distorted (like the first plot I uploaded).

Is there a correct way to handle switching between multiple ephemeris epochs, or a standard method to avoid these distortions?

What could be causing such distortion? Is it a mistake in handling tk, or do I need to interpolate parameters more carefully?

Here's the core part of my code:

def calc_tk(t, toe):
tk = t - toe
if tk > 302400:
tk -= 604800
elif tk < -302400:
tk += 604800
return tk

# satellite position calculation here
toe_list = np.array(list(nav_data[prn].keys()))
toe = toe_list[np.argmin(np.abs(te - toe_list))] #using this calculate nearest epoch


First two ephemeris records from my RINEX file:

Direct link to full RINEX file:
https://github.com/shahriar-nzz/shahriar-repository/blob/main/GODS00USA_R_20240010000_01D_GN.rnx

3.04 N: GNSS NAV DATA G: GPS RINEX VERSION / TYPE
JPS2RIN v.2.1.223 JAVAD GNSS 20240102 000924 UTC PGM / RUN BY / DATE
18 LEAP SECONDS
GODS MARKER NAME COMMENT
40451M128 MARKER NUMBER COMMENT
1130752.1541 -4831349.1034 3994098.9626 COMMENT
This data is provided as a public service by NASA/JPL. COMMENT
No warranty is expressed or implied regarding suitability COMMENT
for use. For further information, contact: COMMENT
ggnops at jpl dot nasa dot gov COMMENT
END OF HEADER
G07 2024 01 01 01 59 44-2.613384276628D-05-9.436007530894D-12 0.000000000000D+00
4.400000000000D+01 1.200000000000D+01 4.854845080927D-09 6.239840220677D-01
8.177012205124D-07 1.773746171966D-02 8.910894393921D-06 5.153683597565D+03
9.358400000000D+04-2.365559339523D-07 1.512483407131D+00 1.322478055954D-07
9.498578426946D-01 2.058125000000D+02-2.159927596900D+00-8.354633718315D-09
2.221521106700D-10 1.000000000000D+00 2.295000000000D+03 0.000000000000D+00
2.000000000000D+00 0.000000000000D+00-1.117587000000D-08 4.400000000000D+01
8.649000000000D+04 4.000000000000D+00
G01 2024 01 01 02 00 00 1.650028862059D-04 7.958078640513D-13 0.000000000000D+00
7.700000000000D+01 7.565625000000D+01 3.961593587703D-09 1.254732735329D+00
3.958120942116D-06 1.305459404830D-02 6.889924407005D-06 5.154026269913D+03
9.360000000000D+04-2.197921276093D-07-1.609739761974D+00 1.154839992523D-07
9.903301560032D-01 2.633125000000D+02 9.986915388941D-01-8.050692486514D-09
1.357199389945D-10 1.000000000000D+00 2.295000000000D+03 0.000000000000D+00
2.000000000000D+00 6.300000000000D+01 5.122274000000D-09 7.700000000000D+01
8.649000000000D+04 4.000000000000D+00


My code (full repo):

https://github.com/shahriar-nzz/shahriar-repository/tree/main

Continue reading...