Maximum Communications Gap¶

import csv
import os
import pathlib
import sys
import time

import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import numpy as np
import plotly.graph_objects as go

import brahe as bh

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)

# Download TLE data for all active satellites as propagators and filter for Umbra
print("Downloading active satellite TLEs from CelesTrak...")
start_time = time.time()
client = bh.celestrak.CelestrakClient()
records = client.get_gp(name="UMBRA")

# Convert to propagators
umbra_props = [r.to_sgp_propagator(60.0) for r in records]
print(f"Found {len(umbra_props)} Umbra satellites")
elapsed = time.time() - start_time

# Load specific KSAT ground stations
print("\nLoading KSAT ground stations...")
start_time = time.time()
all_ksat = bh.datasets.groundstations.load("ksat")

# Filter for the 5 specific stations mentioned in the problem
station_names = ["Svalbard", "Punta Arenas", "Hartebeesthoek", "Awarua", "Athens"]
ksat_stations = [s for s in all_ksat if s.get_name() in station_names]
elapsed = time.time() - start_time
print(f"Loaded {len(ksat_stations)} KSAT ground stations in {elapsed:.2f} seconds:")
for station in ksat_stations:
    print(f"  - {station.get_name()}")

# Create 3D constellation visualization
print("\nPropagating Umbra satellites for one orbit each...")
start_time = time.time()
for prop in umbra_props:
    orbital_period = bh.orbital_period(prop.semi_major_axis)
    prop.propagate_to(prop.epoch + orbital_period)

fig_3d = bh.plot_trajectory_3d(
    [
        {
            "trajectory": prop.trajectory,
            "mode": "markers",
            "size": 2,
            "label": prop.get_name(),
        }
        for prop in umbra_props
    ],
    units="km",
    show_body=True,
    texture="natural_earth_50m",
    backend="plotly",
    view_azimuth=45.0,
    view_elevation=30.0,
    view_distance=2.0,
)
elapsed = time.time() - start_time
print(f"Created 3D visualization in {elapsed:.2f} seconds.")

# Reset propagators and compute 7-day access windows
print("\nComputing 7-day ground contacts...")
start_time = time.time()

# Reset all propagators
for prop in umbra_props:
    prop.reset()

# Define analysis period (7 days from first satellite's epoch)
epoch_start = umbra_props[0].epoch
epoch_end = epoch_start + 7 * 86400.0  # 7 days in seconds

# Propagate all satellites for 7 days
for prop in umbra_props:
    prop.propagate_to(epoch_end)

# Compute access windows with 5 degree minimum elevation
constraint = bh.ElevationConstraint(min_elevation_deg=5.0)
windows = bh.location_accesses(
    ksat_stations, umbra_props, epoch_start, epoch_end, constraint
)
elapsed = time.time() - start_time
print(f"Computed {len(windows)} contact windows in {elapsed:.2f} seconds.")

# Compute communication gaps per spacecraft
print("\nComputing communication gaps...")
start_time = time.time()

# Group windows by spacecraft
spacecraft_windows = {}
for window in windows:
    sat_name = window.satellite_name
    if sat_name not in spacecraft_windows:
        spacecraft_windows[sat_name] = []
    spacecraft_windows[sat_name].append(window)

# Sort each spacecraft's windows by start time
for windows in spacecraft_windows.values():
    windows.sort(key=lambda w: w.start.jd())

# Compute gaps between consecutive contacts
gaps = []
for sat_name, sat_windows in spacecraft_windows.items():
    for i in range(len(sat_windows) - 1):
        gap_start = sat_windows[i].end
        gap_end = sat_windows[i + 1].start
        gap_duration = gap_end - gap_start  # Duration in seconds

        gaps.append(
            {
                "spacecraft": sat_name,
                "gap_start": gap_start,
                "gap_end": gap_end,
                "duration": gap_duration,
                "last_station": sat_windows[i].location_name,
                "next_station": sat_windows[i + 1].location_name,
            }
        )

# Sort gaps by duration (longest first)
gaps.sort(key=lambda g: g["duration"], reverse=True)

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

# Print top 10 gaps
print("\n" + "=" * 100)
print("Top 10 Longest Communication Gaps")
print("=" * 100)
print(
    f"{'Spacecraft':<20} {'Start Time':<25} {'End Time':<25} {'Duration':>25} {'Last→Next Station':<30}"
)
print("-" * 100)
for i, gap in enumerate(gaps[:10]):
    start_str = str(gap["gap_start"]).split(".")[0]  # Remove fractional seconds
    end_str = str(gap["gap_end"]).split(".")[0]
    duration_str = bh.format_time_string(gap["duration"], short=False)
    station_str = f"{gap['last_station']} → {gap['next_station']}"
    print(
        f"{gap['spacecraft']:<20} {start_str:<25} {end_str:<25} {duration_str:>25} {station_str:<30}"
    )
print("=" * 100)

# Export top 10 gaps to CSV for documentation
csv_path = OUTDIR / f"{SCRIPT_NAME}_gaps.csv"
with open(csv_path, "w", newline="") as csvfile:
    writer = csv.writer(csvfile)
    writer.writerow(
        [
            "Spacecraft",
            "Gap Start (UTC)",
            "Gap End (UTC)",
            "Duration",
            "Last Station",
            "Next Station",
        ]
    )
    for gap in gaps[:10]:  # Only export top 10
        start_str = str(gap["gap_start"]).split(".")[0]
        end_str = str(gap["gap_end"]).split(".")[0]
        duration_str = bh.format_time_string(gap["duration"], short=True)
        writer.writerow(
            [
                gap["spacecraft"],
                start_str,
                end_str,
                duration_str,
                gap["last_station"],
                gap["next_station"],
            ]
        )
print(f"\n✓ Exported top 10 gaps to {csv_path}")

# Analyze gap distribution statistics
print("\nGap Distribution Statistics:")
# --8<-- [gap_statistics]
gap_durations = [g["duration"] for g in gaps]
mean_gap = np.mean(gap_durations)
median_gap = np.median(gap_durations)
min_gap = np.min(gap_durations)
max_gap = np.max(gap_durations)

print(f"  Mean: {bh.format_time_string(mean_gap)}")
print(f"  Median: {bh.format_time_string(median_gap)}")
print(f"  Min: {bh.format_time_string(min_gap)}")
print(f"  Max: {bh.format_time_string(max_gap)}")

# Create gap distribution histogram
print("\nCreating gap distribution histogram...")
gap_durations_hours = [d / 3600.0 for d in gap_durations]  # Convert to hours

fig_histogram = go.Figure(
    data=[
        go.Histogram(
            x=gap_durations_hours,
            nbinsx=40,
            marker_color="coral",
            marker_line_color="black",
            marker_line_width=1,
            hovertemplate="Gap Duration: %{x:.1f} hours<br>Count: %{y}<extra></extra>",
        )
    ]
)

fig_histogram.update_layout(
    title="Umbra Constellation Communication Gap Distribution (7-day period)",
    xaxis_title="Gap Duration (hours)",
    yaxis_title="Frequency",
    height=700,
    margin={"l": 60, "r": 40, "t": 80, "b": 60},
    annotations=[
        {
            "text": f"Mean: {bh.format_time_string(mean_gap)}<br>"
            f"Median: {bh.format_time_string(median_gap)}<br>"
            f"Max: {bh.format_time_string(max_gap)}",
            "xref": "paper",
            "yref": "paper",
            "x": 0.95,
            "y": 0.97,
            "xanchor": "right",
            "yanchor": "top",
            "showarrow": False,
            "bordercolor": "grey",
            "borderwidth": 1,
            "borderpad": 8,
        }
    ],
)

# Create cumulative distribution plot
print("\nCreating cumulative distribution plot...")

# Sort gap durations and compute cumulative percentages
sorted_gap_durations_hours = sorted(gap_durations_hours)
cumulative_percentages = [
    (i + 1) / len(sorted_gap_durations_hours) * 100
    for i in range(len(sorted_gap_durations_hours))
]

# Create cumulative distribution figure
fig_cumulative = go.Figure(
    data=[
        go.Scatter(
            x=sorted_gap_durations_hours,
            y=cumulative_percentages,
            mode="lines",
            line={"color": "steelblue", "width": 2.5},
            hovertemplate="Gap Duration: %{x:.2f} hours<br>Cumulative: %{y:.1f}%<extra></extra>",
        )
    ]
)

# Add reference lines for key percentiles
percentile_values = {
    25: sorted_gap_durations_hours[int(len(sorted_gap_durations_hours) * 0.25)],
    50: sorted_gap_durations_hours[int(len(sorted_gap_durations_hours) * 0.50)],
    75: sorted_gap_durations_hours[int(len(sorted_gap_durations_hours) * 0.75)],
    90: sorted_gap_durations_hours[int(len(sorted_gap_durations_hours) * 0.90)],
}

shapes = []
annotations = []
for percentile, value in percentile_values.items():
    # Add horizontal line
    shapes.append(
        {
            "type": "line",
            "x0": 0,
            "x1": value,
            "y0": percentile,
            "y1": percentile,
            "line": {"color": "rgba(128, 128, 128, 0.3)", "width": 1, "dash": "dash"},
        }
    )
    # Add vertical line
    shapes.append(
        {
            "type": "line",
            "x0": value,
            "x1": value,
            "y0": 0,
            "y1": percentile,
            "line": {"color": "rgba(128, 128, 128, 0.3)", "width": 1, "dash": "dash"},
        }
    )
    # Add annotation
    annotations.append(
        {
            "x": value,
            "y": percentile,
            "text": f"P{percentile}: {value:.2f}h",
            "showarrow": False,
            "xanchor": "left",
            "yanchor": "bottom",
            "xshift": 5,
            "yshift": 5,
            "font": {"size": 9, "color": "gray"},
        }
    )

fig_cumulative.update_layout(
    title="Umbra Constellation Cumulative Gap Distribution (7-day period)",
    xaxis_title="Gap Duration (hours)",
    yaxis_title="Cumulative Percentage (%)",
    height=700,
    margin={"l": 60, "r": 40, "t": 80, "b": 60},
    shapes=shapes,
    annotations=annotations,
    yaxis={"range": [0, 105]},
)

# Create ground track visualization for top 3 longest gaps
print("\nCreating ground track visualization for top 3 longest gaps...")
start_time = time.time()


# Get the top 3 gaps
top_3_gaps = gaps[:3]
gap_colors = ["red", "orange", "yellow"]

# For each gap, we need to extract the ground track segment during that period
# We'll need to get the propagator for each gap's spacecraft
gap_segments_all = []

for gap_idx, gap in enumerate(top_3_gaps):
    # Create label with spacecraft name and station transition
    gap_label = f"{gap['spacecraft']}, {gap['last_station']} → {gap['next_station']}"

    # Find the propagator for this gap's spacecraft
    sat_prop = None
    for prop in umbra_props:
        if prop.get_name() == gap["spacecraft"]:
            sat_prop = prop
            break

    if sat_prop is None:
        print(f"Warning: Could not find propagator for {gap['spacecraft']}")
        continue

    # Get states and epochs from the propagator's trajectory
    traj = sat_prop.trajectory
    states = traj.to_matrix()
    epochs = traj.epochs()

    # Extract ground track points during this gap period
    gap_lons = []
    gap_lats = []

    for i, ep in enumerate(epochs):
        if gap["gap_start"] <= ep <= gap["gap_end"]:
            # Convert to geodetic coordinates
            ecef_state = bh.state_eci_to_ecef(ep, states[i])
            lon, lat, alt = bh.position_ecef_to_geodetic(
                ecef_state[:3], bh.AngleFormat.RADIANS
            )
            gap_lons.append(np.degrees(lon))
            gap_lats.append(np.degrees(lat))

    # Split ground track at antimeridian crossings for proper plotting
    segments = bh.split_ground_track_at_antimeridian(gap_lons, gap_lats)

    # Interpolate segments to the edge at ±180° to avoid visual gaps
    # When a track crosses the antimeridian, add interpolated edge points
    interpolated_segments = []
    for seg_idx, (lon_seg, lat_seg) in enumerate(segments):
        lon_list = list(lon_seg)
        lat_list = list(lat_seg)

        # Check if this segment needs edge interpolation
        # If first point is not at edge but previous segment exists, interpolate start
        if seg_idx > 0 and len(lon_list) > 0:
            prev_lon, prev_lat = (
                segments[seg_idx - 1][0][-1],
                segments[seg_idx - 1][1][-1],
            )
            curr_lon, curr_lat = lon_list[0], lat_list[0]

            # Check if there was an antimeridian crossing
            if abs(curr_lon - prev_lon) > 180:
                # Determine which edge we're interpolating to
                if prev_lon > 0:  # Previous segment ends in positive lon
                    edge_lon = 180.0
                    # Linear interpolation to find latitude at 180°
                    t = (edge_lon - prev_lon) / ((curr_lon + 360) - prev_lon)
                    edge_lat = prev_lat + t * (curr_lat - prev_lat)
                    # Prepend edge point to current segment
                    lon_list.insert(0, edge_lon)
                    lat_list.insert(0, edge_lat)
                else:  # Previous segment ends in negative lon
                    edge_lon = -180.0
                    # Linear interpolation to find latitude at -180°
                    t = (edge_lon - prev_lon) / ((curr_lon - 360) - prev_lon)
                    edge_lat = prev_lat + t * (curr_lat - prev_lat)
                    # Prepend edge point to current segment
                    lon_list.insert(0, edge_lon)
                    lat_list.insert(0, edge_lat)

        # Check if segment should be extended to the edge at the end
        if seg_idx < len(segments) - 1 and len(lon_list) > 0:
            curr_lon, curr_lat = lon_list[-1], lat_list[-1]
            next_lon, next_lat = (
                segments[seg_idx + 1][0][0],
                segments[seg_idx + 1][1][0],
            )

            # Check if there will be an antimeridian crossing
            if abs(next_lon - curr_lon) > 180:
                # Determine which edge we're interpolating to
                if curr_lon > 0:  # Current segment ends in positive lon
                    edge_lon = 180.0
                    # Linear interpolation to find latitude at 180°
                    t = (edge_lon - curr_lon) / ((next_lon + 360) - curr_lon)
                    edge_lat = curr_lat + t * (next_lat - curr_lat)
                    # Append edge point to current segment
                    lon_list.append(edge_lon)
                    lat_list.append(edge_lat)
                else:  # Current segment ends in negative lon
                    edge_lon = -180.0
                    # Linear interpolation to find latitude at -180°
                    t = (edge_lon - curr_lon) / ((next_lon - 360) - curr_lon)
                    edge_lat = curr_lat + t * (next_lat - curr_lat)
                    # Append edge point to current segment
                    lon_list.append(edge_lon)
                    lat_list.append(edge_lat)

        interpolated_segments.append((lon_list, lat_list))

    gap_segments_all.append(
        {
            "segments": interpolated_segments,
            "color": gap_colors[gap_idx],
            "label": gap_label,
            "spacecraft": gap["spacecraft"],
            "duration": gap["duration"],
        }
    )

    print(
        f"  Gap {gap_idx + 1}: {len(gap_lons)} points, {len(interpolated_segments)} segments"
    )

# Create base plot with KSAT stations only (matplotlib backend for SVG output)

fig_groundtrack = bh.plot_groundtrack(
    ground_stations=[{"stations": ksat_stations, "color": "blue", "alpha": 0.2}],
    gs_cone_altitude=500e3,  # Approximate Umbra altitude
    gs_min_elevation=5.0,
    basemap="stock",
    show_borders=False,
    show_coastlines=False,
    backend="matplotlib",
)

# Plot each gap segment
ax = fig_groundtrack.get_axes()[0]
for gap_data in gap_segments_all:
    for i, (lon_seg, lat_seg) in enumerate(gap_data["segments"]):
        ax.plot(
            lon_seg,
            lat_seg,
            color=gap_data["color"],
            linewidth=1.5,
            transform=ccrs.Geodetic(),
            zorder=10,
            label=gap_data["label"] if i == 0 else "",
        )

# Add legend with high zorder to render on top
legend = ax.legend(loc="lower left", fontsize=10)
legend.set_zorder(100)

# Add title
ax.set_title(
    "Umbra Constellation: Top 3 Longest Communication Gaps\n"
    "KSAT Network (5 stations, 5° elevation)",
    fontsize=12,
)

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