Source code for lsurf.visualization.absorption_plots

# 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.

"""
Absorption Visualization

Plotting utilities for Beer-Lambert absorption analysis.
"""

from pathlib import Path
from typing import Callable

import numpy as np
import matplotlib.pyplot as plt


[docs] def plot_homogeneous_absorption( alpha: float, distance: float, simulated_intensity: float, simulated_optical_depth: float, initial_intensity: float = 1.0, save_path: str | Path | None = None, show: bool = True, ) -> plt.Figure: """ Plot Beer-Lambert decay for homogeneous absorbing medium. Parameters ---------- alpha : float Absorption coefficient in m^-1 distance : float Total propagation distance in m simulated_intensity : float Simulated final intensity simulated_optical_depth : float Simulated accumulated optical depth initial_intensity : float Initial intensity (default 1.0) save_path : str or Path, optional Path to save figure show : bool Whether to display the figure Returns ------- plt.Figure The matplotlib figure """ distances = np.linspace(0, distance, 100) analytical_intensity = initial_intensity * np.exp(-alpha * distances) analytical_tau = alpha * distances fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5)) # Intensity decay ax1.plot(distances / 1000, analytical_intensity, "b-", label="Analytical") ax1.scatter( [distance / 1000], [simulated_intensity], c="r", s=100, label=f"Simulated ({simulated_intensity:.4f})", zorder=5, ) ax1.set_xlabel("Distance (km)") ax1.set_ylabel("Intensity") ax1.set_title("Beer-Lambert Decay") ax1.legend() ax1.grid(True, alpha=0.3) # Optical depth ax2.plot(distances / 1000, analytical_tau, "b-", label="Analytical") ax2.scatter( [distance / 1000], [simulated_optical_depth], c="r", s=100, label=f"Simulated ({simulated_optical_depth:.4f})", zorder=5, ) ax2.set_xlabel("Distance (km)") ax2.set_ylabel("Optical Depth τ") ax2.set_title("Optical Depth Accumulation") ax2.legend() ax2.grid(True, alpha=0.3) plt.tight_layout() if save_path: fig.savefig(save_path, dpi=150) if show: plt.show() return fig
[docs] def plot_exponential_atmosphere_absorption( alpha_sea_level: float, scale_height: float, final_altitude: float, simulated_intensity: float, simulated_optical_depth: float, max_altitude: float = 100_000.0, save_path: str | Path | None = None, show: bool = True, ) -> plt.Figure: """ Plot absorption profiles for exponential atmosphere. Parameters ---------- alpha_sea_level : float Absorption coefficient at sea level in m^-1 scale_height : float Atmospheric scale height in m final_altitude : float Final ray altitude in m simulated_intensity : float Simulated final intensity simulated_optical_depth : float Simulated accumulated optical depth max_altitude : float Maximum altitude for plots in m save_path : str or Path, optional Path to save figure show : bool Whether to display the figure Returns ------- plt.Figure The matplotlib figure """ altitudes = np.linspace(0, max_altitude, 500) alpha_profile = alpha_sea_level * np.exp(-altitudes / scale_height) tau_profile = ( alpha_sea_level * scale_height * (1 - np.exp(-altitudes / scale_height)) ) intensity_profile = np.exp(-tau_profile) fig, axes = plt.subplots(1, 3, figsize=(15, 5)) # Absorption coefficient profile axes[0].semilogy(alpha_profile * 1e6, altitudes / 1000, "b-") axes[0].set_xlabel("α (×10⁻⁶ m⁻¹)") axes[0].set_ylabel("Altitude (km)") axes[0].set_title("Absorption Coefficient Profile") axes[0].grid(True, alpha=0.3) # Optical depth axes[1].plot(tau_profile, altitudes / 1000, "b-", label="Analytical") axes[1].scatter( [simulated_optical_depth], [final_altitude / 1000], c="r", s=100, label="Simulated", zorder=5, ) axes[1].set_xlabel("Optical Depth τ") axes[1].set_ylabel("Altitude (km)") axes[1].set_title("Accumulated Optical Depth") axes[1].legend() axes[1].grid(True, alpha=0.3) # Intensity axes[2].plot(intensity_profile, altitudes / 1000, "b-", label="Analytical") axes[2].scatter( [simulated_intensity], [final_altitude / 1000], c="r", s=100, label="Simulated", zorder=5, ) axes[2].set_xlabel("Intensity (I/I₀)") axes[2].set_ylabel("Altitude (km)") axes[2].set_title("Transmitted Intensity") axes[2].legend() axes[2].grid(True, alpha=0.3) plt.tight_layout() if save_path: fig.savefig(save_path, dpi=150) if show: plt.show() return fig
[docs] def plot_spectral_atmosphere_absorption( get_extinction: Callable[[float, float], float], get_optical_depth: Callable[[float, float, float], float], earth_radius: float, ray_wavelengths: np.ndarray, ray_positions: np.ndarray, ray_optical_depths: np.ndarray, wavelength_range: tuple[float, float] = (300e-9, 1500e-9), altitude_range: tuple[float, float] = (0, 50_000), test_altitudes: list[float] | None = None, test_wavelengths: list[float] | None = None, save_path: str | Path | None = None, show: bool = True, ) -> plt.Figure: """ Plot spectral absorption characteristics. Parameters ---------- get_extinction : callable Function(altitude, wavelength) -> extinction coefficient get_optical_depth : callable Function(z1, z2, wavelength) -> optical depth earth_radius : float Earth radius in m ray_wavelengths : ndarray Array of ray wavelengths in m ray_positions : ndarray Array of ray positions (N, 3) ray_optical_depths : ndarray Array of simulated optical depths wavelength_range : tuple (min, max) wavelength in m for plots altitude_range : tuple (min, max) altitude in m for plots test_altitudes : list, optional Altitudes for extinction vs wavelength plot test_wavelengths : list, optional Wavelengths for optical depth vs altitude plot save_path : str or Path, optional Path to save figure show : bool Whether to display the figure Returns ------- plt.Figure The matplotlib figure """ if test_altitudes is None: test_altitudes = [5000, 10000, 20000, 30000] if test_wavelengths is None: test_wavelengths = [400e-9, 550e-9, 700e-9, 1000e-9] wavelengths_plot = np.linspace(wavelength_range[0], wavelength_range[1], 100) altitudes_range = np.linspace(altitude_range[0], altitude_range[1], 100) fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5)) # Extinction coefficient vs wavelength at different altitudes for h in test_altitudes: alpha_vals = [get_extinction(h, wl) for wl in wavelengths_plot] ax1.semilogy(wavelengths_plot * 1e9, alpha_vals, label=f"{h/1000:.0f} km") ax1.set_xlabel("Wavelength (nm)") ax1.set_ylabel("Extinction coefficient (m⁻¹)") ax1.set_title("Extinction vs Wavelength") ax1.legend() ax1.grid(True, alpha=0.3) # Optical depth vs altitude for different wavelengths for wl in test_wavelengths: tau_vals = [get_optical_depth(0, h, wl) for h in altitudes_range] ax2.plot(tau_vals, altitudes_range / 1000, label=f"{wl*1e9:.0f} nm") # Add simulated points for i in range(len(ray_wavelengths)): final_h = ( np.sqrt( ray_positions[i, 0] ** 2 + ray_positions[i, 1] ** 2 + ray_positions[i, 2] ** 2 ) - earth_radius ) ax2.scatter(ray_optical_depths[i], final_h / 1000, c="r", s=50, zorder=5) ax2.set_xlabel("Optical Depth τ") ax2.set_ylabel("Altitude (km)") ax2.set_title("Optical Depth vs Altitude") ax2.legend() ax2.grid(True, alpha=0.3) plt.tight_layout() if save_path: fig.savefig(save_path, dpi=150) if show: plt.show() return fig
__all__ = [ "plot_homogeneous_absorption", "plot_exponential_atmosphere_absorption", "plot_spectral_atmosphere_absorption", ]