# The Clear BSD License
#
# Copyright (c) 2026 Tobias Heibges
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted (subject to the limitations in the disclaimer
# below) provided that the following conditions are met:
#
# * Redistributions of source code must retain the above copyright notice,
# this list of conditions and the following disclaimer.
#
# * Redistributions in binary form must reproduce the above copyright
# notice, this list of conditions and the following disclaimer in the
# documentation and/or other materials provided with the distribution.
#
# * Neither the name of the copyright holder nor the names of its
# contributors may be used to endorse or promote products derived from this
# software without specific prior written permission.
#
# NO EXPRESS OR IMPLIED LICENSES TO ANY PARTY'S PATENT RIGHTS ARE GRANTED BY
# THIS LICENSE. THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND
# CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
# PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR
# CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
# EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
# PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR
# BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER
# IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
# POSSIBILITY OF SUCH DAMAGE.
"""
Ring Detector Visualization Module.
Provides plotting functions for constant-size detector ring analysis,
including geometry schematics, intensity heatmaps, timing distributions,
and spread analysis.
Functions are organized into:
- Geometry plots: Ring layout, side views, 3D visualization
- Analysis plots: Intensity/timing heatmaps, radial profiles, distributions
"""
from pathlib import Path
from typing import Optional, Union
import matplotlib.pyplot as plt
import numpy as np
from matplotlib.collections import PatchCollection
from matplotlib.patches import Wedge
from ..detectors.constant_size_rings import ConstantSizeDetectorRings
[docs]
def plot_geometry_analysis(
rings: ConstantSizeDetectorRings,
output_path: Optional[Union[str, Path]] = None,
figsize: tuple = (16, 14),
dpi: int = 150,
) -> plt.Figure:
"""
Create 4-panel geometry analysis figure.
Panels:
(a) Full geometry schematic (Earth + detector sphere + sample rings)
(b) Zoomed view showing no-shadowing geometry with sight lines
(c) Distance vs elevation curve
(d) Angular width and distance vs ring index
Parameters
----------
rings : ConstantSizeDetectorRings
Detector ring configuration
output_path : str or Path, optional
If provided, save figure to this path
figsize : tuple
Figure size in inches
dpi : int
Resolution for saved figure
Returns
-------
Figure
Matplotlib figure object
"""
fig, axes = plt.subplots(2, 2, figsize=figsize)
earth_radius = rings.earth_radius
detector_sphere_radius = rings.detector_sphere_radius
detector_altitude = rings.detector_altitude
n_rings = rings.n_rings
ring_centers_deg = rings.ring_centers_deg
ring_distances = rings.ring_distances
ring_boundaries_deg = rings.ring_boundaries_deg
detector_half_width = rings.detector_half_width
# Compute max horizontal distance
ring_horiz_dists = rings.get_ring_horizontal_distances() / 1000 # km
max_horiz_dist = max(ring_horiz_dists)
# Panel (a): Full geometry schematic
ax_full = axes[0, 0]
theta_arc = np.linspace(-10, 100, 200)
earth_x = earth_radius * np.sin(np.radians(theta_arc)) / 1000
earth_z = earth_radius * np.cos(np.radians(theta_arc)) / 1000 - earth_radius / 1000
ax_full.plot(earth_x, earth_z, "b-", linewidth=2, label="Earth surface")
ax_full.fill_between(earth_x, earth_z, -7000, alpha=0.1, color="blue")
det_x = detector_sphere_radius * np.sin(np.radians(theta_arc)) / 1000
det_z = (
detector_sphere_radius * np.cos(np.radians(theta_arc)) / 1000
- earth_radius / 1000
)
ax_full.plot(
det_x,
det_z,
"g-",
linewidth=1.5,
label=f"Detector sphere ({detector_altitude/1000:.0f} km)",
)
# Draw sample detector rings (every 2nd ring)
for i in range(0, min(n_rings, 12), 2):
theta_c = ring_centers_deg[i]
dist_c = ring_distances[i]
x_c = dist_c * np.cos(np.radians(theta_c)) / 1000
z_c = dist_c * np.sin(np.radians(theta_c)) / 1000
tx = -np.sin(np.radians(theta_c))
tz = np.cos(np.radians(theta_c))
hw_km = detector_half_width / 1000
x1, z1 = x_c - hw_km * tx, z_c - hw_km * tz
x2, z2 = x_c + hw_km * tx, z_c + hw_km * tz
color = plt.cm.viridis(i / min(n_rings, 12))
ax_full.plot([x1, x2], [z1, z2], "-", color=color, linewidth=3)
if i % 4 == 0:
ax_full.text(x_c + 5, z_c + 2, f"R{i}", fontsize=8, color=color)
# Sight lines
for i in range(0, min(4, len(ring_boundaries_deg)), 2):
elev = ring_boundaries_deg[i]
dist = rings.distance_at_elevation(elev)
x_end = dist * np.cos(np.radians(elev)) / 1000
z_end = dist * np.sin(np.radians(elev)) / 1000
ax_full.plot([0, x_end], [0, z_end], "r--", alpha=0.4, linewidth=0.8)
ax_full.scatter([0], [0], c="red", s=100, marker="*", zorder=10, label="Origin")
ax_full.set_xlabel("Horizontal Distance (km)", fontsize=11)
ax_full.set_ylabel("Altitude (km)", fontsize=11)
ax_full.set_title(
"(a) Full Geometry: Earth + Detector Sphere", fontsize=12, fontweight="bold"
)
ax_full.legend(loc="upper right", fontsize=9)
ax_full.set_xlim(-50, max_horiz_dist * 0.3)
ax_full.set_ylim(-50, 150)
ax_full.set_aspect("equal")
ax_full.grid(True, alpha=0.3)
# Panel (b): Zoomed view showing no-shadowing geometry
ax_zoom = axes[0, 1]
for i in range(min(5, n_rings)):
theta_c = ring_centers_deg[i]
dist_c = ring_distances[i]
x_c = dist_c * np.cos(np.radians(theta_c)) / 1000
z_c = dist_c * np.sin(np.radians(theta_c)) / 1000
tx = -np.sin(np.radians(theta_c))
tz = np.cos(np.radians(theta_c))
hw_km = detector_half_width / 1000
x1, z1 = x_c - hw_km * tx, z_c - hw_km * tz
x2, z2 = x_c + hw_km * tx, z_c + hw_km * tz
color = plt.cm.tab10(i)
ax_zoom.plot(
[x1, x2], [z1, z2], "-", color=color, linewidth=4, label=f"Ring {i}"
)
ax_zoom.plot([0, x1], [0, z1], "--", color=color, alpha=0.5, linewidth=1)
ax_zoom.plot([0, x2], [0, z2], "--", color=color, alpha=0.5, linewidth=1)
# Normal arrow (pointing toward origin)
arrow_len = 3
dist_to_origin = np.sqrt(x_c**2 + z_c**2)
if dist_to_origin > 0:
nx, nz = -x_c / dist_to_origin, -z_c / dist_to_origin
else:
nx, nz = 0, -1
ax_zoom.annotate(
"",
xy=(x_c + nx * arrow_len, z_c + nz * arrow_len),
xytext=(x_c, z_c),
arrowprops=dict(arrowstyle="->", color="darkgreen", lw=1.5),
)
ax_zoom.scatter([0], [0], c="red", s=150, marker="*", zorder=10, label="Origin")
ax_zoom.set_xlabel("Horizontal Distance (km)", fontsize=11)
ax_zoom.set_ylabel("Altitude (km)", fontsize=11)
ax_zoom.set_title(
"(b) No-Shadowing Geometry: Adjacent Rings Touch",
fontsize=12,
fontweight="bold",
)
ax_zoom.legend(loc="upper right", fontsize=9)
ax_zoom.set_xlim(-5, 80)
ax_zoom.set_ylim(-5, 50)
ax_zoom.set_aspect("equal")
ax_zoom.grid(True, alpha=0.3)
# Panel (c): Distance vs elevation curve
ax_dist = axes[1, 0]
elev_range = np.linspace(90, ring_boundaries_deg[-1], 200)
distances = [rings.distance_at_elevation(e) / 1000 for e in elev_range]
ax_dist.plot(elev_range, distances, "b-", linewidth=2)
ax_dist.scatter(
ring_centers_deg,
ring_distances / 1000,
c="red",
s=30,
zorder=5,
label="Ring centers",
)
ax_dist.set_xlabel("Elevation Angle (degrees)", fontsize=11)
ax_dist.set_ylabel("Distance from Origin (km)", fontsize=11)
ax_dist.set_title(
"(c) Distance to Detector Sphere vs Elevation", fontsize=12, fontweight="bold"
)
ax_dist.axhline(
y=detector_altitude / 1000,
color="green",
linestyle="--",
alpha=0.7,
label=f"Min distance ({detector_altitude/1000:.0f} km)",
)
ax_dist.legend(fontsize=9)
ax_dist.grid(True, alpha=0.3)
ax_dist.invert_xaxis()
# Panel (d): Angular width and distance vs ring index
ax_ring = axes[1, 1]
angular_widths = np.array([rings.angular_width_at_ring(i) for i in range(n_rings)])
ax_ring.bar(
np.arange(n_rings),
angular_widths,
color="steelblue",
alpha=0.7,
label="Angular width",
)
ax_ring.set_xlabel("Ring Index", fontsize=11)
ax_ring.set_ylabel("Angular Width (degrees)", fontsize=11, color="steelblue")
ax_ring.tick_params(axis="y", labelcolor="steelblue")
ax_ring2 = ax_ring.twinx()
ax_ring2.plot(
np.arange(n_rings), ring_distances / 1000, "r-o", markersize=4, label="Distance"
)
ax_ring2.set_ylabel("Distance from Origin (km)", fontsize=11, color="red")
ax_ring2.tick_params(axis="y", labelcolor="red")
ax_ring.set_title(
f"(d) Ring Properties ({n_rings} rings, {rings.detector_radial_size/1000:.0f} km each)",
fontsize=12,
fontweight="bold",
)
ax_ring.grid(True, alpha=0.3)
lines1, labels1 = ax_ring.get_legend_handles_labels()
lines2, labels2 = ax_ring2.get_legend_handles_labels()
ax_ring.legend(lines1 + lines2, labels1 + labels2, loc="upper right", fontsize=9)
fig.suptitle(
f"Constant-Size Detector Ring Geometry Analysis\n"
f"({rings.detector_radial_size/1000:.0f} km detectors, {detector_altitude/1000:.0f} km altitude, "
f"{n_rings} rings from {ring_boundaries_deg[0]:.1f}\u00b0 to {ring_boundaries_deg[-1]:.1f}\u00b0)",
fontsize=13,
fontweight="bold",
)
fig.tight_layout()
if output_path:
fig.savefig(output_path, dpi=dpi, bbox_inches="tight")
return fig
[docs]
def plot_ring_side_view(
rings: ConstantSizeDetectorRings,
output_path: Optional[Union[str, Path]] = None,
figsize: tuple = (16, 8),
dpi: int = 150,
) -> plt.Figure:
"""
Create side view cross-section of detector rings.
Shows all rings with constant physical size and normals pointing toward origin.
Parameters
----------
rings : ConstantSizeDetectorRings
Detector ring configuration
output_path : str or Path, optional
If provided, save figure to this path
figsize : tuple
Figure size in inches
dpi : int
Resolution for saved figure
Returns
-------
Figure
Matplotlib figure object
"""
fig, (ax_side1, ax_side2) = plt.subplots(1, 2, figsize=figsize)
n_rings = rings.n_rings
ring_centers_deg = rings.ring_centers_deg
ring_distances = rings.ring_distances
ring_boundaries_deg = rings.ring_boundaries_deg
detector_half_width = rings.detector_half_width
for ax_side in [ax_side1, ax_side2]:
for ring_idx in range(n_rings):
theta_c = ring_centers_deg[ring_idx]
dist_c = ring_distances[ring_idx]
x_center = dist_c * np.cos(np.radians(theta_c)) / 1000
z_center = dist_c * np.sin(np.radians(theta_c)) / 1000
dist_to_origin = np.sqrt(x_center**2 + z_center**2)
if dist_to_origin > 0:
nx, nz = -x_center / dist_to_origin, -z_center / dist_to_origin
else:
nx, nz = 0, -1
tx, tz = -nz, nx
hw_km = detector_half_width / 1000
p1x = x_center - hw_km * tx
p1z = z_center - hw_km * tz
p2x = x_center + hw_km * tx
p2z = z_center + hw_km * tz
ax_side.plot([p1x, p2x], [p1z, p2z], "b-", linewidth=2)
line_freq = max(1, n_rings // 10)
if ring_idx % line_freq == 0:
ax_side.plot(
[0, x_center], [0, z_center], "g--", linewidth=0.5, alpha=0.5
)
arrow_freq = max(1, n_rings // 5)
if ring_idx % arrow_freq == arrow_freq // 2:
arrow_len = 3
ax_side.annotate(
"",
xy=(x_center + nx * arrow_len, z_center + nz * arrow_len),
xytext=(x_center, z_center),
arrowprops=dict(arrowstyle="->", color="darkgreen", lw=1.5),
)
# Draw detector sphere arc for reference
arc_elev = np.linspace(
ring_boundaries_deg[-1] - 5, ring_boundaries_deg[0] + 5, 100
)
arc_x_plot = [
rings.distance_at_elevation(e) * np.cos(np.radians(e)) / 1000
for e in arc_elev
]
arc_z_plot = [
rings.distance_at_elevation(e) * np.sin(np.radians(e)) / 1000
for e in arc_elev
]
ax_side.plot(
arc_x_plot,
arc_z_plot,
"c-",
linewidth=0.5,
alpha=0.3,
label="Detector sphere",
)
ax_side.axhline(y=0, color="k", linestyle="-", linewidth=0.5)
ax_side.scatter([0], [0], c="red", s=80, zorder=10, label="Origin (0,0)")
ax_side.set_xlabel("Horizontal Distance (km)", fontsize=11)
ax_side.set_ylabel("Altitude (km)", fontsize=11)
ax_side.grid(True, alpha=0.3)
ax_side.legend(loc="upper right")
ax_side1.set_title(f"Full View: All {n_rings} Rings", fontsize=12)
ax_side2.set_title("Equal Aspect Ratio (First 10 Rings)", fontsize=12)
ax_side2.set_aspect("equal")
zoom_rings = min(10, n_rings)
max_x_zoom = max(
rings.distance_at_elevation(ring_boundaries_deg[zoom_rings])
* np.cos(np.radians(ring_boundaries_deg[zoom_rings]))
/ 1000,
100,
)
max_z_zoom = max(
rings.distance_at_elevation(ring_centers_deg[0])
* np.sin(np.radians(ring_centers_deg[0]))
/ 1000,
50,
)
ax_side2.set_xlim(-5, max_x_zoom * 1.1)
ax_side2.set_ylim(-5, max_z_zoom * 1.2)
fig.suptitle(
f"Vertical Cross-Section: Constant {rings.detector_radial_size/1000:.0f} km Detectors (No Shadowing)",
fontsize=13,
fontweight="bold",
)
fig.tight_layout()
if output_path:
fig.savefig(output_path, dpi=dpi, bbox_inches="tight")
return fig
[docs]
def plot_ring_overview(
rings: ConstantSizeDetectorRings,
ring_intensities: np.ndarray,
detected_positions: Optional[np.ndarray] = None,
output_path: Optional[Union[str, Path]] = None,
figsize: tuple = (10, 10),
dpi: int = 150,
) -> plt.Figure:
"""
Create top-down ring overview with intensity coloring.
Parameters
----------
rings : ConstantSizeDetectorRings
Detector ring configuration
ring_intensities : ndarray
Total intensity per ring (shape: n_rings)
detected_positions : ndarray, optional
Ray detection positions (N, 3) for scatter overlay
output_path : str or Path, optional
If provided, save figure to this path
figsize : tuple
Figure size in inches
dpi : int
Resolution for saved figure
Returns
-------
Figure
Matplotlib figure object
"""
fig, ax = plt.subplots(figsize=figsize)
ring_horiz_dists = rings.get_ring_horizontal_distances() / 1000 # km
max_horiz_dist = max(ring_horiz_dists)
n_rings = rings.n_rings
ring_patches = []
ring_colors = []
for ring_idx in range(n_rings):
inner_r = ring_horiz_dists[ring_idx]
outer_r = ring_horiz_dists[ring_idx + 1]
wedge = Wedge((0, 0), outer_r, 0, 360, width=(outer_r - inner_r))
ring_patches.append(wedge)
ring_colors.append(
ring_intensities[ring_idx] if ring_idx < len(ring_intensities) else 0
)
collection = PatchCollection(ring_patches, cmap="hot", alpha=0.8)
collection.set_array(np.array(ring_colors))
collection.set_edgecolor("gray")
collection.set_linewidth(0.5)
ax.add_collection(collection)
plt.colorbar(collection, ax=ax, label="Total Intensity per Ring")
ax.scatter(
0, 0, c="blue", s=100, marker="x", linewidths=2, label="Zenith", zorder=10
)
if detected_positions is not None:
ax.scatter(
detected_positions[:, 0] / 1000,
detected_positions[:, 1] / 1000,
c="cyan",
s=1,
alpha=0.3,
label=f"Detected rays ({len(detected_positions)})",
)
max_plot_radius = max_horiz_dist * 1.1
ax.set_xlim(-max_plot_radius, max_plot_radius)
ax.set_ylim(-max_plot_radius, max_plot_radius)
ax.set_xlabel("X (km)", fontsize=11)
ax.set_ylabel("Y (km)", fontsize=11)
ax.set_title(
f"Elevation Ring Overview\n(Constant {rings.detector_radial_size/1000:.0f} km detector size, {n_rings} rings)",
fontsize=12,
fontweight="bold",
)
ax.set_aspect("equal")
ax.legend(loc="upper right")
ax.grid(True, alpha=0.3)
if output_path:
fig.savefig(output_path, dpi=dpi, bbox_inches="tight")
return fig
[docs]
def plot_ring_azimuth_heatmap(
rings: ConstantSizeDetectorRings,
intensity_map: np.ndarray,
ray_count_map: np.ndarray,
az_min: float = -10.0,
az_max: float = 10.0,
output_path: Optional[Union[str, Path]] = None,
figsize: tuple = (14, 6),
dpi: int = 150,
) -> plt.Figure:
"""
Create ring-azimuth heatmaps for intensity and ray count.
Parameters
----------
rings : ConstantSizeDetectorRings
Detector ring configuration
intensity_map : ndarray
Intensity per (ring, azimuth_bin), shape (n_rings, n_az_bins)
ray_count_map : ndarray
Ray count per (ring, azimuth_bin), shape (n_rings, n_az_bins)
az_min, az_max : float
Azimuth range in degrees
output_path : str or Path, optional
If provided, save figure to this path
figsize : tuple
Figure size in inches
dpi : int
Resolution for saved figure
Returns
-------
Figure
Matplotlib figure object
"""
fig, (ax_int, ax_count) = plt.subplots(1, 2, figsize=figsize)
n_rings = rings.n_rings
ring_centers_deg = rings.ring_centers_deg
def ring_to_elev(ring_idx):
idx = np.clip(np.asarray(ring_idx), 0, n_rings - 1).astype(int)
return ring_centers_deg[idx]
def elev_to_ring(elev):
return np.interp(elev, ring_centers_deg[::-1], np.arange(n_rings)[::-1])
# Intensity heatmap
im1 = ax_int.imshow(
intensity_map,
aspect="auto",
origin="lower",
cmap="hot",
extent=[az_min, az_max, 0, n_rings],
)
ax_int.axvline(x=0, color="cyan", linestyle="--", linewidth=1, alpha=0.7)
ax_int.set_xlabel("Azimuth (deg) - 0\u00b0 = beam direction", fontsize=11)
ax_int.set_ylabel("Ring Index", fontsize=11)
ax_int.set_title("Intensity by Ring and Azimuth", fontsize=12, fontweight="bold")
plt.colorbar(im1, ax=ax_int, label="Total Intensity")
ax_int_elev = ax_int.secondary_yaxis(
"right", functions=(ring_to_elev, elev_to_ring)
)
ax_int_elev.set_ylabel("Elevation (deg)", fontsize=11)
# Ray count heatmap
im2 = ax_count.imshow(
ray_count_map,
aspect="auto",
origin="lower",
cmap="viridis",
extent=[az_min, az_max, 0, n_rings],
)
ax_count.axvline(x=0, color="red", linestyle="--", linewidth=1, alpha=0.7)
ax_count.set_xlabel("Azimuth (deg) - 0\u00b0 = beam direction", fontsize=11)
ax_count.set_ylabel("Ring Index", fontsize=11)
ax_count.set_title("Ray Count by Ring and Azimuth", fontsize=12, fontweight="bold")
plt.colorbar(im2, ax=ax_count, label="Ray Count")
ax_count_elev = ax_count.secondary_yaxis(
"right", functions=(ring_to_elev, elev_to_ring)
)
ax_count_elev.set_ylabel("Elevation (deg)", fontsize=11)
fig.suptitle(
f"Ring-Azimuth Distribution (\u00b1{int(az_max)}\u00b0 around beam direction)",
fontsize=14,
fontweight="bold",
)
fig.tight_layout()
if output_path:
fig.savefig(output_path, dpi=dpi, bbox_inches="tight")
return fig
[docs]
def plot_timing_heatmap(
rings: ConstantSizeDetectorRings,
first_arrival_map: np.ndarray,
time_spread_map: np.ndarray,
global_first_ns: float,
az_min: float = -10.0,
az_max: float = 10.0,
output_path: Optional[Union[str, Path]] = None,
figsize: tuple = (14, 6),
dpi: int = 150,
) -> plt.Figure:
"""
Create timing heatmaps (first arrival and time spread).
Parameters
----------
rings : ConstantSizeDetectorRings
Detector ring configuration
first_arrival_map : ndarray
First arrival time per bin in ns, shape (n_rings, n_az_bins)
time_spread_map : ndarray
Time spread (10-90%) per bin in ns, shape (n_rings, n_az_bins)
global_first_ns : float
Global first arrival time in ns (for relative timing)
az_min, az_max : float
Azimuth range in degrees
output_path : str or Path, optional
If provided, save figure to this path
figsize : tuple
Figure size in inches
dpi : int
Resolution for saved figure
Returns
-------
Figure
Matplotlib figure object
"""
fig, (ax_first, ax_spread) = plt.subplots(1, 2, figsize=figsize)
n_rings = rings.n_rings
ring_centers_deg = rings.ring_centers_deg
def ring_to_elev(ring_idx):
idx = np.clip(np.asarray(ring_idx), 0, n_rings - 1).astype(int)
return ring_centers_deg[idx]
def elev_to_ring(elev):
return np.interp(elev, ring_centers_deg[::-1], np.arange(n_rings)[::-1])
# First arrival time (relative to global first)
first_arrival_rel = first_arrival_map - global_first_ns
im1 = ax_first.imshow(
first_arrival_rel,
aspect="auto",
origin="lower",
cmap="plasma",
extent=[az_min, az_max, 0, n_rings],
)
ax_first.axvline(x=0, color="cyan", linestyle="--", linewidth=1, alpha=0.7)
ax_first.set_xlabel("Azimuth (deg) - 0\u00b0 = beam direction", fontsize=11)
ax_first.set_ylabel("Ring Index", fontsize=11)
ax_first.set_title("First Arrival Time per Bin", fontsize=12, fontweight="bold")
plt.colorbar(im1, ax=ax_first, label="First Arrival (ns from global first)")
ax_first_elev = ax_first.secondary_yaxis(
"right", functions=(ring_to_elev, elev_to_ring)
)
ax_first_elev.set_ylabel("Elevation (deg)", fontsize=11)
# Time spread
im2 = ax_spread.imshow(
time_spread_map,
aspect="auto",
origin="lower",
cmap="viridis",
extent=[az_min, az_max, 0, n_rings],
)
ax_spread.axvline(x=0, color="red", linestyle="--", linewidth=1, alpha=0.7)
ax_spread.set_xlabel("Azimuth (deg) - 0\u00b0 = beam direction", fontsize=11)
ax_spread.set_ylabel("Ring Index", fontsize=11)
ax_spread.set_title(
"Time Spread (10-90%) by Ring and Azimuth", fontsize=12, fontweight="bold"
)
plt.colorbar(im2, ax=ax_spread, label="Time Spread (ns)")
ax_spread_elev = ax_spread.secondary_yaxis(
"right", functions=(ring_to_elev, elev_to_ring)
)
ax_spread_elev.set_ylabel("Elevation (deg)", fontsize=11)
fig.suptitle(
f"Timing Distribution (\u00b1{int(az_max)}\u00b0 around beam direction)",
fontsize=14,
fontweight="bold",
)
fig.tight_layout()
if output_path:
fig.savefig(output_path, dpi=dpi, bbox_inches="tight")
return fig
[docs]
def plot_per_ring_timing_distributions(
rings: ConstantSizeDetectorRings,
ring_data: list,
az_limit: float = 10.0,
output_path: Optional[Union[str, Path]] = None,
figsize_per_col: float = 4.0,
figsize_per_row: float = 3.5,
dpi: int = 150,
) -> Optional[plt.Figure]:
"""
Create per-ring timing distribution histograms.
Parameters
----------
rings : ConstantSizeDetectorRings
Detector ring configuration
ring_data : list of dict
List of dicts with keys:
- ring_idx: int
- az_center: float (azimuth center in degrees)
- times_rel_ns: ndarray (times relative to bin first arrival, ns)
- intensities: ndarray (ray intensities)
- dist_km: float
az_limit : float
Azimuth limit for title
output_path : str or Path, optional
If provided, save figure to this path
figsize_per_col, figsize_per_row : float
Figure size multipliers
dpi : int
Resolution for saved figure
Returns
-------
Figure or None
Matplotlib figure object, or None if no data
"""
if not ring_data:
return None
n_plots = len(ring_data)
n_cols = min(4, n_plots)
n_rows = (n_plots + n_cols - 1) // n_cols
fig, axes = plt.subplots(
n_rows, n_cols, figsize=(figsize_per_col * n_cols, figsize_per_row * n_rows)
)
if n_plots == 1:
axes = np.array([[axes]])
elif n_rows == 1:
axes = axes.reshape(1, -1)
for i, info in enumerate(ring_data):
row, col = divmod(i, n_cols)
ax = axes[row, col]
times_rel = info["times_rel_ns"]
intensities = info["intensities"]
n_rays = len(times_rel)
if n_rays > 0 and np.sum(intensities) > 0:
counts, edges = np.histogram(times_rel, bins=30, weights=intensities)
centers = (edges[:-1] + edges[1:]) / 2
ax.fill_between(centers, counts, alpha=0.6, color="steelblue", step="mid")
ax.step(centers, counts, where="mid", color="steelblue", linewidth=1.5)
# Compute percentiles
sorted_idx = np.argsort(times_rel)
sorted_t = times_rel[sorted_idx]
sorted_w = intensities[sorted_idx]
cumsum = np.cumsum(sorted_w)
cumsum_norm = cumsum / cumsum[-1]
t10 = sorted_t[np.searchsorted(cumsum_norm, 0.10)]
t90 = sorted_t[np.searchsorted(cumsum_norm, 0.90)]
spread = t90 - t10
ax.axvline(t10, color="red", linestyle="--", linewidth=1.5, alpha=0.8)
ax.axvline(t90, color="red", linestyle="--", linewidth=1.5, alpha=0.8)
ax.axvspan(t10, t90, alpha=0.15, color="red")
ax.set_title(
f'Ring {info["ring_idx"]}, az={info["az_center"]:.1f}\u00b0\n'
f'd={info["dist_km"]:.0f}km, n={n_rays}, \u0394t={spread:.1f}ns',
fontsize=10,
)
else:
ax.set_title(f'Ring {info["ring_idx"]} (no data)', fontsize=10)
ax.set_xlabel("Time (ns)", fontsize=9)
ax.set_ylabel("Intensity", fontsize=9)
ax.grid(True, alpha=0.3)
ax.set_xlim(left=0)
# Hide unused axes
for i in range(n_plots, n_rows * n_cols):
row, col = divmod(i, n_cols)
axes[row, col].axis("off")
fig.suptitle(
f"Per-Ring Timing Distributions (best azimuth bin, \u00b1{az_limit:.0f}\u00b0)",
fontsize=14,
fontweight="bold",
)
fig.tight_layout()
if output_path:
fig.savefig(output_path, dpi=dpi, bbox_inches="tight")
return fig
[docs]
def plot_spread_analysis(
rings: ConstantSizeDetectorRings,
ring_stats: list,
output_path: Optional[Union[str, Path]] = None,
figsize: tuple = (14, 5),
dpi: int = 150,
) -> Optional[plt.Figure]:
"""
Create per-ring spread analysis (azimuthal and timing).
Parameters
----------
rings : ConstantSizeDetectorRings
Detector ring configuration
ring_stats : list of dict
List of dicts with keys:
- ring_idx: int
- n_rays: int
- az_std: float (azimuth std in degrees)
- time_spread_ns: float
- dist_km: float
output_path : str or Path, optional
If provided, save figure to this path
figsize : tuple
Figure size in inches
dpi : int
Resolution for saved figure
Returns
-------
Figure or None
Matplotlib figure object, or None if no data
"""
if not ring_stats:
return None
fig, (ax_az, ax_time) = plt.subplots(1, 2, figsize=figsize)
ring_indices = np.array([s["ring_idx"] for s in ring_stats])
az_spreads = [s["az_std"] for s in ring_stats]
time_spreads = [s["time_spread_ns"] for s in ring_stats]
ray_counts = [s["n_rays"] for s in ring_stats]
colors = np.log10(np.array(ray_counts) + 1)
# Azimuthal spread vs ring index
sc1 = ax_az.scatter(
ring_indices, az_spreads, c=colors, cmap="viridis", s=80, edgecolors="black"
)
ax_az.plot(ring_indices, az_spreads, "b-", alpha=0.3)
ax_az.set_xlabel("Ring Index", fontsize=11)
ax_az.set_ylabel("Azimuthal Spread (std, degrees)", fontsize=11)
ax_az.set_title("Azimuthal Spread vs Ring", fontsize=12, fontweight="bold")
ax_az.grid(True, alpha=0.3)
cbar1 = plt.colorbar(sc1, ax=ax_az)
cbar1.set_label("log\u2081\u2080(ray count)")
# Add elevation axis on top for ax_az
ax_az_top = ax_az.twiny()
ax_az_top.set_xlim(ax_az.get_xlim())
tick_indices = ring_indices[:: max(1, len(ring_indices) // 5)]
tick_elevs = [
rings.ring_centers_deg[int(i)] for i in tick_indices if i < rings.n_rings
]
ax_az_top.set_xticks(tick_indices[: len(tick_elevs)])
ax_az_top.set_xticklabels([f"{e:.1f}°" for e in tick_elevs])
ax_az_top.set_xlabel("Elevation", fontsize=10)
# Time spread vs ring index
sc2 = ax_time.scatter(
ring_indices, time_spreads, c=colors, cmap="viridis", s=80, edgecolors="black"
)
ax_time.plot(ring_indices, time_spreads, "r-", alpha=0.3)
ax_time.set_xlabel("Ring Index", fontsize=11)
ax_time.set_ylabel("Time Spread (10-90%, ns)", fontsize=11)
ax_time.set_title(
"Per-Bin Time Spread vs Ring (median over azimuth)",
fontsize=12,
fontweight="bold",
)
ax_time.grid(True, alpha=0.3)
cbar2 = plt.colorbar(sc2, ax=ax_time)
cbar2.set_label("log\u2081\u2080(ray count)")
# Add elevation axis on top for ax_time
ax_time_top = ax_time.twiny()
ax_time_top.set_xlim(ax_time.get_xlim())
ax_time_top.set_xticks(tick_indices[: len(tick_elevs)])
ax_time_top.set_xticklabels([f"{e:.1f}°" for e in tick_elevs])
ax_time_top.set_xlabel("Elevation", fontsize=10)
fig.suptitle("Per-Ring Spread Analysis", fontsize=14, fontweight="bold")
fig.tight_layout()
if output_path:
fig.savefig(output_path, dpi=dpi, bbox_inches="tight")
return fig
[docs]
def plot_polar_irradiance(
rings: ConstantSizeDetectorRings,
irradiance_map: np.ndarray,
n_azimuth_bins: int,
output_path: Optional[Union[str, Path]] = None,
figsize: tuple = (10, 10),
dpi: int = 150,
) -> plt.Figure:
"""
Create polar heatmap of irradiance.
Parameters
----------
rings : ConstantSizeDetectorRings
Detector ring configuration
irradiance_map : ndarray
Irradiance per (ring, azimuth_bin), shape (n_rings, n_azimuth_bins)
n_azimuth_bins : int
Number of azimuth bins
output_path : str or Path, optional
If provided, save figure to this path
figsize : tuple
Figure size in inches
dpi : int
Resolution for saved figure
Returns
-------
Figure
Matplotlib figure object
"""
fig = plt.figure(figsize=figsize)
ax = fig.add_subplot(111, projection="polar")
ring_horiz_dists = rings.get_ring_horizontal_distances() / 1000 # km
n_rings = rings.n_rings
theta = np.linspace(-np.pi, np.pi, n_azimuth_bins + 1)
r = np.array(list(ring_horiz_dists[: n_rings + 1]))
R, Theta = np.meshgrid(r, theta)
c = ax.pcolormesh(Theta, R, irradiance_map.T, cmap="hot", shading="auto")
plt.colorbar(
c, ax=ax, label="Irradiance (W/m\u00b2)", orientation="horizontal", pad=0.08
)
ax.set_theta_zero_location("E")
ax.set_theta_direction(1)
ax.set_title(
"Irradiance vs Horizontal Distance and Azimuth\n(0\u00b0 = beam direction \u2192)",
fontsize=12,
fontweight="bold",
pad=20,
)
if output_path:
fig.savefig(output_path, dpi=dpi, bbox_inches="tight")
return fig
[docs]
def compute_ring_azimuth_maps(
ray_positions: np.ndarray,
ray_times: np.ndarray,
ray_intensities: np.ndarray,
rings: ConstantSizeDetectorRings,
az_min: float = -10.0,
az_max: float = 10.0,
n_az_bins: int = 40,
) -> dict:
"""
Compute 2D maps of intensity, ray count, and timing by ring and azimuth.
Parameters
----------
ray_positions : ndarray
Ray detection positions (N, 3)
ray_times : ndarray
Ray detection times (N,)
ray_intensities : ndarray
Ray intensities (N,)
rings : ConstantSizeDetectorRings
Detector ring configuration
az_min, az_max : float
Azimuth range in degrees
n_az_bins : int
Number of azimuth bins
Returns
-------
dict
Dictionary with keys:
- intensity_map: (n_rings, n_az_bins) total intensity
- ray_count_map: (n_rings, n_az_bins) ray count
- first_arrival_map: (n_rings, n_az_bins) first arrival time in ns
- time_spread_map: (n_rings, n_az_bins) time spread (10-90%) in ns
- ray_ring_indices: (N,) ring index per ray
- ray_azimuth_deg: (N,) azimuth per ray in degrees
"""
n_rings = rings.n_rings
ring_boundaries_deg = rings.ring_boundaries_deg
# Compute elevation and azimuth for each ray
horizontal_dist = np.sqrt(ray_positions[:, 0] ** 2 + ray_positions[:, 1] ** 2)
vertical_dist = ray_positions[:, 2]
ray_elevation_deg = np.degrees(np.arctan2(vertical_dist, horizontal_dist))
ray_azimuth_deg = np.degrees(np.arctan2(ray_positions[:, 1], ray_positions[:, 0]))
# Assign rays to rings
ray_ring_indices = (
np.searchsorted(-ring_boundaries_deg[:-1], -ray_elevation_deg, side="right") - 1
)
ray_ring_indices = np.clip(ray_ring_indices, 0, n_rings - 1)
# Initialize maps
intensity_map = np.zeros((n_rings, n_az_bins))
ray_count_map = np.zeros((n_rings, n_az_bins))
first_arrival_map = np.full((n_rings, n_az_bins), np.nan)
time_spread_map = np.full((n_rings, n_az_bins), np.nan)
az_bin_width = (az_max - az_min) / n_az_bins
for ring_idx in range(n_rings):
for az_bin in range(n_az_bins):
az_lo = az_min + az_bin * az_bin_width
az_hi = az_lo + az_bin_width
mask = (
(ray_ring_indices == ring_idx)
& (ray_azimuth_deg >= az_lo)
& (ray_azimuth_deg < az_hi)
)
n_in_cell = np.sum(mask)
if n_in_cell > 0:
cell_times = ray_times[mask]
cell_int = ray_intensities[mask]
intensity_map[ring_idx, az_bin] = np.sum(cell_int)
ray_count_map[ring_idx, az_bin] = n_in_cell
bin_first = np.min(cell_times)
first_arrival_map[ring_idx, az_bin] = bin_first * 1e9 # ns
if n_in_cell >= 3 and np.sum(cell_int) > 0:
cell_times_rel = (cell_times - bin_first) * 1e9
sorted_idx = np.argsort(cell_times_rel)
sorted_t = cell_times_rel[sorted_idx]
sorted_w = cell_int[sorted_idx]
cumsum = np.cumsum(sorted_w)
cumsum_norm = cumsum / cumsum[-1]
t10 = sorted_t[np.searchsorted(cumsum_norm, 0.10)]
t90 = sorted_t[np.searchsorted(cumsum_norm, 0.90)]
time_spread_map[ring_idx, az_bin] = t90 - t10
return {
"intensity_map": intensity_map,
"ray_count_map": ray_count_map,
"first_arrival_map": first_arrival_map,
"time_spread_map": time_spread_map,
"ray_ring_indices": ray_ring_indices,
"ray_azimuth_deg": ray_azimuth_deg,
}
def plot_time_spread_comparison(
rings: ConstantSizeDetectorRings,
ring_stats: list,
source_position: tuple,
beam_direction: tuple,
divergence_angle_rad: float,
surface=None,
output_path: Optional[Union[str, Path]] = None,
figsize: tuple = (14, 5),
dpi: int = 150,
) -> Optional[plt.Figure]:
"""
Compare ray-traced time spread with geometric estimate.
The geometric estimate computes the path difference for rays at the edge
of the divergence cone reflecting through the beam footprint on the surface.
Parameters
----------
rings : ConstantSizeDetectorRings
Detector ring configuration
ring_stats : list of dict
Per-ring statistics with time_spread_ns and dist_km
source_position : tuple
Source (x, y, z) position in meters
beam_direction : tuple
Beam direction unit vector
divergence_angle_rad : float
Beam half-angle divergence in radians
surface : Surface, optional
Ocean surface for geometric estimate. If None, uses flat surface at z=0.
output_path : str or Path, optional
If provided, save figure to this path
figsize : tuple
Figure size in inches
dpi : int
Resolution for saved figure
Returns
-------
Figure or None
Matplotlib figure object, or None if no data
"""
from ..utilities import estimate_time_spread
if not ring_stats:
return None
fig, (ax_abs, ax_ratio) = plt.subplots(1, 2, figsize=figsize)
# Collect data
distances = []
raytraced_spreads = []
geometric_spreads = []
ring_indices = []
print("\n Time Spread Comparison Debug:")
print(f" Source: {source_position}")
print(f" Beam dir: {beam_direction}")
print(f" Divergence: {np.degrees(divergence_angle_rad):.2f} deg")
for stats in ring_stats:
ring_idx = stats["ring_idx"]
dist_km = stats["dist_km"]
raytraced_ns = stats["time_spread_ns"]
# Skip rings with zero or invalid time spread
if raytraced_ns <= 0 or ring_idx >= rings.n_rings:
continue
# Get detector position for this ring center
theta_c = rings.ring_centers_deg[ring_idx]
dist_m = rings.ring_distances[ring_idx]
# Detector position (azimuth = 0)
det_x = dist_m * np.cos(np.radians(theta_c))
det_y = 0.0
det_z = dist_m * np.sin(np.radians(theta_c))
detector_position = (det_x, det_y, det_z)
# Compute geometric estimate
result = estimate_time_spread(
source_position=source_position,
beam_direction=beam_direction,
divergence_angle=divergence_angle_rad,
detector_position=detector_position,
surface=surface,
)
# Debug output for first few rings
if len(distances) < 5:
print(
f" Ring {ring_idx}: dist={dist_km:.0f}km, "
f"raytraced={raytraced_ns:.2f}ns, geometric={result.time_spread_ns:.2f}ns, "
f"path_spread={result.path_spread:.1f}m"
)
distances.append(dist_km)
raytraced_spreads.append(raytraced_ns)
geometric_spreads.append(result.time_spread_ns)
ring_indices.append(ring_idx)
if not distances:
return None
distances = np.array(distances)
raytraced_spreads = np.array(raytraced_spreads)
geometric_spreads = np.array(geometric_spreads)
ring_indices = np.array(ring_indices)
# Left panel: Absolute comparison
ax_abs.scatter(
ring_indices,
raytraced_spreads,
c="blue",
s=60,
label="Ray-traced (10-90%)",
alpha=0.7,
)
ax_abs.scatter(
ring_indices,
geometric_spreads,
c="red",
s=60,
marker="^",
label="Geometric estimate",
alpha=0.7,
)
ax_abs.plot(ring_indices, geometric_spreads, "r--", alpha=0.5)
ax_abs.set_xlabel("Ring Index", fontsize=11)
ax_abs.set_ylabel("Time Spread (ns)", fontsize=11)
ax_abs.set_title(
"Time Spread: Ray-Traced vs Geometric", fontsize=12, fontweight="bold"
)
ax_abs.legend(loc="upper left", fontsize=9)
ax_abs.grid(True, alpha=0.3)
ax_abs.set_yscale("log")
# Add elevation axis on top for ax_abs
ax_abs_top = ax_abs.twiny()
ax_abs_top.set_xlim(ax_abs.get_xlim())
tick_indices = ring_indices[:: max(1, len(ring_indices) // 5)]
tick_elevs = [
rings.ring_centers_deg[int(i)] for i in tick_indices if i < rings.n_rings
]
ax_abs_top.set_xticks(tick_indices[: len(tick_elevs)])
ax_abs_top.set_xticklabels([f"{e:.1f}°" for e in tick_elevs])
ax_abs_top.set_xlabel("Elevation", fontsize=10)
# Right panel: Ratio
ratio = raytraced_spreads / geometric_spreads
valid_mask = geometric_spreads > 0
if np.any(valid_mask):
sc = ax_ratio.scatter(
ring_indices[valid_mask],
ratio[valid_mask],
c=distances[valid_mask],
cmap="viridis",
s=60,
alpha=0.7,
)
ax_ratio.axhline(
y=1.0, color="red", linestyle="--", linewidth=1.5, label="Perfect agreement"
)
# Add colorbar for distance
cbar = plt.colorbar(sc, ax=ax_ratio)
cbar.set_label("Distance (km)")
ax_ratio.set_xlabel("Ring Index", fontsize=11)
ax_ratio.set_ylabel("Ray-Traced / Geometric", fontsize=11)
ax_ratio.set_title("Ratio of Time Spreads", fontsize=12, fontweight="bold")
ax_ratio.legend(loc="upper right", fontsize=9)
ax_ratio.grid(True, alpha=0.3)
# Add elevation axis on top for ax_ratio
ax_ratio_top = ax_ratio.twiny()
ax_ratio_top.set_xlim(ax_ratio.get_xlim())
ax_ratio_top.set_xticks(tick_indices[: len(tick_elevs)])
ax_ratio_top.set_xticklabels([f"{e:.1f}°" for e in tick_elevs])
ax_ratio_top.set_xlabel("Elevation", fontsize=10)
fig.suptitle(
"Time Spread Comparison: Ray-Traced vs Single-Bounce Geometric Estimate",
fontsize=14,
fontweight="bold",
)
fig.tight_layout()
if output_path:
fig.savefig(output_path, dpi=dpi, bbox_inches="tight")
return fig
[docs]
def compute_constant_size_bin_stats(
ray_positions: np.ndarray,
ray_times: np.ndarray,
ray_intensities: np.ndarray,
rings: ConstantSizeDetectorRings,
az_bin_size_m: float = 10000.0,
az_range_deg: float = 10.0,
min_rays_per_bin: int = 3,
) -> list[dict]:
"""
Compute timing statistics for constant physical size bins.
Each bin has approximately the same physical size:
- Radial: detector_radial_size (from rings configuration)
- Azimuthal: az_bin_size_m (varies in angular size with distance)
Parameters
----------
ray_positions : ndarray
Ray detection positions (N, 3)
ray_times : ndarray
Ray detection times (N,)
ray_intensities : ndarray
Ray intensities (N,)
rings : ConstantSizeDetectorRings
Detector ring configuration
az_bin_size_m : float
Physical azimuthal bin size in meters (default: 10 km)
az_range_deg : float
Azimuth range in degrees (±this value)
min_rays_per_bin : int
Minimum rays required to compute statistics
Returns
-------
list of dict
Per-bin statistics with keys:
- ring_idx, az_bin_idx, n_az_bins
- az_center_deg, elev_center_deg
- distance_km
- n_rays, total_intensity
- first_arrival_ns (relative to global first)
- time_spread_ns (10-90 percentile)
- bin_area_m2
"""
n_rings = rings.n_rings
ring_boundaries_deg = rings.ring_boundaries_deg
# Compute elevation and azimuth for each ray
horizontal_dist = np.sqrt(ray_positions[:, 0] ** 2 + ray_positions[:, 1] ** 2)
vertical_dist = ray_positions[:, 2]
ray_elevation_deg = np.degrees(np.arctan2(vertical_dist, horizontal_dist))
ray_azimuth_deg = np.degrees(np.arctan2(ray_positions[:, 1], ray_positions[:, 0]))
# Assign rays to rings
ray_ring_indices = (
np.searchsorted(-ring_boundaries_deg[:-1], -ray_elevation_deg, side="right") - 1
)
ray_ring_indices = np.clip(ray_ring_indices, 0, n_rings - 1)
# Global first arrival for reference
global_first = np.min(ray_times)
# Get constant-size grid
grid = rings.get_constant_size_grid(az_bin_size_m, az_range_deg)
results = []
for bin_spec in grid:
ring_idx = bin_spec["ring_idx"]
az_lo = bin_spec["az_lo_deg"]
az_hi = bin_spec["az_hi_deg"]
# Find rays in this bin
mask = (
(ray_ring_indices == ring_idx)
& (ray_azimuth_deg >= az_lo)
& (ray_azimuth_deg < az_hi)
)
n_rays = np.sum(mask)
if n_rays >= min_rays_per_bin:
bin_times = ray_times[mask]
bin_int = ray_intensities[mask]
# First arrival relative to global first
bin_first = np.min(bin_times)
first_arrival_ns = (bin_first - global_first) * 1e9
# Time spread within bin (10-90 percentile, intensity-weighted)
bin_times_rel = (bin_times - bin_first) * 1e9 # ns from bin's first
time_spread_ns = 0.0
if np.sum(bin_int) > 0 and n_rays >= 3:
sorted_idx = np.argsort(bin_times_rel)
sorted_t = bin_times_rel[sorted_idx]
sorted_w = bin_int[sorted_idx]
cumsum = np.cumsum(sorted_w)
cumsum_norm = cumsum / cumsum[-1]
t10 = sorted_t[np.searchsorted(cumsum_norm, 0.10)]
t90 = sorted_t[np.searchsorted(cumsum_norm, 0.90)]
time_spread_ns = t90 - t10
results.append(
{
"ring_idx": ring_idx,
"az_bin_idx": bin_spec["az_bin_idx"],
"n_az_bins": bin_spec["n_az_bins"],
"az_center_deg": bin_spec["az_center_deg"],
"elev_center_deg": bin_spec["elev_center_deg"],
"distance_km": bin_spec["distance_m"] / 1000,
"n_rays": n_rays,
"total_intensity": np.sum(bin_int),
"first_arrival_ns": first_arrival_ns,
"time_spread_ns": time_spread_ns,
"bin_area_m2": bin_spec["bin_area_m2"],
}
)
return results
[docs]
def plot_constant_size_timing_analysis(
rings: ConstantSizeDetectorRings,
bin_stats: list[dict],
output_path: Optional[Union[str, Path]] = None,
figsize: tuple = (16, 10),
dpi: int = 150,
) -> Optional[plt.Figure]:
"""
Plot timing analysis for constant physical size bins.
Creates a 2x2 figure:
- (a) Time spread vs distance (each point is one bin)
- (b) First arrival vs distance
- (c) Time spread heatmap (ring vs azimuth, scaled to show structure)
- (d) Number of azimuth bins per ring (shows scaling with distance)
Parameters
----------
rings : ConstantSizeDetectorRings
Detector ring configuration
bin_stats : list of dict
Results from compute_constant_size_bin_stats
output_path : str or Path, optional
If provided, save figure to this path
figsize : tuple
Figure size in inches
dpi : int
Resolution for saved figure
Returns
-------
Figure or None
Matplotlib figure object, or None if no data
"""
if not bin_stats:
return None
fig, axes = plt.subplots(2, 2, figsize=figsize)
# Extract data
_distances = np.array([s["distance_km"] for s in bin_stats]) # noqa: F841
time_spreads = np.array([s["time_spread_ns"] for s in bin_stats])
first_arrivals = np.array([s["first_arrival_ns"] for s in bin_stats])
n_rays = np.array([s["n_rays"] for s in bin_stats])
ring_indices = np.array([s["ring_idx"] for s in bin_stats])
az_centers = np.array([s["az_center_deg"] for s in bin_stats])
# Color by log ray count
colors = np.log10(n_rays + 1)
# Get tick positions for elevation axis (use unique rings)
unique_ring_indices = np.unique(ring_indices)
tick_ring_indices = unique_ring_indices[:: max(1, len(unique_ring_indices) // 5)]
tick_elevs = [
rings.ring_centers_deg[int(i)] for i in tick_ring_indices if i < rings.n_rings
]
# (a) Time spread vs ring index
ax = axes[0, 0]
sc = ax.scatter(
ring_indices, time_spreads, c=colors, cmap="viridis", s=20, alpha=0.7
)
ax.set_xlabel("Ring Index", fontsize=11)
ax.set_ylabel("Time Spread (10-90%, ns)", fontsize=11)
ax.set_title("(a) Time Spread vs Ring", fontsize=12, fontweight="bold")
ax.grid(True, alpha=0.3)
cbar = plt.colorbar(sc, ax=ax)
cbar.set_label("log\u2081\u2080(ray count)")
# Add elevation axis on top
ax_top = ax.twiny()
ax_top.set_xlim(ax.get_xlim())
ax_top.set_xticks(tick_ring_indices[: len(tick_elevs)])
ax_top.set_xticklabels([f"{e:.1f}°" for e in tick_elevs])
ax_top.set_xlabel("Elevation", fontsize=10)
# (b) First arrival vs ring index
ax = axes[0, 1]
sc = ax.scatter(
ring_indices, first_arrivals, c=colors, cmap="viridis", s=20, alpha=0.7
)
ax.set_xlabel("Ring Index", fontsize=11)
ax.set_ylabel("First Arrival (ns from global first)", fontsize=11)
ax.set_title("(b) First Arrival Time vs Ring", fontsize=12, fontweight="bold")
ax.grid(True, alpha=0.3)
cbar = plt.colorbar(sc, ax=ax)
cbar.set_label("log\u2081\u2080(ray count)")
# Add elevation axis on top
ax_top = ax.twiny()
ax_top.set_xlim(ax.get_xlim())
ax_top.set_xticks(tick_ring_indices[: len(tick_elevs)])
ax_top.set_xticklabels([f"{e:.1f}°" for e in tick_elevs])
ax_top.set_xlabel("Elevation", fontsize=10)
# (c) Time spread as scatter (ring_idx vs az_center, colored by time spread)
ax = axes[1, 0]
valid_mask = time_spreads > 0
if np.any(valid_mask):
sc = ax.scatter(
az_centers[valid_mask],
ring_indices[valid_mask],
c=time_spreads[valid_mask],
cmap="plasma",
s=30,
alpha=0.7,
)
ax.set_xlabel("Azimuth (deg) - 0\u00b0 = beam direction", fontsize=11)
ax.set_ylabel("Ring Index", fontsize=11)
ax.set_title("(c) Time Spread by Position", fontsize=12, fontweight="bold")
ax.grid(True, alpha=0.3)
cbar = plt.colorbar(sc, ax=ax)
cbar.set_label("Time Spread (ns)")
# Secondary y-axis for elevation
def ring_to_elev(ring_idx):
idx = np.clip(np.asarray(ring_idx), 0, rings.n_rings - 1).astype(int)
return rings.ring_centers_deg[idx]
def elev_to_ring(elev):
return np.interp(
elev, rings.ring_centers_deg[::-1], np.arange(rings.n_rings)[::-1]
)
ax_elev = ax.secondary_yaxis("right", functions=(ring_to_elev, elev_to_ring))
ax_elev.set_ylabel("Elevation (deg)", fontsize=11)
# (d) Number of azimuth bins per ring
ax = axes[1, 1]
# Get unique rings and their bin counts
unique_rings = sorted(set(s["ring_idx"] for s in bin_stats))
n_az_bins_per_ring = []
ring_distances = []
for ring_idx in unique_rings:
ring_bins = [s for s in bin_stats if s["ring_idx"] == ring_idx]
if ring_bins:
n_az_bins_per_ring.append(ring_bins[0]["n_az_bins"])
ring_distances.append(ring_bins[0]["distance_km"])
ax.bar(unique_rings, n_az_bins_per_ring, color="steelblue", alpha=0.7)
ax.set_xlabel("Ring Index", fontsize=11)
ax.set_ylabel("Number of Azimuth Bins", fontsize=11)
ax.set_title(
"(d) Azimuth Bins per Ring (constant physical size)",
fontsize=12,
fontweight="bold",
)
ax.grid(True, alpha=0.3)
# Add elevation on secondary x-axis
ax2 = ax.twiny()
ax2.set_xlim(ax.get_xlim())
tick_rings = unique_rings[:: max(1, len(unique_rings) // 5)]
tick_elevs_d = [
rings.ring_centers_deg[int(r)] for r in tick_rings if r < rings.n_rings
]
ax2.set_xticks(tick_rings[: len(tick_elevs_d)])
ax2.set_xticklabels([f"{e:.1f}°" for e in tick_elevs_d])
ax2.set_xlabel("Elevation", fontsize=10)
bin_size_km = rings.detector_radial_size / 1000
az_bin_size_km = (
bin_stats[0]["bin_area_m2"] / rings.detector_radial_size / 1000
if bin_stats
else 10
)
fig.suptitle(
f"Constant Physical Size Bin Analysis\n"
f"(Radial: {bin_size_km:.0f} km, Azimuthal: ~{az_bin_size_km:.0f} km, "
f"{len(bin_stats)} bins with data)",
fontsize=13,
fontweight="bold",
)
fig.tight_layout()
if output_path:
fig.savefig(output_path, dpi=dpi, bbox_inches="tight")
return fig
[docs]
def plot_constant_size_timing_distributions(
rings: ConstantSizeDetectorRings,
bin_stats: list[dict],
ray_positions: np.ndarray,
ray_times: np.ndarray,
ray_intensities: np.ndarray,
n_top: int = 12,
output_path: Optional[Union[str, Path]] = None,
figsize_per_col: float = 4.0,
figsize_per_row: float = 3.5,
dpi: int = 150,
) -> Optional[plt.Figure]:
"""
Plot timing distributions for top bins by ray count.
Parameters
----------
rings : ConstantSizeDetectorRings
Detector ring configuration
bin_stats : list of dict
Results from compute_constant_size_bin_stats
ray_positions, ray_times, ray_intensities : ndarray
Raw ray data for histogram computation
n_top : int
Number of top bins to plot
output_path : str or Path, optional
If provided, save figure to this path
figsize_per_col, figsize_per_row : float
Figure size multipliers
dpi : int
Resolution for saved figure
Returns
-------
Figure or None
Matplotlib figure object, or None if no data
"""
if not bin_stats:
return None
# Sort by ray count and take top N
sorted_bins = sorted(bin_stats, key=lambda x: x["n_rays"], reverse=True)
top_bins = sorted_bins[:n_top]
if not top_bins:
return None
n_plots = len(top_bins)
n_cols = min(4, n_plots)
n_rows = (n_plots + n_cols - 1) // n_cols
fig, axes = plt.subplots(
n_rows, n_cols, figsize=(figsize_per_col * n_cols, figsize_per_row * n_rows)
)
if n_plots == 1:
axes = np.array([[axes]])
elif n_rows == 1:
axes = axes.reshape(1, -1)
# Precompute ray assignments
ring_boundaries_deg = rings.ring_boundaries_deg
n_rings = rings.n_rings
horizontal_dist = np.sqrt(ray_positions[:, 0] ** 2 + ray_positions[:, 1] ** 2)
vertical_dist = ray_positions[:, 2]
ray_elevation_deg = np.degrees(np.arctan2(vertical_dist, horizontal_dist))
ray_azimuth_deg = np.degrees(np.arctan2(ray_positions[:, 1], ray_positions[:, 0]))
ray_ring_indices = (
np.searchsorted(-ring_boundaries_deg[:-1], -ray_elevation_deg, side="right") - 1
)
ray_ring_indices = np.clip(ray_ring_indices, 0, n_rings - 1)
for i, bin_info in enumerate(top_bins):
row, col = divmod(i, n_cols)
ax = axes[row, col]
ring_idx = bin_info["ring_idx"]
n_az_bins = bin_info["n_az_bins"]
az_width = 20.0 / n_az_bins # Assuming ±10° range
az_center = bin_info["az_center_deg"]
az_lo = az_center - az_width / 2
az_hi = az_center + az_width / 2
# Find rays in this bin
mask = (
(ray_ring_indices == ring_idx)
& (ray_azimuth_deg >= az_lo)
& (ray_azimuth_deg < az_hi)
)
bin_times = ray_times[mask]
bin_int = ray_intensities[mask]
n_rays = len(bin_times)
if n_rays > 0 and np.sum(bin_int) > 0:
# Times relative to bin's first arrival
bin_first = np.min(bin_times)
times_rel = (bin_times - bin_first) * 1e9
# Intensity-weighted histogram
counts, edges = np.histogram(times_rel, bins=30, weights=bin_int)
centers = (edges[:-1] + edges[1:]) / 2
ax.fill_between(centers, counts, alpha=0.6, color="steelblue", step="mid")
ax.step(centers, counts, where="mid", color="steelblue", linewidth=1.5)
# Compute percentiles
sorted_idx = np.argsort(times_rel)
sorted_t = times_rel[sorted_idx]
sorted_w = bin_int[sorted_idx]
cumsum = np.cumsum(sorted_w)
cumsum_norm = cumsum / cumsum[-1]
t10 = sorted_t[np.searchsorted(cumsum_norm, 0.10)]
t90 = sorted_t[np.searchsorted(cumsum_norm, 0.90)]
spread = t90 - t10
ax.axvline(t10, color="red", linestyle="--", linewidth=1.5, alpha=0.8)
ax.axvline(t90, color="red", linestyle="--", linewidth=1.5, alpha=0.8)
ax.axvspan(t10, t90, alpha=0.15, color="red")
ax.set_title(
f"R{ring_idx} az={az_center:.1f}\u00b0\n"
f"d={bin_info['distance_km']:.0f}km n={n_rays} \u0394t={spread:.1f}ns",
fontsize=9,
)
else:
ax.set_title(f"R{ring_idx} az={az_center:.1f}\u00b0 (no data)", fontsize=9)
ax.set_xlabel("Time (ns)", fontsize=9)
ax.set_ylabel("Intensity", fontsize=9)
ax.grid(True, alpha=0.3)
ax.set_xlim(left=0)
# Hide unused axes
for i in range(n_plots, n_rows * n_cols):
row, col = divmod(i, n_cols)
axes[row, col].axis("off")
bin_size_km = rings.detector_radial_size / 1000
fig.suptitle(
f"Per-Bin Timing Distributions (top {n_top} by ray count)\n"
f"Constant physical size: {bin_size_km:.0f} km radial",
fontsize=13,
fontweight="bold",
)
fig.tight_layout()
if output_path:
fig.savefig(output_path, dpi=dpi, bbox_inches="tight")
return fig
[docs]
def compute_ring_spread_stats(
ray_positions: np.ndarray,
ray_times: np.ndarray,
ray_intensities: np.ndarray,
rings: ConstantSizeDetectorRings,
az_limit: float = 10.0,
n_az_bins: int = 40,
) -> list:
"""
Compute per-ring spread statistics (azimuthal and timing).
Parameters
----------
ray_positions : ndarray
Ray detection positions (N, 3)
ray_times : ndarray
Ray detection times (N,)
ray_intensities : ndarray
Ray intensities (N,)
rings : ConstantSizeDetectorRings
Detector ring configuration
az_limit : float
Azimuth limit in degrees (filter to +/- this range)
n_az_bins : int
Number of azimuth bins within the range
Returns
-------
list of dict
Per-ring statistics with keys: ring_idx, n_rays, az_std,
time_spread_ns, dist_km, total_intensity
"""
n_rings = rings.n_rings
ring_boundaries_deg = rings.ring_boundaries_deg
ring_distances = rings.ring_distances
# Compute elevation and azimuth
horizontal_dist = np.sqrt(ray_positions[:, 0] ** 2 + ray_positions[:, 1] ** 2)
vertical_dist = ray_positions[:, 2]
ray_elevation_deg = np.degrees(np.arctan2(vertical_dist, horizontal_dist))
ray_azimuth_deg = np.degrees(np.arctan2(ray_positions[:, 1], ray_positions[:, 0]))
# Assign rays to rings
ray_ring_indices = (
np.searchsorted(-ring_boundaries_deg[:-1], -ray_elevation_deg, side="right") - 1
)
ray_ring_indices = np.clip(ray_ring_indices, 0, n_rings - 1)
az_bin_width = 2 * az_limit / n_az_bins
ring_stats = []
for ring_idx in range(n_rings):
ring_mask = (ray_ring_indices == ring_idx) & (
np.abs(ray_azimuth_deg) <= az_limit
)
n_rays_in_ring = np.sum(ring_mask)
if n_rays_in_ring >= 3:
ring_az = ray_azimuth_deg[ring_mask]
ring_int = ray_intensities[ring_mask]
az_std = np.std(ring_az)
# Compute time spread per azimuth bin, then take median
bin_time_spreads = []
for az_bin in range(n_az_bins):
az_lo = -az_limit + az_bin * az_bin_width
az_hi = az_lo + az_bin_width
bin_mask = (
(ray_ring_indices == ring_idx)
& (ray_azimuth_deg >= az_lo)
& (ray_azimuth_deg < az_hi)
)
n_in_bin = np.sum(bin_mask)
if n_in_bin >= 3:
bin_times_abs = ray_times[bin_mask]
bin_int = ray_intensities[bin_mask]
bin_first = np.min(bin_times_abs)
bin_times_rel = (bin_times_abs - bin_first) * 1e9
if np.sum(bin_int) > 0:
sorted_idx = np.argsort(bin_times_rel)
sorted_t = bin_times_rel[sorted_idx]
sorted_w = bin_int[sorted_idx]
cumsum = np.cumsum(sorted_w)
cumsum_norm = cumsum / cumsum[-1]
t10 = sorted_t[np.searchsorted(cumsum_norm, 0.10)]
t90 = sorted_t[np.searchsorted(cumsum_norm, 0.90)]
bin_time_spreads.append(t90 - t10)
time_spread = np.median(bin_time_spreads) if bin_time_spreads else 0
ring_stats.append(
{
"ring_idx": ring_idx,
"n_rays": n_rays_in_ring,
"az_std": az_std,
"time_spread_ns": time_spread,
"dist_km": ring_distances[ring_idx] / 1000,
"total_intensity": np.sum(ring_int),
}
)
return ring_stats