Calculating Doppler Compensation¶

In this example we'll compute the Doppler compensation required for a ground station to maintain communications with a satellite in low Earth orbit (LEO). We'll simulate the International Space Station (ISS) passing over the Cape Canaveral ground station from the NASA Near Earth Network (NEN) and calculate the frequency shift due to the relative motion between the satellite and ground station during the access window.

To accomplish this we'll define a custom access property computer that calculates the Doppler shift in Hz for S-band uplink communications (2.2 GHz) and X-Band downlink communications (8.4 GHz) based on the relative velocity between the satellite and ground station along the line of sight. We'll assume that the satellite uses a fixed carrier frequency for both uplink and downlink, and the ground station must adjust its local oscillator frequency to compensate for the Doppler shift in order to maintain a stable communication link.

Setup¶

First, we'll import the necessary libraries, initialize Earth orientation parameters, download the TLE for the ISS (25544) from CelesTrak, and load the NASA Near Earth Network ground station network, and select the Cape Canaveral ground station.

import time
import csv
import os
import pathlib
import sys
import brahe as bh
import numpy as np
import plotly.graph_objects as go
from plotly.subplots import make_subplots

bh.initialize_eop()

We download the ISS TLE directly by NORAD ID and load all NASA NEN ground stations:

iss = bh.datasets.celestrak.get_tle_by_id_as_propagator(25544, 60.0)

Then select Cape Canaveral from the loaded stations:

nen_stations = bh.datasets.groundstations.load("nasa nen")

# Select Cape Canaveral ground station
cape_canaveral = None
for station in nen_stations:
    if "Merrit Island" in station.get_name():
        cape_canaveral = station
        break

Custom Doppler Shift Property Computer¶

Next, we'll define a custom access property computer that calculates the Doppler shift compensation required for maintaining stable communication links.

Doppler Shift Physics¶

For electromagnetic signals traveling between a moving satellite and a stationary ground station, the observed frequency differs from the transmitted frequency due to relative motion along the line-of-sight from the transmitter to the receiver. Using the non-relativistic Doppler approximation (valid when \(|v_{los}| \ll c\)), the received frequency \(f_r\) given a transmitted frequency \(f_x\) is:

\[ f_r = f_x \left( 1 - \frac{v_{los}}{c} \right) \]

where \(f_x\) is the transmitted carrier frequency, \(v_{los}\) is the velocity along the line of sight from the station to the spacecraft (negative when approaching, positive when receding), and \(c\) is the speed of light.

Sign Convention:

\(v_{los} < 0\): satellite approaching → frequency increases and wavelength decreases (blueshift)
\(v_{los} > 0\): satellite receeding → frequency decreases and wavelength increases (redshift)

Compensation Strategy¶

To maintain a stable communication link, the ground station must compensate for the Doppler shift that arises from the relative line-of-sight velocity \(v_{los}\) between the satellite and the station. Both the satellite and ground station have design frequencies that they expect to operate at for uplink. The challenge is determining how much the ground station must adjust its transmit and receive frequencies to account for the Doppler effect.

Downlink Compensation (Satellite → Ground)¶

The spacecraft transmits at a fixed design frequency \(f_x^d\) (e.g., 8.4 GHz for X-band). The ground station's initial receiver frequency \(f_r^0\) is initially tuned to \(f_r^0 = f_x^d\), expecting to receive at the design frequency. The challenge then is to determine the required frequency adjustment \(\Delta f_r = f_r - f_r^0\) that is added to base the receiver frequency to correctly receive the Doppler-shifted signal.

Starting from the Doppler relation:

\[ f_r = f_x \left( 1 - \frac{v_{los}}{c} \right) \]

where \(f_r\) is the actual received frequency, \(f_x\) is the transmitted frequency, and \(v_{los}\) is the line-of-sight velocity (negative when approaching, positive when receding).

Since the spacecraft transmits at the fixed design frequency \(f_x = f_x^d\):

\[ f_r = f_x^d \left( 1 - \frac{v_{los}}{c} \right) \]

The receiver compensation term is the difference between the actual received frequency and the initial tuning:

\[ \Delta f_r = f_r - f_r^0 = f_x^d \left( 1 - \frac{v_{los}}{c} \right) - f_r^0 \]

Since \(f_r^0 = f_x^d\):

\[ \begin{align} \Delta f_r &= f_x^d \left( 1 - \frac{v_{los}}{c} \right) - f_r^0 \\ &= f_x^d \left( 1 - \frac{v_{los}}{c} \right) - f_x^d \\ &= f_x^d \left( 1 - \frac{v_{los}}{c} - 1 \right) \\ &= -f_x^d \frac{v_{los}}{c} \end{align} \]

Therefore, the downlink compensation is:

\[ \boxed{\Delta f_r = -f_x^d \frac{v_{los}}{c}} \]

Sign interpretation: - When approaching (\(v_{los} < 0\)): \(\Delta f_r > 0\) → tune receiver higher (blueshift) - When receding (\(v_{los} > 0\)): \(\Delta f_r < 0\) → tune receiver lower (redshift)

Uplink Compensation (Ground → Satellite)¶

For uplink, the spacecraft expects to receive at its design frequency \(f_r^d\) (e.g., 2.2 GHz for S-band). The ground station initially transmits at \(f_x^0 = f_r^d\), but must pre-compensate so that after Doppler shift, the spacecraft receives exactly \(f_r^d\). We must determine the required frequency adjustment \(\Delta f_x = f_x - f_x^0\) that the ground station must apply to its transmit frequency.

We again start from the Doppler relation:

\[ \begin{align} f_r^d &= (f_x + \Delta f_x) \left( 1 - \frac{v_{los}}{c} \right) &= f_x^0 \left( 1 - \frac{v_{los}}{c} \right) + \Delta f_x \left( 1 - \frac{v_{los}}{c} \right) \end{align} \]

Solving for \(\Delta f_x\):

\[ \Delta f_x \left( 1 - \frac{v_{los}}{c} \right) = f_r^d - f_x^0 \left( 1 - \frac{v_{los}}{c} \right) \]

\[ \Delta f_x = \frac{f_r^d - f_x^0 \left( 1 - \frac{v_{los}}{c} \right)}{1 - \frac{v_{los}}{c}} \]

Since we know that \(f_x^0 = f_r^d\), we can substitute and simplify:

\[ \begin{align} \Delta f_x &= \frac{f_r^d - f_r^d \left( 1 - \frac{v_{los}}{c} \right)}{1 - \frac{v_{los}}{c}} \\ &= \frac{f_r^d \left( 1 - \left( 1 - \frac{v_{los}}{c} \right) \right)}{1 - \frac{v_{los}}{c}} \\ &= f_r^d \frac{\frac{v_{los}}{c}}{1 - \frac{v_{los}}{c}} \\ &= f_x^0 \frac{\frac{v_{los}}{c}}{\frac{c - v_{los}}{c}} \\ &= f_x^0 \frac{v_{los}}{c - v_{los}} \end{align} \]

And so we find the uplink compensation is:

\[ \boxed{\Delta f_x = f_x^0 \frac{v_{los}}{c - v_{los}}} \]

Frequency Scaling

The magnitude of Doppler compensation scales with carrier frequency. X-band (8.4 GHz) experiences about \(8.4/2.2 \approx 3.8\times\) the frequency shift of S-band (2.2 GHz) for the same line-of-sight velocity.

Implementation¶

We create a custom property computer that calculates the line-of-sight velocity from the satellite state and computes Doppler compensation for both S-band uplink (2.2 GHz) and X-band downlink (8.4 GHz):

# Communication frequency bands
S_BAND_FREQ = 2.2e9  # Hz (uplink)
X_BAND_FREQ = 8.4e9  # Hz (downlink)


class DopplerComputer(bh.AccessPropertyComputer):
    """Compute Doppler shift time series for S-band and X-band communications."""

    def sampling_config(self):
        """Configure sampling at 10 Hz during access windows."""
        interval_hz = 10.0  # 10 samples per second
        interval_s = 1.0 / interval_hz
        return bh.SamplingConfig.fixed_interval(interval_s, 0.0)

    def compute(
        self, window, sample_times, sample_states_ecef, location_ecef, location_geodetic
    ):
        """Calculate Doppler compensation frequencies at each sample point.

        Args:
            window: AccessWindow with timing information
            sample_times: Sample epochs in MJD [N]
            sample_states_ecef: Satellite states [N x 6] in ECEF (m, m/s)
            location_ecef: Location position [3] in ECEF (m)
            location_geodetic: Location geodetic coords [lon, lat, alt] in (degrees, degrees, m)

        Returns:
            dict: Time series of Doppler compensation frequencies in Hz
        """
        loc_pos = np.array(location_ecef)

        # Compute Doppler at each sample point
        doppler_s_band_list = []
        doppler_x_band_list = []
        v_los_list = []

        for state in sample_states_ecef:
            # Extract satellite position and velocity
            sat_pos = state[:3]
            sat_vel = state[3:6]

            # Compute line-of-sight vector (from ground station to satellite)
            los_vec = sat_pos - loc_pos
            los_unit = los_vec / np.linalg.norm(los_vec)

            # Compute line-of-sight velocity (negative when approaching, positive when receding)
            v_los = np.dot(sat_vel, los_unit)

            # Compute Doppler compensation from first principles
            # Uplink (S-band): Δf_x = f_x^0 × v_los / (c - v_los)
            #   Ground pre-compensates transmit frequency so spacecraft receives design frequency
            # Downlink (X-band): Δf_r = -f_x^d × v_los / c
            #   Ground adjusts receive frequency to match Doppler-shifted spacecraft transmission
            doppler_s = S_BAND_FREQ * v_los / (bh.C_LIGHT - v_los)  # Uplink
            doppler_x = -X_BAND_FREQ * v_los / bh.C_LIGHT  # Downlink

            doppler_s_band_list.append(doppler_s)
            doppler_x_band_list.append(doppler_x)
            v_los_list.append(v_los)

        # Convert sample times to relative times (seconds from window open)
        relative_times = (sample_times - window.window_open.mjd()) * 86400.0

        # Return as time series
        return {
            "doppler_s_band": {
                "times": relative_times.tolist(),
                "values": doppler_s_band_list,
            },
            "doppler_x_band": {
                "times": relative_times.tolist(),
                "values": doppler_x_band_list,
            },
            "v_los": {
                "times": relative_times.tolist(),
                "values": v_los_list,
            },
        }

    def property_names(self):
        """List properties this computer provides."""
        return ["doppler_s_band", "doppler_x_band", "v_los"]

The property computer extracts the satellite velocity from the provided ECEF state, projects it onto the line-of-sight unit vector to get the line-of-sight velocity, and applies the Doppler formula for both frequency bands.

Access Computation¶

Next, we'll compute the access windows between the ISS and the Cape Canaveral ground station over a 72-hour period with our custom Doppler shift property computer to calculate the Doppler compensation required to establish S-band and X-band communications during each access.

epoch_start = iss.epoch
epoch_end = epoch_start + 72 * 3600.0  # 72 hours in seconds

# Propagate for full 72-hour period
iss.propagate_to(epoch_end)

# Compute access windows with 5 degree minimum elevation and Doppler properties
constraint = bh.ElevationConstraint(min_elevation_deg=5.0)
doppler_computer = DopplerComputer()

windows = bh.location_accesses(
    [cape_canaveral],
    [iss],
    epoch_start,
    epoch_end,
    constraint,
    property_computers=[doppler_computer],
)

This computes all access windows where the ISS rises above 5° elevation as viewed from Cape Canaveral, and automatically calculates the Doppler compensation frequencies at the midpoint of each window using our custom property computer.

Ground Track Visualization¶

We first visualize the ISS ground track over one orbital period, showing Cape Canaveral's location and its communication cone based on the 5° minimum elevation constraint:

orbital_period = bh.orbital_period(iss.semi_major_axis)
iss.propagate_to(iss.epoch + orbital_period)


# Create ground track visualization
fig_groundtrack = bh.plot_groundtrack(
    trajectories=[
        {
            "trajectory": iss.trajectory,
            "color": "red",
            "line_width": 2.0,
            "track_length": 3,
            "track_units": "orbits",
        }
    ],
    ground_stations=[
        {
            "stations": [cape_canaveral],
            "color": "blue",
            "alpha": 0.15,
            "point_size": 8.0,
        }
    ],
    gs_cone_altitude=iss.semi_major_axis - bh.R_EARTH,
    gs_min_elevation=5.0,
    show_borders=True,
    show_coastlines=True,
    show_legend=False,
    backend="plotly",
)

The resulting plot shows the ISS ground track in red and Cape Canaveral with its communication cone in blue:

Doppler Compensation Analysis¶

For detailed analysis, we select one access window extract the computed doppler compensation value. This provides a high-resolution profile showing how the Doppler compensation varies as the satellite approaches, reaches closest approach, and recedes:

The top panel shows \(v_{los}\) (line-of-sight velocity from ground station to satellite), confirming our sign convention: negative when approaching (start of pass), crossing through zero at closest approach, then positive when receding (end of pass).

The middle and bottom panels show the Doppler compensation frequencies with opposite signs as required:

S-band uplink (middle): Uses \(\Delta f_x = f_x^0 v_{los}/(c - v_{los})\), so it's negative when approaching (transmit lower) and positive when receding (transmit higher)
X-band downlink (bottom): Uses \(\Delta f_r = -f_x^d v_{los}/c\), so it's positive when approaching (receive higher) and negative when receding (receive lower)

Note that X-band requires approximately 3.8× more compensation than S-band due to its higher carrier frequency (8.4 GHz vs 2.2 GHz). The uplink compensation includes the \((c - v_{los})\) denominator term for proper pre-compensation.

# Convert times to minutes relative to window start for better readability
times_rel_min = [t / 60.0 for t in doppler_s_ts["times"]]
doppler_s_band_khz = [d / 1000.0 for d in doppler_s_ts["values"]]
doppler_x_band_khz = [d / 1000.0 for d in doppler_x_ts["values"]]
v_los_km_s = [v / 1000.0 for v in v_los_ts["values"]]  # Convert m/s to km/s

# Create figure with three subplots
fig_doppler = make_subplots(
    rows=3,
    cols=1,
    subplot_titles=(
        "Line-of-Sight Velocity (negative = approaching, positive = receding)",
        "S-band Uplink Doppler Compensation (2.2 GHz)",
        "X-band Downlink Doppler Compensation (8.4 GHz)",
    ),
    vertical_spacing=0.12,
)

# v_los plot
fig_doppler.add_trace(
    go.Scatter(
        x=times_rel_min,
        y=v_los_km_s,
        mode="lines",
        name="v_los",
        line=dict(color="green", width=2),
        hovertemplate="Time: %{x:.2f} min<br>v_los: %{y:.2f} km/s<extra></extra>",
    ),
    row=1,
    col=1,
)

# S-band plot
fig_doppler.add_trace(
    go.Scatter(
        x=times_rel_min,
        y=doppler_s_band_khz,
        mode="lines",
        name="S-band (2.2 GHz)",
        line=dict(color="steelblue", width=2),
        hovertemplate="Time: %{x:.2f} min<br>Doppler: %{y:.2f} kHz<extra></extra>",
    ),
    row=2,
    col=1,
)

# X-band plot
fig_doppler.add_trace(
    go.Scatter(
        x=times_rel_min,
        y=doppler_x_band_khz,
        mode="lines",
        name="X-band (8.4 GHz)",
        line=dict(color="coral", width=2),
        hovertemplate="Time: %{x:.2f} min<br>Doppler: %{y:.2f} kHz<extra></extra>",
    ),
    row=3,
    col=1,
)

# Add zero reference lines
fig_doppler.add_hline(
    y=0, line_dash="dash", line_color="gray", opacity=0.5, row=1, col=1
)
fig_doppler.add_hline(
    y=0, line_dash="dash", line_color="gray", opacity=0.5, row=2, col=1
)
fig_doppler.add_hline(
    y=0, line_dash="dash", line_color="gray", opacity=0.5, row=3, col=1
)

# Update axes labels with smaller fonts (matching plot_cartesian_trajectory pattern)
axis_config = {
    "title_font": {"size": 11},
    "tickfont": {"size": 10},
}
fig_doppler.update_xaxes(
    title_text="Time from AOS (minutes)", row=3, col=1, **axis_config
)
fig_doppler.update_yaxes(title_text="v_los (km/s)", row=1, col=1, **axis_config)
fig_doppler.update_yaxes(title_text="Doppler (kHz)", row=2, col=1, **axis_config)
fig_doppler.update_yaxes(title_text="Doppler (kHz)", row=3, col=1, **axis_config)

# Make subplot titles smaller
for annotation in fig_doppler.layout.annotations:
    annotation.font.size = 11

# Update layout - NO explicit width/height for responsive sizing
fig_doppler.update_layout(
    title="Doppler Compensation for ISS Pass over Cape Canaveral",
    showlegend=False,
)

Doppler Compensation Data¶

Below is a table of sampled Doppler compensation values during the access window. The compensation frequency indicates the adjustment needed for the ground station equipment:

Time (UTC)	S-band Doppler (kHz)	X-band Doppler (kHz)
2026-01-03 10:12:11	-43.57	166.36
2026-01-03 10:12:56	-40.66	155.27
2026-01-03 10:13:41	-35.89	137.05
2026-01-03 10:14:27	-28.09	107.25
2026-01-03 10:15:12	-16.05	61.29
2026-01-03 10:15:57	-0.21	0.8
2026-01-03 10:16:43	15.7	-59.95
2026-01-03 10:17:28	27.88	-106.44
2026-01-03 10:18:13	35.81	-136.71
2026-01-03 10:18:59	40.66	-155.25
2026-01-03 10:19:44	43.63	-166.6

Full Code Example¶

doppler_compensation.py
import time
import csv
import os
import pathlib
import sys
import brahe as bh
import numpy as np
import plotly.graph_objects as go
from plotly.subplots import make_subplots

bh.initialize_eop()

# Configuration for output files
SCRIPT_NAME = pathlib.Path(__file__).stem
OUTDIR = pathlib.Path(os.getenv("BRAHE_FIGURE_OUTPUT_DIR", "./docs/figures/"))
os.makedirs(OUTDIR, exist_ok=True)


# Communication frequency bands
S_BAND_FREQ = 2.2e9  # Hz (uplink)
X_BAND_FREQ = 8.4e9  # Hz (downlink)


class DopplerComputer(bh.AccessPropertyComputer):
    """Compute Doppler shift time series for S-band and X-band communications."""

    def sampling_config(self):
        """Configure sampling at 10 Hz during access windows."""
        interval_hz = 10.0  # 10 samples per second
        interval_s = 1.0 / interval_hz
        return bh.SamplingConfig.fixed_interval(interval_s, 0.0)

    def compute(
        self, window, sample_times, sample_states_ecef, location_ecef, location_geodetic
    ):
        """Calculate Doppler compensation frequencies at each sample point.

        Args:
            window: AccessWindow with timing information
            sample_times: Sample epochs in MJD [N]
            sample_states_ecef: Satellite states [N x 6] in ECEF (m, m/s)
            location_ecef: Location position [3] in ECEF (m)
            location_geodetic: Location geodetic coords [lon, lat, alt] in (degrees, degrees, m)

        Returns:
            dict: Time series of Doppler compensation frequencies in Hz
        """
        loc_pos = np.array(location_ecef)

        # Compute Doppler at each sample point
        doppler_s_band_list = []
        doppler_x_band_list = []
        v_los_list = []

        for state in sample_states_ecef:
            # Extract satellite position and velocity
            sat_pos = state[:3]
            sat_vel = state[3:6]

            # Compute line-of-sight vector (from ground station to satellite)
            los_vec = sat_pos - loc_pos
            los_unit = los_vec / np.linalg.norm(los_vec)

            # Compute line-of-sight velocity (negative when approaching, positive when receding)
            v_los = np.dot(sat_vel, los_unit)

            # Compute Doppler compensation from first principles
            # Uplink (S-band): Δf_x = f_x^0 × v_los / (c - v_los)
            #   Ground pre-compensates transmit frequency so spacecraft receives design frequency
            # Downlink (X-band): Δf_r = -f_x^d × v_los / c
            #   Ground adjusts receive frequency to match Doppler-shifted spacecraft transmission
            doppler_s = S_BAND_FREQ * v_los / (bh.C_LIGHT - v_los)  # Uplink
            doppler_x = -X_BAND_FREQ * v_los / bh.C_LIGHT  # Downlink

            doppler_s_band_list.append(doppler_s)
            doppler_x_band_list.append(doppler_x)
            v_los_list.append(v_los)

        # Convert sample times to relative times (seconds from window open)
        relative_times = (sample_times - window.window_open.mjd()) * 86400.0

        # Return as time series
        return {
            "doppler_s_band": {
                "times": relative_times.tolist(),
                "values": doppler_s_band_list,
            },
            "doppler_x_band": {
                "times": relative_times.tolist(),
                "values": doppler_x_band_list,
            },
            "v_los": {
                "times": relative_times.tolist(),
                "values": v_los_list,
            },
        }

    def property_names(self):
        """List properties this computer provides."""
        return ["doppler_s_band", "doppler_x_band", "v_los"]




# Download TLE data for ISS from CelesTrak
# ISS (International Space Station)
# NORAD ID: 25544
print("Downloading ISS TLE from CelesTrak...")
start_time = time.time()
iss = bh.datasets.celestrak.get_tle_by_id_as_propagator(25544, 60.0)
elapsed = time.time() - start_time
print(f"Downloaded ISS TLE in {elapsed:.2f} seconds.")
print(f"Epoch: {iss.epoch}")
print(f"Semi-major axis: {iss.semi_major_axis / 1000:.1f} km")

# Load NASA Near Earth Network ground stations
print("\nLoading NASA Near Earth Network ground stations...")
start_time = time.time()
nen_stations = bh.datasets.groundstations.load("nasa nen")

# Select Cape Canaveral ground station
cape_canaveral = None
for station in nen_stations:
    if "Merrit Island" in station.get_name():
        cape_canaveral = station
        break

if cape_canaveral is None:
    print("Error: Cape Canaveral ground station not found.")
    sys.exit(1)

elapsed = time.time() - start_time
print(f"Loaded {len(nen_stations)} NASA NEN ground stations in {elapsed:.2f} seconds.")


# Propagate ISS for one orbit to create ground track visualization
start_time = time.time()
print("\nCreating ground track visualization...")
orbital_period = bh.orbital_period(iss.semi_major_axis)
iss.propagate_to(iss.epoch + orbital_period)


# Create ground track visualization
fig_groundtrack = bh.plot_groundtrack(
    trajectories=[
        {
            "trajectory": iss.trajectory,
            "color": "red",
            "line_width": 2.0,
            "track_length": 3,
            "track_units": "orbits",
        }
    ],
    ground_stations=[
        {
            "stations": [cape_canaveral],
            "color": "blue",
            "alpha": 0.15,
            "point_size": 8.0,
        }
    ],
    gs_cone_altitude=iss.semi_major_axis - bh.R_EARTH,
    gs_min_elevation=5.0,
    show_borders=True,
    show_coastlines=True,
    show_legend=False,
    backend="plotly",
)
elapsed = time.time() - start_time
print(f"Created ground track visualization in {elapsed:.2f} seconds.")

# Reset propagator and compute 72-hour access windows with Doppler properties
start_time = time.time()
print("\nComputing 72-hour access windows with Doppler compensation...")
iss.reset()


epoch_start = iss.epoch
epoch_end = epoch_start + 72 * 3600.0  # 72 hours in seconds

# Propagate for full 72-hour period
iss.propagate_to(epoch_end)

# Compute access windows with 5 degree minimum elevation and Doppler properties
constraint = bh.ElevationConstraint(min_elevation_deg=5.0)
doppler_computer = DopplerComputer()

windows = bh.location_accesses(
    [cape_canaveral],
    [iss],
    epoch_start,
    epoch_end,
    constraint,
    property_computers=[doppler_computer],
)

elapsed = time.time() - start_time
print(f"Computed {len(windows)} access windows in {elapsed:.2f} seconds.")

if len(windows) == 0:
    print("No access windows found. Exiting.")
    sys.exit(0)

# Print sample of access windows with Doppler properties
print("\n" + "=" * 100)
print("Sample Access Windows with Doppler Time Series (first 5)")
print("=" * 100)
print(
    f"{'Start Time':<25} {'Duration':>10} {'Max Elev':>10} {'# Samples':>10} {'Peak S-band':>13} {'Peak X-band':>13}"
)
print("-" * 100)
for i, window in enumerate(windows[:5]):
    duration_min = window.duration / 60.0
    max_elev = window.properties.elevation_max

    # Extract time series data
    doppler_s_ts = window.properties.additional["doppler_s_band"]
    doppler_x_ts = window.properties.additional["doppler_x_band"]

    n_samples = len(doppler_s_ts["values"])
    peak_s = max(abs(v) for v in doppler_s_ts["values"]) / 1000.0
    peak_x = max(abs(v) for v in doppler_x_ts["values"]) / 1000.0

    start_str = str(window.start).split(".")[0]  # Remove fractional seconds
    print(
        f"{start_str:<25} {duration_min:>8.1f} m {max_elev:>8.1f}° {n_samples:>10} {peak_s:>10.2f} kHz {peak_x:>10.2f} kHz"
    )
print("=" * 100)

# Select a window with good duration for detailed analysis
selected_window = None
for window in windows:
    if window.duration > 300:  # At least 5 minutes
        selected_window = window
        break

if selected_window is None:
    # Fallback to first window if none are long enough
    selected_window = windows[0]

print(
    f"\nAnalyzing detailed Doppler profile for window starting at {selected_window.start}"
)
print(f"Window duration: {selected_window.duration / 60:.1f} minutes")

# Extract precomputed time series from window properties
print("Extracting precomputed Doppler time series...")
start_time = time.time()

doppler_s_ts = selected_window.properties.additional["doppler_s_band"]
doppler_x_ts = selected_window.properties.additional["doppler_x_band"]
v_los_ts = selected_window.properties.additional["v_los"]

# Convert relative times to UTC epochs
times_utc = [selected_window.start + t for t in doppler_s_ts["times"]]

num_samples = len(doppler_s_ts["times"])

elapsed = time.time() - start_time
print(f"Extracted {num_samples} samples (0.1s interval) in {elapsed:.2f} seconds.")

# Export ~10 evenly-spaced samples to CSV
csv_path = OUTDIR / f"{SCRIPT_NAME}_data.csv"
num_csv_samples = 11
csv_indices = np.linspace(0, len(times_utc) - 1, num_csv_samples, dtype=int)

with open(csv_path, "w", newline="") as csvfile:
    writer = csv.writer(csvfile)
    writer.writerow(["Time (UTC)", "S-band Doppler (kHz)", "X-band Doppler (kHz)"])
    for idx in csv_indices:
        time_str = str(times_utc[idx]).split(".")[0]  # Remove fractional seconds
        writer.writerow(
            [
                time_str,
                f"{doppler_s_ts['values'][idx] / 1000.0:.2f}",
                f"{doppler_x_ts['values'][idx] / 1000.0:.2f}",
            ]
        )

print(f"✓ Exported {num_csv_samples} samples to {csv_path}")

# Create Doppler vs time visualization
print("\nCreating Doppler compensation visualization...")
start_time = time.time()


# Convert times to minutes relative to window start for better readability
times_rel_min = [t / 60.0 for t in doppler_s_ts["times"]]
doppler_s_band_khz = [d / 1000.0 for d in doppler_s_ts["values"]]
doppler_x_band_khz = [d / 1000.0 for d in doppler_x_ts["values"]]
v_los_km_s = [v / 1000.0 for v in v_los_ts["values"]]  # Convert m/s to km/s

# Create figure with three subplots
fig_doppler = make_subplots(
    rows=3,
    cols=1,
    subplot_titles=(
        "Line-of-Sight Velocity (negative = approaching, positive = receding)",
        "S-band Uplink Doppler Compensation (2.2 GHz)",
        "X-band Downlink Doppler Compensation (8.4 GHz)",
    ),
    vertical_spacing=0.12,
)

# v_los plot
fig_doppler.add_trace(
    go.Scatter(
        x=times_rel_min,
        y=v_los_km_s,
        mode="lines",
        name="v_los",
        line=dict(color="green", width=2),
        hovertemplate="Time: %{x:.2f} min<br>v_los: %{y:.2f} km/s<extra></extra>",
    ),
    row=1,
    col=1,
)

# S-band plot
fig_doppler.add_trace(
    go.Scatter(
        x=times_rel_min,
        y=doppler_s_band_khz,
        mode="lines",
        name="S-band (2.2 GHz)",
        line=dict(color="steelblue", width=2),
        hovertemplate="Time: %{x:.2f} min<br>Doppler: %{y:.2f} kHz<extra></extra>",
    ),
    row=2,
    col=1,
)

# X-band plot
fig_doppler.add_trace(
    go.Scatter(
        x=times_rel_min,
        y=doppler_x_band_khz,
        mode="lines",
        name="X-band (8.4 GHz)",
        line=dict(color="coral", width=2),
        hovertemplate="Time: %{x:.2f} min<br>Doppler: %{y:.2f} kHz<extra></extra>",
    ),
    row=3,
    col=1,
)

# Add zero reference lines
fig_doppler.add_hline(
    y=0, line_dash="dash", line_color="gray", opacity=0.5, row=1, col=1
)
fig_doppler.add_hline(
    y=0, line_dash="dash", line_color="gray", opacity=0.5, row=2, col=1
)
fig_doppler.add_hline(
    y=0, line_dash="dash", line_color="gray", opacity=0.5, row=3, col=1
)

# Update axes labels with smaller fonts (matching plot_cartesian_trajectory pattern)
axis_config = {
    "title_font": {"size": 11},
    "tickfont": {"size": 10},
}
fig_doppler.update_xaxes(
    title_text="Time from AOS (minutes)", row=3, col=1, **axis_config
)
fig_doppler.update_yaxes(title_text="v_los (km/s)", row=1, col=1, **axis_config)
fig_doppler.update_yaxes(title_text="Doppler (kHz)", row=2, col=1, **axis_config)
fig_doppler.update_yaxes(title_text="Doppler (kHz)", row=3, col=1, **axis_config)

# Make subplot titles smaller
for annotation in fig_doppler.layout.annotations:
    annotation.font.size = 11

# Update layout - NO explicit width/height for responsive sizing
fig_doppler.update_layout(
    title="Doppler Compensation for ISS Pass over Cape Canaveral",
    showlegend=False,
)

elapsed = time.time() - start_time
print(f"Created Doppler visualization in {elapsed:.2f} seconds.")

# ============================================================================
# Plot Output Section (for documentation generation)
# ============================================================================

# Add plots directory to path for importing brahe_theme
sys.path.insert(0, str(pathlib.Path(__file__).parent.parent.parent / "plots"))
from brahe_theme import save_themed_html  # noqa: E402

# Save the ground track figure as themed HTML
light_path, dark_path = save_themed_html(
    fig_groundtrack, OUTDIR / f"{SCRIPT_NAME}_groundtrack"
)
print(f"\n✓ Generated {light_path}")
print(f"✓ Generated {dark_path}")

# Save Doppler figure as themed HTML
light_path, dark_path = save_themed_html(fig_doppler, OUTDIR / f"{SCRIPT_NAME}_doppler")
print(f"✓ Generated {light_path}")
print(f"✓ Generated {dark_path}")

print("\nDoppler Compensation Analysis Complete!")
print(f"Peak S-band Doppler: {max(abs(d) for d in doppler_s_band_khz):.2f} kHz")
print(f"Peak X-band Doppler: {max(abs(d) for d in doppler_x_band_khz):.2f} kHz")