# 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.
"""
Polarization Visualization
Functions for plotting polarization vector components, ellipses,
and polarization-resolved reflectance analysis.
"""
from typing import TYPE_CHECKING
import matplotlib.pyplot as plt
import numpy as np
from matplotlib.figure import Figure
from numpy.typing import NDArray
if TYPE_CHECKING:
from ..utilities.recording_sphere import RecordedRays
from .common import save_figure
[docs]
def get_ray_coordinates(
recorded_rays: "RecordedRays", projection: str = "angular"
) -> tuple[np.ndarray, np.ndarray, str, str]:
"""Convert recorded rays to coordinates for binning/plotting.
Parameters
----------
recorded_rays : RecordedRays
Recorded rays at the detection sphere.
projection : str
Type of projection for binning:
- 'angular': Use elevation and azimuth angles
- 'spatial': Use X and Y positions on detection surface
Returns
-------
x_coord : ndarray
X coordinates for plotting/binning.
y_coord : ndarray
Y coordinates for plotting/binning.
xlabel : str
Label for x-axis.
ylabel : str
Label for y-axis.
"""
if projection == "angular":
# Use ray direction-based angles
directions = recorded_rays.directions
# Elevation: angle from horizontal plane (arcsin of z component)
y_coord = np.degrees(np.arcsin(directions[:, 2]))
# Azimuth: angle in X-Y plane from X axis
x_coord = np.degrees(np.arctan2(directions[:, 1], directions[:, 0]))
xlabel = "Azimuth (degrees)"
ylabel = "Ray Angle from Horizontal (degrees)"
else: # spatial
x_coord = recorded_rays.positions[:, 0]
y_coord = recorded_rays.positions[:, 1]
xlabel = "X Position (m)"
ylabel = "Y Position (m)"
return x_coord, y_coord, xlabel, ylabel
[docs]
def plot_polarization_vector_components(
recorded_rays: "RecordedRays",
bins: int = 50,
figsize: tuple[float, float] = (18, 6),
save_path: str | None = None,
vmin: float = None,
vmax: float = None,
cmap: str = "RdBu_r",
projection: str = "angular",
) -> Figure:
"""
Plot 3D polarization vector components at the detection surface.
Creates three subfigures showing the average X (horizontal), Y (vertical),
and Z (depth) components of the polarization vector binned by position.
Parameters
----------
recorded_rays : RecordedRays
Recorded rays at the detection sphere, must have polarization_vectors.
bins : int
Number of bins in each dimension for the 2D histogram.
figsize : tuple
Figure size (width, height).
save_path : str, optional
Path to save figure.
vmin : float, optional
Minimum value for colormap. Default is symmetric around 0.
vmax : float, optional
Maximum value for colormap. Default is symmetric around 0.
cmap : str
Colormap name. Default 'RdBu_r' is diverging (good for signed values).
projection : str
Type of projection for binning:
- 'angular': Use elevation and azimuth angles
- 'spatial': Use X and Y positions on detection surface
Returns
-------
Figure
Matplotlib figure with three subplots.
Notes
-----
The polarization vector represents the electric field direction at each ray.
For unpolarized light scattered from a wavy surface:
- X component: Horizontal polarization
- Y component: Vertical polarization
- Z component: Along depth/propagation direction
Each bin shows the intensity-weighted average polarization component.
"""
if recorded_rays.polarization_vectors is None:
raise ValueError(
"RecordedRays does not contain polarization_vectors. "
"Enable track_polarization_vector=True in process_surface_interaction."
)
n_rays = len(recorded_rays.positions)
if n_rays == 0:
raise ValueError("No rays to plot")
pol_vectors = recorded_rays.polarization_vectors
intensities = recorded_rays.intensities
fig, axes = plt.subplots(1, 3, figsize=figsize, constrained_layout=True)
# Get coordinates for binning
x_coord, y_coord, xlabel, ylabel = get_ray_coordinates(recorded_rays, projection)
# Component labels and data
components = [
("X (Horizontal)", pol_vectors[:, 0]),
("Y (Vertical)", pol_vectors[:, 1]),
("Z (Depth)", pol_vectors[:, 2]),
]
# Compute intensity-weighted average for each bin
for ax, (comp_name, comp_data) in zip(axes, components, strict=False):
# Create 2D histogram
weighted_sum, x_edges, y_edges = np.histogram2d(
x_coord,
y_coord,
bins=bins,
weights=comp_data * intensities,
)
intensity_sum, _, _ = np.histogram2d(
x_coord,
y_coord,
bins=bins,
weights=intensities,
)
# Compute average (avoid division by zero)
with np.errstate(divide="ignore", invalid="ignore"):
avg_component = np.where(
intensity_sum > 0, weighted_sum / intensity_sum, np.nan
)
# Set symmetric colormap limits if not specified
if vmin is None or vmax is None:
max_abs = np.nanmax(np.abs(avg_component))
if max_abs == 0:
max_abs = 1.0
_vmin = -max_abs if vmin is None else vmin
_vmax = max_abs if vmax is None else vmax
else:
_vmin, _vmax = vmin, vmax
# Plot
x_centers = 0.5 * (x_edges[:-1] + x_edges[1:])
y_centers = 0.5 * (y_edges[:-1] + y_edges[1:])
X, Y = np.meshgrid(x_centers, y_centers)
im = ax.pcolormesh(
X,
Y,
avg_component.T,
cmap=cmap,
vmin=_vmin,
vmax=_vmax,
shading="auto",
)
plt.colorbar(im, ax=ax, label=f"Mean {comp_name}")
ax.set_xlabel(xlabel, fontsize=11)
ax.set_ylabel(ylabel, fontsize=11)
ax.set_title(f"Polarization {comp_name}", fontsize=12, fontweight="bold")
ax.grid(True, alpha=0.3)
fig.suptitle(
"3D Polarization Vector Components at Detection Surface",
fontsize=14,
fontweight="bold",
)
if save_path:
save_figure(fig, save_path)
return fig
[docs]
def plot_polarization_ellipse(
recorded_rays: "RecordedRays",
bins: int = 30,
figsize: tuple[float, float] = (12, 10),
save_path: str | None = None,
projection: str = "angular",
arrow_scale: float = 1.0,
) -> Figure:
"""
Plot polarization state as arrows/ellipses on the detection surface.
Shows the average polarization direction in each bin as an arrow.
Parameters
----------
recorded_rays : RecordedRays
Recorded rays at the detection sphere, must have polarization_vectors.
bins : int
Number of bins in each dimension.
figsize : tuple
Figure size.
save_path : str, optional
Path to save figure.
projection : str
Type of projection: 'angular' or 'spatial'.
arrow_scale : float
Scale factor for arrow lengths.
Returns
-------
Figure
Matplotlib figure with polarization arrows.
"""
if recorded_rays.polarization_vectors is None:
raise ValueError("RecordedRays does not contain polarization_vectors.")
n_rays = len(recorded_rays.positions)
if n_rays == 0:
raise ValueError("No rays to plot")
pol_vectors = recorded_rays.polarization_vectors
intensities = recorded_rays.intensities
fig, ax = plt.subplots(figsize=figsize)
# Get coordinates for binning
x_coord, y_coord, xlabel, ylabel = get_ray_coordinates(recorded_rays, projection)
# Compute intensity-weighted average polarization in each bin
pol_x_sum, x_edges, y_edges = np.histogram2d(
x_coord, y_coord, bins=bins, weights=pol_vectors[:, 0] * intensities
)
pol_y_sum, _, _ = np.histogram2d(
x_coord, y_coord, bins=bins, weights=pol_vectors[:, 1] * intensities
)
intensity_sum, _, _ = np.histogram2d(
x_coord, y_coord, bins=bins, weights=intensities
)
# Compute bin centers
x_centers = 0.5 * (x_edges[:-1] + x_edges[1:])
y_centers = 0.5 * (y_edges[:-1] + y_edges[1:])
# Background: intensity distribution
im = ax.pcolormesh(
x_edges,
y_edges,
intensity_sum.T,
cmap="gray_r",
alpha=0.5,
shading="auto",
)
plt.colorbar(im, ax=ax, label="Intensity sum", shrink=0.7)
# Draw polarization arrows
for i in range(len(x_centers)):
for j in range(len(y_centers)):
if intensity_sum[i, j] > 0:
avg_pol_x = pol_x_sum[i, j] / intensity_sum[i, j]
avg_pol_y = pol_y_sum[i, j] / intensity_sum[i, j]
# Arrow length proportional to polarization magnitude
pol_mag = np.sqrt(avg_pol_x**2 + avg_pol_y**2)
if pol_mag > 0.01: # Only draw significant polarization
# Scale arrows to fit in bins
bin_size = min(x_edges[1] - x_edges[0], y_edges[1] - y_edges[0])
scale = bin_size * 0.4 * arrow_scale / max(pol_mag, 0.1)
ax.arrow(
x_centers[i],
y_centers[j],
avg_pol_x * scale,
avg_pol_y * scale,
head_width=bin_size * 0.1,
head_length=bin_size * 0.05,
fc="red",
ec="darkred",
alpha=0.8,
linewidth=0.5,
)
ax.set_xlabel(xlabel, fontsize=12)
ax.set_ylabel(ylabel, fontsize=12)
ax.set_title(
"Polarization Direction at Detection Surface\n"
"(arrows show X-Y polarization component)",
fontsize=13,
fontweight="bold",
)
ax.grid(True, alpha=0.3)
ax.set_aspect("equal")
if save_path:
save_figure(fig, save_path)
return fig
[docs]
def plot_polarization_vs_elevation(
recorded_rays: "RecordedRays",
bins: int = 50,
figsize: tuple[float, float] = (14, 5),
save_path: str | None = None,
) -> Figure:
"""
Plot polarization degree as a function of ray elevation angle.
Averages over all azimuth angles to show how polarization varies with
the ray angle from horizontal.
Parameters
----------
recorded_rays : RecordedRays
Recorded rays at the detection sphere, must have polarization_vectors.
bins : int
Number of elevation angle bins.
figsize : tuple
Figure size (width, height).
save_path : str, optional
Path to save figure.
Returns
-------
Figure
Matplotlib figure with polarization vs elevation plots.
Notes
-----
The polarization degree is computed as:
DoP = (E_x² - E_y²) / (E_x² + E_y²)
Where E_x is the horizontal component and E_y is the vertical component
of the polarization vector. DoP = +1 means fully horizontal polarization,
DoP = -1 means fully vertical polarization, DoP = 0 means unpolarized
or 45° linear polarization.
"""
if recorded_rays.polarization_vectors is None:
raise ValueError("RecordedRays does not contain polarization_vectors.")
n_rays = len(recorded_rays.positions)
if n_rays == 0:
raise ValueError("No rays to plot")
pol_vectors = recorded_rays.polarization_vectors
intensities = recorded_rays.intensities
directions = recorded_rays.directions
# Compute elevation angle from ray direction
elevation = np.degrees(np.arcsin(directions[:, 2]))
# Create figure with 3 subplots
fig, axes = plt.subplots(1, 3, figsize=figsize, constrained_layout=True)
# Bin edges
elev_min, elev_max = elevation.min(), elevation.max()
bin_edges = np.linspace(elev_min, elev_max, bins + 1)
bin_centers = 0.5 * (bin_edges[:-1] + bin_edges[1:])
# Compute intensity-weighted statistics in each bin
Ex_weighted = np.zeros(bins)
Ey_weighted = np.zeros(bins)
Ez_weighted = np.zeros(bins)
Ex2_weighted = np.zeros(bins)
Ey2_weighted = np.zeros(bins)
intensity_sum = np.zeros(bins)
ray_count = np.zeros(bins)
bin_indices = np.digitize(elevation, bin_edges) - 1
bin_indices = np.clip(bin_indices, 0, bins - 1)
for i in range(bins):
mask = bin_indices == i
if np.any(mask):
weights = intensities[mask]
total_weight = np.sum(weights)
if total_weight > 0:
Ex_weighted[i] = np.sum(pol_vectors[mask, 0] * weights) / total_weight
Ey_weighted[i] = np.sum(pol_vectors[mask, 1] * weights) / total_weight
Ez_weighted[i] = np.sum(pol_vectors[mask, 2] * weights) / total_weight
Ex2_weighted[i] = (
np.sum(pol_vectors[mask, 0] ** 2 * weights) / total_weight
)
Ey2_weighted[i] = (
np.sum(pol_vectors[mask, 1] ** 2 * weights) / total_weight
)
intensity_sum[i] = total_weight
ray_count[i] = np.sum(mask)
# Compute polarization degree: (E_x² - E_y²) / (E_x² + E_y²)
with np.errstate(divide="ignore", invalid="ignore"):
polarization_degree = np.where(
(Ex2_weighted + Ey2_weighted) > 0,
(Ex2_weighted - Ey2_weighted) / (Ex2_weighted + Ey2_weighted),
np.nan,
)
# Plot 1: Mean polarization components
ax1 = axes[0]
valid = intensity_sum > 0
ax1.plot(
bin_centers[valid],
Ex_weighted[valid],
"r-",
linewidth=2,
label="E_x (horizontal)",
)
ax1.plot(
bin_centers[valid],
Ey_weighted[valid],
"b-",
linewidth=2,
label="E_y (vertical)",
)
ax1.plot(
bin_centers[valid], Ez_weighted[valid], "g-", linewidth=2, label="E_z (depth)"
)
ax1.axhline(0, color="k", linestyle="--", alpha=0.3)
ax1.set_xlabel("Ray Angle from Horizontal (degrees)", fontsize=11)
ax1.set_ylabel("Mean Polarization Component", fontsize=11)
ax1.set_title("Polarization Components vs Elevation", fontweight="bold")
ax1.legend(loc="best")
ax1.grid(True, alpha=0.3)
# Plot 2: Polarization degree
ax2 = axes[1]
ax2.plot(bin_centers[valid], polarization_degree[valid], "purple", linewidth=2)
ax2.axhline(0, color="k", linestyle="--", alpha=0.3)
ax2.axhline(1, color="r", linestyle=":", alpha=0.5, label="Fully horizontal")
ax2.axhline(-1, color="b", linestyle=":", alpha=0.5, label="Fully vertical")
ax2.set_xlabel("Ray Angle from Horizontal (degrees)", fontsize=11)
ax2.set_ylabel(
r"Polarization Degree $(E_x^2 - E_y^2)/(E_x^2 + E_y^2)$", fontsize=11
)
ax2.set_title("Degree of Linear Polarization", fontweight="bold")
ax2.set_ylim(-1.1, 1.1)
ax2.legend(loc="best")
ax2.grid(True, alpha=0.3)
# Plot 3: Intensity distribution and ray count
ax3 = axes[2]
ax3_twin = ax3.twinx()
ax3.bar(
bin_centers,
intensity_sum,
width=(elev_max - elev_min) / bins * 0.8,
alpha=0.6,
color="orange",
label="Total intensity",
)
ax3_twin.plot(bin_centers, ray_count, "k-", linewidth=1.5, label="Ray count")
ax3.set_xlabel("Ray Angle from Horizontal (degrees)", fontsize=11)
ax3.set_ylabel("Total Intensity", fontsize=11, color="orange")
ax3_twin.set_ylabel("Ray Count", fontsize=11)
ax3.set_title("Intensity Distribution", fontweight="bold")
ax3.tick_params(axis="y", labelcolor="orange")
ax3.grid(True, alpha=0.3)
# Combined legend
lines1, labels1 = ax3.get_legend_handles_labels()
lines2, labels2 = ax3_twin.get_legend_handles_labels()
ax3.legend(lines1 + lines2, labels1 + labels2, loc="best")
fig.suptitle(
"Polarization vs Ray Elevation (averaged over azimuth)",
fontsize=14,
fontweight="bold",
)
if save_path:
save_figure(fig, save_path)
return fig
[docs]
def plot_fresnel_reflectance_curves(
n1: float = 1.0,
n2: float = 1.33,
angles_deg: NDArray[np.float64] | None = None,
brewster_angle_deg: float | None = None,
figsize: tuple[float, float] = (8, 6),
save_path: str | None = None,
) -> Figure:
"""
Plot Fresnel reflection coefficients R_s, R_p, and R_unpolarized vs elevation.
Creates a single-subplot figure showing how the reflection coefficients
for s-polarization, p-polarization, and unpolarized light vary with
elevation angle (angle above horizontal). This convention matches the
measured reflectance plots for easy comparison.
Parameters
----------
n1 : float, optional
Refractive index of incident medium (default: 1.0 for air)
n2 : float, optional
Refractive index of transmitting medium (default: 1.33 for water)
angles_deg : ndarray, optional
Elevation angles in degrees to compute (default: 0 to 90 in 0.5° steps)
brewster_angle_deg : float, optional
Brewster angle (as incidence angle) to mark on plot. If None, computed from n1, n2.
figsize : tuple, optional
Figure size in inches
save_path : str, optional
Path to save figure
Returns
-------
Figure
Matplotlib figure
Notes
-----
The x-axis shows elevation angle (angle above horizontal), which relates to
incidence angle as: elevation = 90° - incidence_angle.
- Low elevation (grazing) → high incidence angle → high reflectance
- High elevation (steep) → low incidence angle → low reflectance
Fresnel equations for reflection:
- R_s = |r_s|² where r_s = (n1*cos_i - n2*cos_t) / (n1*cos_i + n2*cos_t)
- R_p = |r_p|² where r_p = (n2*cos_i - n1*cos_t) / (n2*cos_i + n1*cos_t)
- R_unpolarized = (R_s + R_p) / 2
Brewster angle: θ_B = arctan(n2/n1), where R_p = 0
"""
from ..utilities.fresnel import fresnel_coefficients
# Elevation angles (angle above horizontal)
if angles_deg is None:
elevation_deg = np.linspace(1, 90, 179) # 1° to 90° elevation
else:
elevation_deg = angles_deg
# Convert elevation to incidence angle: incidence = 90° - elevation
incidence_deg = 90 - elevation_deg
incidence_rad = np.radians(incidence_deg)
# Compute Fresnel coefficients at each incidence angle
R_s = np.zeros(len(elevation_deg))
R_p = np.zeros(len(elevation_deg))
for i, angle_rad in enumerate(incidence_rad):
cos_theta = np.cos(angle_rad)
R_s_val, _ = fresnel_coefficients(n1, n2, cos_theta, "s")
R_p_val, _ = fresnel_coefficients(n1, n2, cos_theta, "p")
R_s[i] = float(R_s_val.item() if hasattr(R_s_val, "item") else R_s_val)
R_p[i] = float(R_p_val.item() if hasattr(R_p_val, "item") else R_p_val)
# Unpolarized is average of s and p
R_unpol = (R_s + R_p) / 2
# Compute polarization fractions: what fraction of reflected light is s vs p
R_total = R_s + R_p
with np.errstate(divide="ignore", invalid="ignore"):
F_s = np.where(R_total > 0, R_s / R_total, 0.5) # s-fraction
F_p = np.where(R_total > 0, R_p / R_total, 0.5) # p-fraction
# Compute Brewster angle if not provided (as incidence angle)
if brewster_angle_deg is None:
brewster_angle_deg = np.degrees(np.arctan(n2 / n1))
# Convert Brewster angle to elevation
brewster_elevation_deg = 90 - brewster_angle_deg
# Create figure with two subplots
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(figsize[0] * 1.8, figsize[1]))
# Left plot: Reflectance coefficients (decreasing with elevation)
ax1.plot(elevation_deg, R_s, "r-", linewidth=2, label=r"$R_s$ (s-polarization)")
ax1.plot(elevation_deg, R_p, "b-", linewidth=2, label=r"$R_p$ (p-polarization)")
ax1.plot(
elevation_deg, R_unpol, "k--", linewidth=2, label=r"$R_{unpol}$ (unpolarized)"
)
ax1.axvline(
brewster_elevation_deg,
color="green",
linestyle=":",
linewidth=1.5,
alpha=0.8,
label=f"Brewster ({brewster_elevation_deg:.1f}°)",
)
ax1.plot(brewster_elevation_deg, 0, "go", markersize=8, zorder=10)
ax1.set_xlabel("Elevation Angle (degrees)", fontsize=12)
ax1.set_ylabel("Reflectance Coefficient", fontsize=12)
ax1.set_title("Fresnel Reflectance", fontweight="bold")
ax1.set_xlim(0, 90)
ax1.set_ylim(0, 1)
ax1.legend(loc="upper right", fontsize=9)
ax1.grid(True, alpha=0.3)
# Right plot: Polarization fractions (shows what fraction of reflected light is s vs p)
ax2.plot(
elevation_deg, F_s, "r-", linewidth=2, label=r"$R_s/(R_s+R_p)$ (s-fraction)"
)
ax2.plot(
elevation_deg, F_p, "b-", linewidth=2, label=r"$R_p/(R_s+R_p)$ (p-fraction)"
)
ax2.axvline(
brewster_elevation_deg,
color="green",
linestyle=":",
linewidth=1.5,
alpha=0.8,
label=f"Brewster ({brewster_elevation_deg:.1f}°)",
)
# At Brewster, F_s = 1.0
ax2.plot(brewster_elevation_deg, 1.0, "go", markersize=8, zorder=10)
ax2.annotate(
r"$F_s = 1$",
xy=(brewster_elevation_deg, 1.0),
xytext=(brewster_elevation_deg + 5, 0.85),
fontsize=10,
arrowprops={"arrowstyle": "->", "color": "green", "alpha": 0.7},
)
ax2.set_xlabel("Elevation Angle (degrees)", fontsize=12)
ax2.set_ylabel("Polarization Fraction", fontsize=12)
ax2.set_title("Polarization of Reflected Light", fontweight="bold")
ax2.set_xlim(0, 90)
ax2.set_ylim(0, 1)
ax2.legend(loc="center right", fontsize=9)
ax2.grid(True, alpha=0.3)
fig.suptitle(
f"Fresnel Reflection: Air ($n$ = {n1:.3f}) → Water ($n$ = {n2:.3f})",
fontsize=13,
fontweight="bold",
)
plt.tight_layout()
if save_path:
save_figure(fig, save_path)
return fig
[docs]
def plot_measured_polarization_reflectance(
recorded_rays: "RecordedRays",
bins: int = 50,
figsize: tuple[float, float] = (10, 6),
save_path: str | None = None,
) -> Figure:
"""
Plot measured reflected intensity decomposed into s and p polarization.
Analyzes recorded rays to show how the reflected intensity is distributed
between s-polarization (horizontal) and p-polarization (vertical) as a
function of ray elevation angle. This treats the surface as an unknown
reflector and measures its polarization behavior.
Parameters
----------
recorded_rays : RecordedRays
Recorded rays with polarization vectors
bins : int, optional
Number of elevation angle bins (default: 50)
figsize : tuple, optional
Figure size in inches
save_path : str, optional
Path to save figure
Returns
-------
Figure
Matplotlib figure
Notes
-----
For each ray, the polarization vector is decomposed into:
- s-component: horizontal (perpendicular to vertical plane containing ray)
- p-component: vertical (in the vertical plane containing ray)
The intensity is then weighted by |E_s|² and |E_p|² to get the
intensity in each polarization state.
"""
if recorded_rays.polarization_vectors is None:
raise ValueError("Recorded rays must have polarization vectors")
directions = recorded_rays.directions
polarizations = recorded_rays.polarization_vectors
intensities = recorded_rays.intensities
# Compute elevation angle from ray direction
# Elevation = angle from horizontal plane = arcsin(z)
elevation_deg = np.degrees(np.arcsin(directions[:, 2]))
# Compute s and p basis vectors for each ray
# s = horizontal = ray × Z (perpendicular to vertical plane)
# p = vertical in plane = ray × s
z_axis = np.array([0, 0, 1], dtype=np.float32)
s_vectors = np.cross(directions, z_axis)
s_norms = np.linalg.norm(s_vectors, axis=1, keepdims=True)
# Handle rays parallel to Z (vertical rays)
parallel_mask = s_norms.squeeze() < 1e-6
if np.any(parallel_mask):
s_vectors[parallel_mask] = np.array([1, 0, 0], dtype=np.float32)
s_norms[parallel_mask] = 1.0
s_vectors = s_vectors / np.maximum(s_norms, 1e-10)
# p = ray × s (vertical component perpendicular to ray)
p_vectors = np.cross(directions, s_vectors)
p_norms = np.linalg.norm(p_vectors, axis=1, keepdims=True)
p_vectors = p_vectors / np.maximum(p_norms, 1e-10)
# Project polarization onto s and p
E_s = np.sum(polarizations * s_vectors, axis=1) # s-component magnitude
E_p = np.sum(polarizations * p_vectors, axis=1) # p-component magnitude
# Intensity in each polarization state
I_s = intensities * E_s**2 # Intensity with s-polarization character
I_p = intensities * E_p**2 # Intensity with p-polarization character
I_total = intensities
# Bin by elevation angle
elevation_min, elevation_max = elevation_deg.min(), elevation_deg.max()
bin_edges = np.linspace(elevation_min, elevation_max, bins + 1)
bin_centers = (bin_edges[:-1] + bin_edges[1:]) / 2
# Compute binned statistics
I_s_binned = np.zeros(bins)
I_p_binned = np.zeros(bins)
I_total_binned = np.zeros(bins)
ray_counts = np.zeros(bins)
bin_indices = np.digitize(elevation_deg, bin_edges) - 1
bin_indices = np.clip(bin_indices, 0, bins - 1)
for i in range(bins):
mask = bin_indices == i
if np.any(mask):
I_s_binned[i] = np.sum(I_s[mask])
I_p_binned[i] = np.sum(I_p[mask])
I_total_binned[i] = np.sum(I_total[mask])
ray_counts[i] = np.sum(mask)
# Normalize to show relative reflectance (per ray or as fraction)
# Option 1: Show absolute intensity
# Option 2: Normalize by ray count to show mean intensity per ray
# Option 3: Normalize by total to show as fraction
# Use mean intensity per ray for a cleaner comparison
valid_bins = ray_counts > 0
I_s_mean = np.zeros(bins)
I_p_mean = np.zeros(bins)
I_total_mean = np.zeros(bins)
I_s_mean[valid_bins] = I_s_binned[valid_bins] / ray_counts[valid_bins]
I_p_mean[valid_bins] = I_p_binned[valid_bins] / ray_counts[valid_bins]
I_total_mean[valid_bins] = I_total_binned[valid_bins] / ray_counts[valid_bins]
# Compute polarization fractions for each bin
I_sp_sum = I_s_mean + I_p_mean
F_s_measured = np.zeros(bins)
F_p_measured = np.zeros(bins)
valid_fraction = (valid_bins) & (I_sp_sum > 0)
F_s_measured[valid_fraction] = I_s_mean[valid_fraction] / I_sp_sum[valid_fraction]
F_p_measured[valid_fraction] = I_p_mean[valid_fraction] / I_sp_sum[valid_fraction]
# Create figure with two subplots
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(figsize[0] * 1.8, figsize[1]))
# Left plot: Mean intensity per ray
ax1.plot(
bin_centers[valid_bins],
I_s_mean[valid_bins],
"r-",
linewidth=2,
label=r"$I_s$ (s-polarization)",
marker="o",
markersize=3,
)
ax1.plot(
bin_centers[valid_bins],
I_p_mean[valid_bins],
"b-",
linewidth=2,
label=r"$I_p$ (p-polarization)",
marker="s",
markersize=3,
)
ax1.plot(
bin_centers[valid_bins],
I_total_mean[valid_bins],
"k--",
linewidth=2,
label=r"$I_{total}$",
alpha=0.7,
)
ax1.set_xlabel("Ray Elevation Angle (degrees)", fontsize=12)
ax1.set_ylabel("Mean Reflected Intensity", fontsize=12)
ax1.set_title("Mean Intensity per Ray", fontweight="bold")
ax1.legend(loc="best", fontsize=9)
ax1.grid(True, alpha=0.3)
# Right plot: Polarization fractions (comparable to theoretical Fresnel)
ax2.plot(
bin_centers[valid_fraction],
F_s_measured[valid_fraction],
"r-",
linewidth=2,
label=r"$I_s/(I_s+I_p)$ (s-fraction)",
marker="o",
markersize=3,
)
ax2.plot(
bin_centers[valid_fraction],
F_p_measured[valid_fraction],
"b-",
linewidth=2,
label=r"$I_p/(I_s+I_p)$ (p-fraction)",
marker="s",
markersize=3,
)
ax2.set_xlabel("Ray Elevation Angle (degrees)", fontsize=12)
ax2.set_ylabel("Polarization Fraction", fontsize=12)
ax2.set_title("Polarization of Reflected Light", fontweight="bold")
ax2.set_ylim(0, 1)
ax2.legend(loc="best", fontsize=9)
ax2.grid(True, alpha=0.3)
fig.suptitle(
"Measured Polarization-Resolved Reflectance (from recorded rays)",
fontsize=13,
fontweight="bold",
)
# Add secondary axis showing ray count on left plot
ax1_twin = ax1.twinx()
ax1_twin.fill_between(
bin_centers,
0,
ray_counts,
alpha=0.15,
color="gray",
label="Ray count",
)
ax1_twin.set_ylabel("Ray Count per Bin", fontsize=10, color="gray")
ax1_twin.tick_params(axis="y", labelcolor="gray")
plt.tight_layout()
if save_path:
save_figure(fig, save_path)
return fig