feat: add --use-obs option to limit skyplot to receiver observation period

utils/skyplot/README.md

@@ -16,12 +16,14 @@ showing satellite visibility over time.
 ## Features
 
 - Processes RINEX navigation files.
+- Optionally uses an OBS file to limit plot to the receiver active pricessing time (--use-obs).
+  - When enabled, the tool looks for a matching file by replacing the last character of the NAV filename with O/o and uses it if found.
 - Calculates satellite positions using broadcast ephemeris.
 - Plots satellite tracks in azimuth-elevation coordinates.
 - Color-codes satellites by constellation (GPS, Galileo, GLONASS, BeiDou).
 - Customizable observer location.
 - Outputs high-quality image in PDF format.
-- Non-interative mode for CI jobs (with `--no-show` flag).
+- Non-interactive mode for CI jobs (with `--no-show` flag).
 
 ## Requirements
 
@@ -35,18 +37,19 @@ showing satellite visibility over time.
 ### Basic Command
 
 ```
-./skyplot.py <RINEX_FILE> [LATITUDE] [LONGITUDE] [ALTITUDE] [--no-show]
+./skyplot.py <RINEX_FILE> [LATITUDE] [LONGITUDE] [ALTITUDE] [--use-obs] [--no-show]
 ```
 
 ### Arguments
 
-| Argument     | Type     | Units       | Description            | Default  |
-| ------------ | -------- | ----------- | ---------------------- | -------- |
-| `RINEX_FILE` | Required | -           | RINEX nav file path    | -        |
-| `LATITUDE`   | Optional | degrees (°) | North/South position   | 41.275°N |
-| `LONGITUDE`  | Optional | degrees (°) | East/West position     | 1.9876°E |
-| `ALTITUDE`   | Optional | meters (m)  | Height above sea level | 80.0 m   |
-| `--no-show`  | Optional | -           | Do not show plot       | -        |
+| Argument         | Type     | Units       | Description            | Default  |
+| ---------------- | -------- | ----------- | ---------------------- | -------- |
+| `RINEX_NAV_FILE` | Required | -           | RINEX nav file path    | -        |
+| `LATITUDE`       | Optional | degrees (°) | North/South position   | 41.275°N |
+| `LONGITUDE`      | Optional | degrees (°) | East/West position     | 1.9876°E |
+| `ALTITUDE`       | Optional | meters (m)  | Height above sea level | 80.0 m   |
+| `--use-obs`      | Optional | -           | Use RINEX obs data     | -        |
+| `--no-show`      | Optional | -           | Do not show plot       | -        |
 
 ### Examples
 
@@ -58,9 +61,9 @@ showing satellite visibility over time.
   ```
   ./skyplot.py brdc0010.22n 40.7128 -74.0060 10.0
   ```
-- Skyplot from custom location (Santiago, Chile):
+- Skyplot from custom location (Santiago, Chile) using receiver active window:
   ```
-  ./skyplot.py brdc0010.22n -33.4592 -70.6453 520.0
+  ./skyplot.py brdc0010.22n -33.4592 -70.6453 520.0 --use-obs
   ```
 - Non-interactive mode (for CI jobs):
   ```
@@ -73,6 +76,8 @@ The script generates a PDF file named `skyplot_<RINEX_FILE>.pdf` (with dots in
 `<RINEX_FILE>` replaced by `_`) with:
 
 - Satellite trajectories over all epochs in the file.
+  - Entire NAV file span (default)
+  - Receiver active span if `--use-obs` is specified and OBS file is found
 - Color-coded by constellation.
 - Observer location in title.
 - Time range in footer.

utils/skyplot/skyplot.py

@@ -1,20 +1,21 @@
 #!/usr/bin/env python
 """
- skyplot.py
+skyplot.py
 
- Reads a RINEX navigation file and generates a skyplot
+Reads a RINEX navigation file and generates a skyplot. Optionally, a RINEX observation file can
+also be read to match the skyplot to the receiver processing time.
 
- Usage: python skyplot.py <RINEX_NAV_FILE> [observer_lat] [observer_lon] [observer_alt]
+Usage: python skyplot.py <RINEX_NAV_FILE> [observer_lat] [observer_lon] [observer_alt] [--use-obs]
 
- -----------------------------------------------------------------------------
+-----------------------------------------------------------------------------
 
- GNSS-SDR is a Global Navigation Satellite System software-defined receiver.
- This file is part of GNSS-SDR.
+GNSS-SDR is a Global Navigation Satellite System software-defined receiver.
+This file is part of GNSS-SDR.
 
- SPDX-FileCopyrightText: 2025 Carles Fernandez-Prades cfernandez(at)cttc.es
- SPDX-License-Identifier: GPL-3.0-or-later
+SPDX-FileCopyrightText: 2025 Carles Fernandez-Prades cfernandez(at)cttc.es
+SPDX-License-Identifier: GPL-3.0-or-later
 
- -----------------------------------------------------------------------------
+-----------------------------------------------------------------------------
 """
 
 import argparse
@@ -28,6 +29,65 @@ import matplotlib.pyplot as plt
 import numpy as np
 
 
+def _read_obs_time_bounds(obs_path):
+    """
+    Return (start_time, end_time) from a RINEX 2/3 OBS file by scanning epoch lines.
+    If parsing fails or the file is not OBS, return (None, None).
+    """
+    try:
+        with open(obs_path, "r", encoding="utf-8", errors="ignore") as f:
+            # Detect OBS in header and skip to END OF HEADER
+            is_obs = False
+            for line in f:
+                if "RINEX VERSION / TYPE" in line:
+                    # Robust OBS detection:
+                    # - In RINEX 2/3 the file-type letter at col 21 (0-based idx 20)
+                    tchar = line[20:21]
+                    if tchar == "O" or "OBSERVATION DATA" in line:
+                        is_obs = True
+                if "END OF HEADER" in line:
+                    break
+            if not is_obs:
+                return (None, None)
+
+            start_time = None
+            end_time = None
+
+            for line in f:
+                if not line.strip():
+                    continue
+                if line.startswith(">"):  # RINEX 3 epoch line
+                    yyyy = int(line[2:6])
+                    mm = int(line[7:9])
+                    dd = int(line[10:12])
+                    hh = int(line[13:15])
+                    mi = int(line[16:18])
+                    ss = float(line[19:29])
+                    epoch = datetime(yyyy, mm, dd, hh, mi, int(ss), int((ss % 1) * 1e6))
+                else:
+                    # RINEX 2 epoch line
+                    try:
+                        yy = int(line[1:3])
+                        mm = int(line[4:6])
+                        dd = int(line[7:9])
+                        hh = int(line[10:12])
+                        mi = int(line[13:15])
+                        ss = float(line[15:26])
+                        yyyy = 1900 + yy if yy >= 80 else 2000 + yy
+                        epoch = datetime(
+                            yyyy, mm, dd, hh, mi, int(ss), int((ss % 1) * 1e6)
+                        )
+                    except Exception:
+                        continue
+                if start_time is None or epoch < start_time:
+                    start_time = epoch
+                if end_time is None or epoch > end_time:
+                    end_time = epoch
+            return (start_time, end_time)
+    except Exception:
+        return (None, None)
+
+
 def parse_rinex_float(s):
     """Parse RINEX formatted float string which may contain D or E exponent and compact spacing"""
     # Handle empty string
@@ -35,22 +95,22 @@ def parse_rinex_float(s):
         return 0.0
 
     # Replace D exponent with E (some RINEX files use D instead of E)
-    s = s.replace('D', 'E').replace('d', 'e')
+    s = s.replace("D", "E").replace("d", "e")
 
     # Handle cases where exponent lacks E (e.g., "12345-3")
-    if re.match(r'[+-]?\d+[+-]\d+', s.strip()):
-        s = s.replace('+', 'E+').replace('-', 'E-')
+    if re.match(r"[+-]?\d+[+-]\d+", s.strip()):
+        s = s.replace("+", "E+").replace("-", "E-")
 
     try:
         return float(s)
     except ValueError:
         # Handle cases where the number runs into the next field
         # Try to split at the exponent if present
-        if 'E' in s:
-            base, exp = s.split('E')[:2]
+        if "E" in s:
+            base, exp = s.split("E")[:2]
             # Take first character of exponent if needed
-            if exp and exp[0] in '+-' and len(exp) > 1:
-                return float(base + 'E' + exp[0] + exp[1:].split()[0])
+            if exp and exp[0] in "+-" and len(exp) > 1:
+                return float(base + "E" + exp[0] + exp[1:].split()[0])
         return 0.0  # Default if parsing fails
 
 
@@ -58,7 +118,7 @@ def read_rinex_nav(filename):
     """Read RINEX v3.0 navigation file"""
     satellites = {}
     line_number = 0
-    with open(filename, 'r', encoding='utf-8') as f:
+    with open(filename, "r", encoding="utf-8") as f:
         # Skip header
         while True:
             line = f.readline()
@@ -94,7 +154,7 @@ def read_rinex_nav(filename):
 
                 # Read the next lines
                 lines = [current_line]
-                line_count = 4 if system == 'R' else 7
+                line_count = 4 if system == "R" else 7
                 for _ in range(line_count):
                     next_line = f.readline()
                     line_number += 1
@@ -107,61 +167,64 @@ def read_rinex_nav(filename):
                     line_number += 1
                     continue
 
-                if system == 'R':  # GLONASS specific parsing
+                if system == "R":  # GLONASS specific parsing
                     ephemeris = {
-                        'prn': prn,
-                        'epoch': epoch,
-                        'sv_clock_bias': parse_rinex_float(lines[0][23:42]),
-                        'sv_relative_freq_bias': parse_rinex_float(lines[0][42:61]),
-                        'message_frame_time': parse_rinex_float(lines[0][61:80]),
-                        'x': parse_rinex_float(lines[1][4:23]),  # Position (km)
-                        'x_vel': parse_rinex_float(lines[1][23:42]),  # Velocity (km/s)
-                        'x_acc': parse_rinex_float(lines[1][42:61]),
-                        'health': parse_rinex_float(lines[1][61:80]),
-                        'y': parse_rinex_float(lines[2][4:23]),
-                        'y_vel': parse_rinex_float(lines[2][23:42]),
-                        'y_acc': parse_rinex_float(lines[2][42:61]),
-                        'freq_num': parse_rinex_float(lines[2][61:80]),
-                        'z': parse_rinex_float(lines[3][4:23]),
-                        'z_vel': parse_rinex_float(lines[3][23:42]),
-                        'z_acc': parse_rinex_float(lines[3][42:61]),
-                        'age': parse_rinex_float(lines[3][61:80])
+                        "prn": prn,
+                        "epoch": epoch,
+                        "sv_clock_bias": parse_rinex_float(lines[0][23:42]),
+                        "sv_relative_freq_bias": parse_rinex_float(lines[0][42:61]),
+                        "message_frame_time": parse_rinex_float(lines[0][61:80]),
+                        "x": parse_rinex_float(lines[1][4:23]),  # Position (km)
+                        "x_vel": parse_rinex_float(lines[1][23:42]),  # Velocity (km/s)
+                        "x_acc": parse_rinex_float(lines[1][42:61]),
+                        "health": parse_rinex_float(lines[1][61:80]),
+                        "y": parse_rinex_float(lines[2][4:23]),
+                        "y_vel": parse_rinex_float(lines[2][23:42]),
+                        "y_acc": parse_rinex_float(lines[2][42:61]),
+                        "freq_num": parse_rinex_float(lines[2][61:80]),
+                        "z": parse_rinex_float(lines[3][4:23]),
+                        "z_vel": parse_rinex_float(lines[3][23:42]),
+                        "z_acc": parse_rinex_float(lines[3][42:61]),
+                        "age": parse_rinex_float(lines[3][61:80]),
                     }
                 else:
                     # Parse all ephemeris parameters
                     ephemeris = {
-                        'prn': prn,
-                        'epoch': epoch,
-                        'sv_clock_bias': parse_rinex_float(lines[0][23:42]),
-                        'sv_clock_drift': parse_rinex_float(lines[0][42:61]),
-                        'sv_clock_drift_rate': parse_rinex_float(lines[0][61:80]),
-                        'iode': parse_rinex_float(lines[1][4:23]),
-                        'crs': parse_rinex_float(lines[1][23:42]),
-                        'delta_n': parse_rinex_float(lines[1][42:61]),
-                        'm0': parse_rinex_float(lines[1][61:80]),
-                        'cuc': parse_rinex_float(lines[2][4:23]),
-                        'ecc': parse_rinex_float(lines[2][23:42]),
-                        'cus': parse_rinex_float(lines[2][42:61]),
-                        'sqrt_a': parse_rinex_float(lines[2][61:80]),
-                        'toe': parse_rinex_float(lines[3][4:23]),
-                        'cic': parse_rinex_float(lines[3][23:42]),
-                        'omega0': parse_rinex_float(lines[3][42:61]),
-                        'cis': parse_rinex_float(lines[3][61:80]),
-                        'i0': parse_rinex_float(lines[4][4:23]),
-                        'crc': parse_rinex_float(lines[4][23:42]),
-                        'omega': parse_rinex_float(lines[4][42:61]),
-                        'omega_dot': parse_rinex_float(lines[4][61:80]),
-                        'idot': parse_rinex_float(lines[5][4:23]),
-                        'codes_l2': parse_rinex_float(lines[5][23:42]),
-                        'gps_week': parse_rinex_float(lines[5][42:61]),
-                        'l2p_flag': parse_rinex_float(lines[5][61:80]),
-                        'sv_accuracy': parse_rinex_float(lines[6][4:23]),
-                        'sv_health': parse_rinex_float(lines[6][23:42]),
-                        'tgd': parse_rinex_float(lines[6][42:61]),
-                        'iodc': parse_rinex_float(lines[6][61:80]),
-                        'transmission_time': parse_rinex_float(lines[7][4:23]),
-                        'fit_interval': (
-                            parse_rinex_float(lines[7][23:42])) if len(lines[7]) > 23 else 0.0
+                        "prn": prn,
+                        "epoch": epoch,
+                        "sv_clock_bias": parse_rinex_float(lines[0][23:42]),
+                        "sv_clock_drift": parse_rinex_float(lines[0][42:61]),
+                        "sv_clock_drift_rate": parse_rinex_float(lines[0][61:80]),
+                        "iode": parse_rinex_float(lines[1][4:23]),
+                        "crs": parse_rinex_float(lines[1][23:42]),
+                        "delta_n": parse_rinex_float(lines[1][42:61]),
+                        "m0": parse_rinex_float(lines[1][61:80]),
+                        "cuc": parse_rinex_float(lines[2][4:23]),
+                        "ecc": parse_rinex_float(lines[2][23:42]),
+                        "cus": parse_rinex_float(lines[2][42:61]),
+                        "sqrt_a": parse_rinex_float(lines[2][61:80]),
+                        "toe": parse_rinex_float(lines[3][4:23]),
+                        "cic": parse_rinex_float(lines[3][23:42]),
+                        "omega0": parse_rinex_float(lines[3][42:61]),
+                        "cis": parse_rinex_float(lines[3][61:80]),
+                        "i0": parse_rinex_float(lines[4][4:23]),
+                        "crc": parse_rinex_float(lines[4][23:42]),
+                        "omega": parse_rinex_float(lines[4][42:61]),
+                        "omega_dot": parse_rinex_float(lines[4][61:80]),
+                        "idot": parse_rinex_float(lines[5][4:23]),
+                        "codes_l2": parse_rinex_float(lines[5][23:42]),
+                        "gps_week": parse_rinex_float(lines[5][42:61]),
+                        "l2p_flag": parse_rinex_float(lines[5][61:80]),
+                        "sv_accuracy": parse_rinex_float(lines[6][4:23]),
+                        "sv_health": parse_rinex_float(lines[6][23:42]),
+                        "tgd": parse_rinex_float(lines[6][42:61]),
+                        "iodc": parse_rinex_float(lines[6][61:80]),
+                        "transmission_time": parse_rinex_float(lines[7][4:23]),
+                        "fit_interval": (
+                            (parse_rinex_float(lines[7][23:42]))
+                            if len(lines[7]) > 23
+                            else 0.0
+                        ),
                     }
 
                 if prn not in satellites:
@@ -189,52 +252,60 @@ def read_rinex_nav(filename):
 
 def calculate_satellite_position(ephemeris, transmit_time):
     """Calculate satellite position in ECEF coordinates at given transmission time"""
-    system = ephemeris['prn'][0]
+    system = ephemeris["prn"][0]
 
-    if system == 'R':  # GLONASS - use position + velocity * time
+    if system == "R":  # GLONASS - use position + velocity * time
         dt = transmit_time
         # Convert km to meters
-        xk = (ephemeris['x'] + ephemeris['x_vel'] * dt + 0.5 * ephemeris['x_acc'] * dt**2) * 1000
-        yk = (ephemeris['y'] + ephemeris['y_vel'] * dt + 0.5 * ephemeris['y_acc'] * dt**2) * 1000
-        zk = (ephemeris['z'] + ephemeris['z_vel'] * dt + 0.5 * ephemeris['z_acc'] * dt**2) * 1000
+        xk = (
+            ephemeris["x"] + ephemeris["x_vel"] * dt + 0.5 * ephemeris["x_acc"] * dt**2
+        ) * 1000
+        yk = (
+            ephemeris["y"] + ephemeris["y_vel"] * dt + 0.5 * ephemeris["y_acc"] * dt**2
+        ) * 1000
+        zk = (
+            ephemeris["z"] + ephemeris["z_vel"] * dt + 0.5 * ephemeris["z_acc"] * dt**2
+        ) * 1000
     else:
         # Constants
         mu = 3.986005e14  # Earth's gravitational constant (m^3/s^2)
         omega_e_dot = 7.2921151467e-5  # Earth rotation rate (rad/s)
 
         # Semi-major axis
-        a = ephemeris['sqrt_a'] ** 2
+        a = ephemeris["sqrt_a"] ** 2
 
         # Corrected mean motion
-        n0 = sqrt(mu / (a ** 3))
-        n = n0 + ephemeris['delta_n']
+        n0 = sqrt(mu / (a**3))
+        n = n0 + ephemeris["delta_n"]
 
         # Mean anomaly
-        mk = ephemeris['m0'] + n * transmit_time
+        mk = ephemeris["m0"] + n * transmit_time
 
         # Solve Kepler's equation for eccentric anomaly (Ek)
         ek = mk
         for _ in range(10):
             ek_old = ek
-            ek = mk + ephemeris['ecc'] * sin(ek)
+            ek = mk + ephemeris["ecc"] * sin(ek)
             if abs(ek - ek_old) < 1e-12:
                 break
 
         # True anomaly
-        nu_k = atan2(sqrt(1 - ephemeris['ecc']**2) * sin(ek), cos(ek) - ephemeris['ecc'])
+        nu_k = atan2(
+            sqrt(1 - ephemeris["ecc"] ** 2) * sin(ek), cos(ek) - ephemeris["ecc"]
+        )
 
         # Argument of latitude
-        phi_k = nu_k + ephemeris['omega']
+        phi_k = nu_k + ephemeris["omega"]
 
         # Second harmonic perturbations
-        delta_uk = ephemeris['cus'] * sin(2 * phi_k) + ephemeris['cuc'] * cos(2 * phi_k)
-        delta_rk = ephemeris['crs'] * sin(2 * phi_k) + ephemeris['crc'] * cos(2 * phi_k)
-        delta_ik = ephemeris['cis'] * sin(2 * phi_k) + ephemeris['cic'] * cos(2 * phi_k)
+        delta_uk = ephemeris["cus"] * sin(2 * phi_k) + ephemeris["cuc"] * cos(2 * phi_k)
+        delta_rk = ephemeris["crs"] * sin(2 * phi_k) + ephemeris["crc"] * cos(2 * phi_k)
+        delta_ik = ephemeris["cis"] * sin(2 * phi_k) + ephemeris["cic"] * cos(2 * phi_k)
 
         # Corrected argument of latitude, radius and inclination
         uk = phi_k + delta_uk
-        rk = a * (1 - ephemeris['ecc'] * cos(ek)) + delta_rk
-        ik = ephemeris['i0'] + delta_ik + ephemeris['idot'] * transmit_time
+        rk = a * (1 - ephemeris["ecc"] * cos(ek)) + delta_rk
+        ik = ephemeris["i0"] + delta_ik + ephemeris["idot"] * transmit_time
 
         # Positions in orbital plane
         xk_prime = rk * cos(uk)
@@ -242,9 +313,9 @@ def calculate_satellite_position(ephemeris, transmit_time):
 
         # Corrected longitude of ascending node
         omega_k = (
-            ephemeris['omega0']
-            + (ephemeris['omega_dot'] - omega_e_dot) * transmit_time
-            - omega_e_dot * ephemeris['toe']
+            ephemeris["omega0"]
+            + (ephemeris["omega_dot"] - omega_e_dot) * transmit_time
+            - omega_e_dot * ephemeris["toe"]
         )
 
         # Earth-fixed coordinates
@@ -261,10 +332,10 @@ def calculate_satellite_positions(ephemeris, start_time, end_time, step_min=5):
     """
     positions = []
     current_time = start_time
-    system = ephemeris['prn'][0]
-    max_valid_time = 1800 if system == 'R' else 14400
+    system = ephemeris["prn"][0]
+    max_valid_time = 1800 if system == "R" else 14400
     while current_time <= end_time:
-        transmit_time = (current_time - ephemeris['epoch']).total_seconds()
+        transmit_time = (current_time - ephemeris["epoch"]).total_seconds()
 
         if abs(transmit_time) <= max_valid_time:
             x, y, z = calculate_satellite_position(ephemeris, transmit_time)
@@ -282,7 +353,7 @@ def ecef_to_az_el(x, y, z, obs_lat, obs_lon, obs_alt):
     e_sq = 6.69437999014e-3  # first eccentricity squared
 
     # Convert geodetic coordinates to ECEF
-    n = a / sqrt(1 - e_sq * sin(obs_lat)**2)
+    n = a / sqrt(1 - e_sq * sin(obs_lat) ** 2)
     obs_x = (n + obs_alt) * cos(obs_lat) * cos(obs_lon)
     obs_y = (n + obs_alt) * cos(obs_lat) * sin(obs_lon)
     obs_z = (n * (1 - e_sq) + obs_alt) * sin(obs_lat)
@@ -294,8 +365,16 @@ def ecef_to_az_el(x, y, z, obs_lat, obs_lon, obs_alt):
 
     # Convert to local ENU (East, North, Up) coordinates
     enu_x = -sin(obs_lon) * dx + cos(obs_lon) * dy
-    enu_y = -sin(obs_lat) * cos(obs_lon) * dx - sin(obs_lat) * sin(obs_lon) * dy + cos(obs_lat) * dz
-    enu_z = cos(obs_lat) * cos(obs_lon) * dx + cos(obs_lat) * sin(obs_lon) * dy + sin(obs_lat) * dz
+    enu_y = (
+        -sin(obs_lat) * cos(obs_lon) * dx
+        - sin(obs_lat) * sin(obs_lon) * dy
+        + cos(obs_lat) * dz
+    )
+    enu_z = (
+        cos(obs_lat) * cos(obs_lon) * dx
+        + cos(obs_lat) * sin(obs_lon) * dy
+        + sin(obs_lat) * dz
+    )
 
     # Calculate azimuth and elevation
     azimuth = atan2(enu_x, enu_y)
@@ -308,44 +387,52 @@ def ecef_to_az_el(x, y, z, obs_lat, obs_lon, obs_alt):
     return azimuth, elevation
 
 
-def plot_satellite_tracks(satellites, obs_lat, obs_lon, obs_alt,
-                         footer_text=None, filename=None,
-                         show_plot=True):
+def plot_satellite_tracks(
+    satellites,
+    obs_lat,
+    obs_lon,
+    obs_alt,
+    footer_text=None,
+    filename=None,
+    show_plot=True,
+    start_time=None,
+    end_time=None,
+):
     """Plot trajectories for all visible satellites"""
     plt.rcParams["font.family"] = "Times New Roman"
     fig = plt.figure(figsize=(8, 8))
-    ax = fig.add_subplot(111, projection='polar')
+    ax = fig.add_subplot(111, projection="polar")
     ax.tick_params(labelsize=16, pad=7)
 
     # Polar plot setup
-    ax.set_theta_zero_location('N')
+    ax.set_theta_zero_location("N")
     ax.set_theta_direction(-1)
     ax.set_ylim(0, 90)
 
     # Elevation ticks
     ax.set_yticks(range(0, 91, 15))
-    ax.set_yticklabels(['90°', '', '60°', '', '30°', '', '0°'], fontsize=14)
+    ax.set_yticklabels(["90°", "", "60°", "", "30°", "", "0°"], fontsize=14)
 
     # Color scheme by constellation
     system_colors = {
-        'G': 'blue',    # GPS
-        'E': 'green',   # Galileo
-        'R': 'red',     # GLONASS
-        'C': 'orange',  # BeiDou
-        'J': 'brown',   # QZSS
-        'I': 'pink',    # IRNSS
-        'S': 'gray'     # SBAS
+        "G": "blue",  # GPS
+        "E": "green",  # Galileo
+        "R": "red",  # GLONASS
+        "C": "orange",  # BeiDou
+        "J": "brown",  # QZSS
+        "I": "pink",  # IRNSS
+        "S": "gray",  # SBAS
     }
 
     # System names mapping
     system_names = {
-        'G': 'GPS',
-        'E': 'Galileo',
-        'R': 'GLONASS',
-        'C': 'BeiDou',
-        'J': 'QZSS',
-        'I': 'IRNSS',
-        'S': 'SBAS'
+        "G": "GPS",
+        "E": "Galileo",
+        "R": "GLONASS",
+        "C": "BeiDou",
+        "J": "QZSS",
+        "I": "IRNSS",
+        "S": "SBAS",
     }
 
     # Find which systems are actually present
@@ -353,17 +440,35 @@ def plot_satellite_tracks(satellites, obs_lat, obs_lon, obs_alt,
 
     # Plot each satellite
     for prn, ephemeris_list in satellites.items():
-        color = system_colors.get(prn[0], 'purple')  # Default to purple for unknown systems
+        color = system_colors.get(
+            prn[0], "purple"
+        )  # Default to purple for unknown systems
 
         # Get the most recent ephemeris
         if not ephemeris_list:
             continue
-        ephemeris = max(ephemeris_list, key=lambda x: x['epoch'])
+
+        mid_time = start_time + (end_time - start_time) / 2
+        prev_eph = [
+            e for e in ephemeris_list if e["epoch"] <= mid_time
+        ]  # Previous ephemeris
+
+        if prev_eph:  # Pick the most recent ephemeris before or at the midpoint
+            ephemeris = max(prev_eph, key=lambda e: e["epoch"])
+        else:  # pick the ephemeris whose epoch is closest in time to the midpoint
+            ephemeris = min(
+                ephemeris_list,
+                key=lambda e: abs((e["epoch"] - mid_time).total_seconds()),
+            )
 
         # Calculate trajectory
-        all_epochs = sorted({e['epoch'] for prn_data in satellites.values() for e in prn_data})
-        start_time = min(all_epochs)
-        end_time = max(all_epochs)
+        if start_time is None or end_time is None:
+            all_epochs = sorted(
+                {e["epoch"] for prn_data in satellites.values() for e in prn_data}
+            )
+            start_time = min(all_epochs)
+            end_time = max(all_epochs)
+
         positions = calculate_satellite_positions(ephemeris, start_time, end_time)
 
         az = []
@@ -380,7 +485,7 @@ def plot_satellite_tracks(satellites, obs_lat, obs_lon, obs_alt,
             r = 90 - np.array(el)
 
             # Plot trajectory
-            ax.plot(theta, r, '-', color=color, alpha=0.7, linewidth=2.5)
+            ax.plot(theta, r, "-", color=color, alpha=0.7, linewidth=2.5)
 
             # Add arrow at last point
             if len(theta) >= 2:  # Need at least 2 points for direction
@@ -391,30 +496,43 @@ def plot_satellite_tracks(satellites, obs_lat, obs_lon, obs_alt,
                 arrow_length_factor = 1.3
                 extended_theta = theta[-2] + dx * arrow_length_factor
                 extended_r = r[-2] + dy * arrow_length_factor
-                ax.annotate('',
-                           xytext=(theta[-1], r[-1]),
-                           xy=(extended_theta, extended_r),
-                           arrowprops={
-                            'arrowstyle': '->',
-                            'color': color,
-                            'alpha': 0.9,
-                            'linewidth': 1.5,
-                            'shrinkA': 0,
-                            'shrinkB': 0
-                          })
+                ax.annotate(
+                    "",
+                    xytext=(theta[-1], r[-1]),
+                    xy=(extended_theta, extended_r),
+                    arrowprops={
+                        "arrowstyle": "->",
+                        "color": color,
+                        "alpha": 0.9,
+                        "linewidth": 1.5,
+                        "shrinkA": 0,
+                        "shrinkB": 0,
+                    },
+                )
 
             # Label at midpoint
-            mid_idx = len(theta)//2
-            ax.text(theta[mid_idx], r[mid_idx], prn,
-                    fontsize=12, ha='center', va='center',
-                    bbox={"facecolor": "white", "alpha": 0.8, "pad": 2})
+            mid_idx = len(theta) // 2
+            ax.text(
+                theta[mid_idx],
+                r[mid_idx],
+                prn,
+                fontsize=12,
+                ha="center",
+                va="center",
+                bbox={"facecolor": "white", "alpha": 0.8, "pad": 2},
+            )
 
     # Create legend elements only for present systems
     legend_elements = [
-        plt.Line2D([0], [0], marker='o', color='w',
-                  label=f'{system_names[sys]} ({sys})',
-                  markerfacecolor=system_colors[sys],
-                  markersize=10)
+        plt.Line2D(
+            [0],
+            [0],
+            marker="o",
+            color="w",
+            label=f"{system_names[sys]} ({sys})",
+            markerfacecolor=system_colors[sys],
+            markersize=10,
+        )
         for sys in present_systems
     ]
 
@@ -422,36 +540,36 @@ def plot_satellite_tracks(satellites, obs_lat, obs_lon, obs_alt,
     if legend_elements:
         ax.legend(
             handles=legend_elements,
-            loc='upper right',
+            loc="upper right",
             bbox_to_anchor=(1.3, 1.1),
-            fontsize=14
+            fontsize=14,
         )
 
     lat_deg = np.degrees(obs_lat)
     lon_deg = np.degrees(obs_lon)
-    lat_hemisphere = 'N' if lat_deg >= 0 else 'S'
-    lon_hemisphere = 'E' if lon_deg >= 0 else 'W'
+    lat_hemisphere = "N" if lat_deg >= 0 else "S"
+    lon_hemisphere = "E" if lon_deg >= 0 else "W"
 
     plt.title(
         f"GNSS skyplot from {abs(lat_deg):.2f}° {lat_hemisphere}, "
         f"{abs(lon_deg):.2f}° {lon_hemisphere}",
         pad=25,
-        fontsize=20
+        fontsize=20,
     )
 
     if footer_text:
-        fig.text(0.42, 0.05, footer_text, ha='center', va='center', fontsize=16)
+        fig.text(0.42, 0.05, footer_text, ha="center", va="center", fontsize=16)
 
     plt.tight_layout()
 
     if filename:
         filename_no_path = Path(filename).name
-        filename_no_dots = filename_no_path.replace('.', '_')
+        filename_no_dots = filename_no_path.replace(".", "_")
         output_name = f"skyplot_{filename_no_dots}.pdf"
     else:
         output_name = "skyplot.pdf"
 
-    plt.savefig(output_name, format='pdf', bbox_inches='tight')
+    plt.savefig(output_name, format="pdf", bbox_inches="tight")
     print(f"Image saved as {output_name}")
     if show_plot:
         plt.show()
@@ -463,28 +581,33 @@ def main():
     """Generate the skyplot"""
     # Set up argument parser
     parser = argparse.ArgumentParser(
-        description='Generate GNSS skyplot from RINEX navigation file',
-        add_help=False
+        description="Generate GNSS skyplot from RINEX navigation file", add_help=False
     )
 
-    # Add only the no-show flag
+    # Add the no-show flag
     parser.add_argument(
-        '--no-show', 
-        action='store_true',
-        help='Run without displaying plot window'
+        "--no-show", action="store_true", help="Run without displaying plot window"
+    )
+    # Add the observation filename
+    parser.add_argument(
+        "--use-obs",
+        action="store_true",
+        help="Use corresponding RINEX observation file to bound the skyplot to the receiver time window",
     )
 
     # Parse known args (this ignores other positional args)
     args, remaining_args = parser.parse_known_args()
 
     # Handle help manually
-    if '-h' in remaining_args or '--help' in remaining_args:
-        print("""
-Usage: python skyplot.py <RINEX_FILE> [LATITUDE] [LONGITUDE] [ALTITUDE] [--no-show]
+    if "-h" in remaining_args or "--help" in remaining_args:
+        print(
+            """
+Usage: python skyplot.py <RINEX_FILE> [LATITUDE] [LONGITUDE] [ALTITUDE] [--use-obs] [--no-show]
 
 Example:
-  python skyplot.py brdc0010.22n 41.275 1.9876 80.0 --no-show
-""")
+  python skyplot.py brdc0010.22n 41.275 1.9876 80.0 --use-obs --no-show
+"""
+        )
         sys.exit(0)
 
     if len(remaining_args) < 1:
@@ -521,9 +644,15 @@ Example:
         return
 
     # Print summary information
-    all_epochs = sorted(list(set(
-        e['epoch'] for prn, ephemerides in satellites.items() for e in ephemerides
-    )))
+    all_epochs = sorted(
+        list(
+            set(
+                e["epoch"]
+                for prn, ephemerides in satellites.items()
+                for e in ephemerides
+            )
+        )
+    )
     print("\nFile contains:")
     print(f"- {len(satellites)} unique satellites")
     print(f"- {len(all_epochs)} unique epochs")
@@ -537,17 +666,50 @@ Example:
 
     print("\nSatellite systems found:")
     for system, count in sorted(system_counts.items()):
-        system_name = {
-            'G': 'GPS',
-            'R': 'GLONASS',
-            'E': 'Galileo',
-            'C': 'BeiDou'
-        }.get(system, 'Unknown')
+        system_name = {"G": "GPS", "R": "GLONASS", "E": "Galileo", "C": "BeiDou"}.get(
+            system, "Unknown"
+        )
         print(f"- {system_name} ({system}): {count} satellites")
 
+    # Generate the combined skyplot
+    # Time window: OBS bounds if provided; else NAV span
+    use_start, use_end = all_epochs[0], all_epochs[-1]
+    if args.use_obs:
+        tried = []
+        obs_path = None
+        stem = filename[:-1]
+
+        for s in ("O", "o"):  # try uppercase then lowercase
+            candidate = stem + s
+            tried.append(candidate)
+            if Path(candidate).exists():
+                obs_path = candidate
+                break
+
+        if obs_path:
+            obs_start, obs_end = _read_obs_time_bounds(obs_path)
+            if obs_start and obs_end:
+                use_start, use_end = obs_start, obs_end
+                print(
+                    f"\nObservation window detected from {obs_path}: {use_start} → {use_end}"
+                )
+            else:
+                print(
+                    f"\nWarning: Could not read valid times from {obs_path}. Using NAV span instead."
+                )
+        else:
+            print(
+                f"\nOBS file not found. Tried: {', '.join(tried)}. Using NAV span instead."
+            )
+
+    # Ensure at least two samples with the default 5-minute step
+    if (use_end - use_start) < timedelta(minutes=5):
+        use_end = use_start + timedelta(minutes=5)
+
     # Generate the combined skyplot
     print("\nGenerating skyplot...")
-    footer = f"From {all_epochs[0]} to {all_epochs[-1]} UTC"
+    footer = f"From {use_start} to {use_end} UTC"
+
     plot_satellite_tracks(
         satellites,
         obs_lat,
@@ -555,7 +717,9 @@ Example:
         obs_alt,
         footer_text=footer,
         filename=filename,
-        show_plot=not args.no_show
+        show_plot=not args.no_show,
+        start_time=use_start,
+        end_time=use_end,
     )