From 2cbf762b8636f2957c942de3783f3bfc47e2c681 Mon Sep 17 00:00:00 2001 From: kopecn Date: Mon, 1 Dec 2025 14:03:23 -0800 Subject: [PATCH 01/22] humph... --- .env | 2 + .gitignore | 157 ++++++ .pylintrc | 30 ++ src/pyMathTools/generators/__init__.py | 0 .../generators/sphericalGenerators.py | 115 +++++ src/pyMathTools/hints.py | 14 + src/pyMathTools/transforms/__init__.py | 0 .../transforms/sphericalTransforms.py | 81 ++++ .../plotUnitSphericalSmallCircles.py | 458 +++++++++--------- 9 files changed, 630 insertions(+), 227 deletions(-) create mode 100644 .env create mode 100644 .gitignore create mode 100644 .pylintrc create mode 100644 src/pyMathTools/generators/__init__.py create mode 100644 src/pyMathTools/generators/sphericalGenerators.py create mode 100644 src/pyMathTools/hints.py create mode 100644 src/pyMathTools/transforms/__init__.py create mode 100644 src/pyMathTools/transforms/sphericalTransforms.py diff --git a/.env b/.env new file mode 100644 index 0000000..5567480 --- /dev/null +++ b/.env @@ -0,0 +1,2 @@ +PYTHON=python3 +VENV=/tmp/pipTest diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..d609a52 --- /dev/null +++ b/.gitignore @@ -0,0 +1,157 @@ +# Byte-compiled / optimized / DLL files +__pycache__/ +*.py[cod] +*$py.class + +# Cython and C extensions +*.c +*.so +*.pyd +*.dll + +# Distribution / packaging +.Python +build/ +develop-eggs/ +dist/ +downloads/ +eggs/ +.eggs/ +lib/ +lib64/ +parts/ +sdist/ +var/ +*.egg-info/ +.installed.cfg +*.egg +MANIFEST + +# PyInstaller +# Usually contains a dist/ folder and build/ artifacts +*.manifest +*.spec + +# PyBuilder +target/ + +# Virtual environments +.env/ +.venv/ +env/ +venv/ +ENV/ +env.bak/ +venv.bak/ + +# pipenv +Pipfile.lock + +# poetry +poetry.lock + +# hatch +*.hatch/ + +# pdm +__pypackages__/ + +# tox and nox +.tox/ +.nox/ + +# pytest +.cache/ +.pytest_cache/ + +# mypy +.mypy_cache/ +.dmypy.json +dmypy.json + +# coverage and test reports +.coverage +.coverage.* +nosetests.xml +coverage.xml +*.cover +.hypothesis/ + +# Pyre type checker +.pyre/ + +# CMake files +CMakeFiles/ +CMakeCache.txt +cmake_install.cmake + + +# Translations +*.mo +*.pot + +# Django stuff +*.log +local_settings.py +db.sqlite3 + +# Flask stuff +instance/ +.webassets-cache + +# Scrapy stuff +.scrapy + +# IPython / Jupyter +.ipynb_checkpoints/ +profile + +# macOS system files +*.DS_Store +.AppleDouble +.LSOverride + +# Thumbnails and metadata +._* + +# Files that might appear in root of a volume +.DocumentRevisions-V100 +.fseventsd +.Spotlight-V100 +.TemporaryItems +.Trashes +.VolumeIcon.icns +.com.apple.timemachine.donotpresent + +# Directories created on remote AFP share +.AppleDB +.AppleDesktop +Network Trash Folder +Temporary Items +.apdisk + +# Installer logs +pip-log.txt +pip-delete-this-directory.txt + +# Unit test / coverage reports +htmlcov/ + +# Sphinx documentation build +docs/_build/ + +# IDEs and editors +# IntelliJ IDEA family +.idea/ +*.iml +*.ipr +*.iws +.idea_modules/ + +# Visual Studio Code +.vscode/ + +# Cookiecutter +output/ +python_boilerplate/ +cookiecutter-pypackage-env/ diff --git a/.pylintrc b/.pylintrc new file mode 100644 index 0000000..2dc5dc7 --- /dev/null +++ b/.pylintrc @@ -0,0 +1,30 @@ +# https://pylint.readthedocs.io/en/v2.12.2/user_guide/options.html + +[BASIC] +module-naming-style=any +const-naming-style=any +class-naming-style=any +function-naming-style=any +method-naming-style=any +attr-naming-style=any +argument-naming-style=any +variable-naming-style=any +class-attribute-naming-style=any +class-const-naming-style=any +inlinevar-naming-style=any + +# Custom regexes to allow both snake_case and camelCase +module-rgx=(?:(?P[a-z_][a-z0-9_]+)|(?P[a-z]+(?:[A-Z][a-z0-9]*)*))$ +const-rgx=(?:(?P[A-Z_][A-Z0-9_]*)|(?P[A-Z][a-zA-Z0-9]*))$ +class-rgx=(?:(?P[A-Z][a-zA-Z0-9]*)|(?P[a-z_][a-z0-9_]+))$ +function-rgx=(?:(?P[a-z_][a-z0-9_]*)|(?P[a-z]+(?:[A-Z][a-z0-9]*)*))$ +method-rgx=(?:(?P[a-z_][a-z0-9_]*)|(?P[a-z]+(?:[A-Z][a-z0-9]*)*))$ +attr-rgx=(?:(?P[a-z_][a-z0-9_]*)|(?P[a-z]+(?:[A-Z][a-z0-9]*)*))$ +argument-rgx=(?:(?P[a-z_][a-z0-9_]*)|(?P[a-z]+(?:[A-Z][a-z0-9]*)*))$ +variable-rgx=(?:(?P[a-z_][a-z0-9_]*)|(?P[a-z]+(?:[A-Z][a-z0-9]*)*))$ +class-attribute-rgx=(?:(?P[a-z_][a-z0-9_]*)|(?P[a-z]+(?:[A-Z][a-z0-9]*)*))$ +class-const-rgx=(?:(?P[A-Z_][A-Z0-9_]*)|(?P[A-Z][a-zA-Z0-9]*))$ +inlinevar-rgx=(?:(?P[a-z_][a-z0-9_]*)|(?P[a-z]+(?:[A-Z][a-z0-9]*)*))$ + +[design] +max-locals = 24 diff --git a/src/pyMathTools/generators/__init__.py b/src/pyMathTools/generators/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/src/pyMathTools/generators/sphericalGenerators.py b/src/pyMathTools/generators/sphericalGenerators.py new file mode 100644 index 0000000..3ecbe52 --- /dev/null +++ b/src/pyMathTools/generators/sphericalGenerators.py @@ -0,0 +1,115 @@ +from typing import Tuple +from numpy import ( + clip, + pi, + arctan2, + cos, + sin, + array, + cross, + linspace, + zeros, + sqrt, + arcsin, + dot, +) +from numpy.linalg import norm + + +from pyMathTools.hints import FloatOrNDArray, FloatNDArray + + +def generate_spherical_small_circle_points( + azimuth_center: FloatOrNDArray, + polar_center: FloatOrNDArray, + radius_angle: FloatOrNDArray, + num_points: int = 120, +) -> Tuple[FloatNDArray, FloatNDArray]: + """ + Generate points on a true small circle on a unit sphere. + + A small circle is the set of all points at a fixed angular distance (radius_angle) + from a center point on the sphere. + + Parameters: + ----------- + azimuth_center : float or NDArray + Azimuth angle of circle center (radians, 0 to 2π) + polar_center : float or NDArray + Polar angle of circle center (radians, -π/2 to π/2) + radius_angle : float or NDArray + Angular radius of the circle (radians) + num_points : int + Number of points to generate around the circle + + Returns: + -------- + azimuth_points : NDArray + Azimuth angles of points on the circle + polar_points : NDArray + Polar angles of points on the circle + """ + # Convert center to Cartesian coordinates + x_center = cos(polar_center) * cos(azimuth_center) + y_center = cos(polar_center) * sin(azimuth_center) + z_center = sin(polar_center) + center = array([x_center, y_center, z_center]) + + # Create orthonormal basis at the center point + # North pole direction + north = array([0, 0, 1]) + + # Tangent vector 1 (east direction) + if abs(z_center) < 0.999: # Not at poles + east = cross(north, center) + east = east / norm(east) + else: # At poles, use arbitrary east direction + east = array([1, 0, 0]) + + # Tangent vector 2 (north direction in tangent plane) + north_tangent = cross(center, east) + north_tangent = north_tangent / norm(north_tangent) + + # Generate points around the circle + theta = linspace(0, 2 * pi, num_points) + + # Points on the small circle using rotation + # The small circle is at angular distance radius_angle from center + azimuth_points = zeros(num_points) + polar_points = zeros(num_points) + + for i, t in enumerate(theta): + # Local coordinates: move radius_angle away from center in direction t + local_direction = cos(t) * east + sin(t) * north_tangent + + # Rotate center point toward local_direction by radius_angle + # Using Rodrigues' rotation formula + axis = cross(center, local_direction) + axis_norm = norm(axis) + + if axis_norm > 1e-10: + axis = axis / axis_norm + # Rotate center around axis by radius_angle + cos_angle = cos(radius_angle) + sin_angle = sin(radius_angle) + + point = ( + cos_angle * center + + sin_angle * cross(axis, center) + + (1 - cos_angle) * dot(axis, center) * axis + ) + else: + # If no rotation needed (shouldn't happen) + point = center + + # Convert back to spherical coordinates + x, y, z = point + r = sqrt(x**2 + y**2 + z**2) + x, y, z = x / r, y / r, z / r # Normalize to unit sphere + + polar_points[i] = arcsin(clip(z, -1, 1)) + azimuth_points[i] = arctan2(y, x) + if azimuth_points[i] < 0: + azimuth_points[i] += 2 * pi + + return azimuth_points, polar_points diff --git a/src/pyMathTools/hints.py b/src/pyMathTools/hints.py new file mode 100644 index 0000000..1033936 --- /dev/null +++ b/src/pyMathTools/hints.py @@ -0,0 +1,14 @@ +""" +Type hinting utilities for complex mathematical operations. + +This module provides reusable type hints for NumPy arrays and scalars, +enabling strict and clear typing in mathematical code throughout the package. +""" + +from typing import Union +from numpy import float64 +from numpy.typing import NDArray + +FloatOrNDArray = Union[float, NDArray[float64]] + +FloatNDArray = NDArray[float64] diff --git a/src/pyMathTools/transforms/__init__.py b/src/pyMathTools/transforms/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/src/pyMathTools/transforms/sphericalTransforms.py b/src/pyMathTools/transforms/sphericalTransforms.py new file mode 100644 index 0000000..79b7b1b --- /dev/null +++ b/src/pyMathTools/transforms/sphericalTransforms.py @@ -0,0 +1,81 @@ +""" +Transforms for spherical coordinate systems. + +This module provides functions for converting to and from spherical coordinates, +including projections and other mathematical transformations commonly used in +geometry, mapping, and scientific computing. +""" + +from numpy import log, clip, pi, tan, arctan2, cos, sin, arctanh +from pyMathTools.hints import FloatOrNDArray + + +def mercator_transform( + azimuth: FloatOrNDArray, polar: FloatOrNDArray +) -> tuple[FloatOrNDArray, FloatOrNDArray]: + """ + Apply Mercator projection to spherical coordinates. + + Parameters: + ----------- + azimuth : float or ndarray + Azimuth angle(s) in radians (0 to 2π) + polar : float or ndarray + Polar angle(s) in radians (-π/2 to π/2) + + Returns: + -------- + - x : float or ndarray + Mercator x coordinate (longitude) + - y : float or ndarray + Polar (radians) coordinate (latitude projection) + """ + x = azimuth + # Mercator projection: y = ln(tan(π/4 + lat/2)) + # Clamp polar to avoid infinities at poles + polar_clamped = clip(polar, -pi / 2 + 0.001, pi / 2 - 0.001) + y = log(tan(pi / 4 + polar_clamped / 2)) + return x, y + + +def gauss_kruger_transform( + azimuth: FloatOrNDArray, + polar: FloatOrNDArray, + central_meridian: float = pi, +) -> tuple[FloatOrNDArray, FloatOrNDArray]: + """ + Apply Gauss-Krüger (Transverse Mercator) projection to spherical coordinates. + + Parameters: + ----------- + azimuth : float or ndarray + Azimuth angle(s) in radians (0 to 2π) + polar : float or ndarray + Polar angle(s) in radians (-π/2 to π/2) + central_meridian : float, optional + Central meridian in radians (default: π) + + Returns: + -------- + easting : float or ndarray + Gauss-Krüger easting coordinate + northing : float or ndarray + Gauss-Krüger northing coordinate + """ + # Adjust azimuth relative to central meridian + lambda_rel = azimuth - central_meridian + + # Normalize to [-π, π] + lambda_rel = arctan2(sin(lambda_rel), cos(lambda_rel)) + + # Clamp polar to avoid numerical issues + polar_clamped = clip(polar, -pi / 2 + 0.001, pi / 2 - 0.001) + + # Transverse Mercator formulas for a sphere + # Easting + easting = arctanh(cos(polar_clamped) * sin(lambda_rel)) + + # Northing + northing = arctan2(sin(polar_clamped), cos(polar_clamped) * cos(lambda_rel)) + + return easting, northing diff --git a/src/pyMathToolsPlotHelpers/plotUnitSphericalSmallCircles.py b/src/pyMathToolsPlotHelpers/plotUnitSphericalSmallCircles.py index 7213cb5..47be1ce 100644 --- a/src/pyMathToolsPlotHelpers/plotUnitSphericalSmallCircles.py +++ b/src/pyMathToolsPlotHelpers/plotUnitSphericalSmallCircles.py @@ -1,131 +1,27 @@ -from typing import List, Union +from typing import List, Union, Tuple, Optional import numpy as np import matplotlib.pyplot as plt -from mpl_toolkits.mplot3d import Axes3D -from foundationTypes.mathTypes.UnitSphericalSmallCircle import UnitSphericalSmallCircle - - -def generate_spherical_small_circle_points(azimuth_center, polar_center, radius_angle, num_points=120): - """ - Generate points on a true small circle on a unit sphere. - - A small circle is the set of all points at a fixed angular distance (radius_angle) - from a center point on the sphere. - - Parameters: - ----------- - azimuth_center : float - Azimuth angle of circle center (radians, 0 to 2π) - polar_center : float - Polar angle of circle center (radians, -π/2 to π/2) - radius_angle : float - Angular radius of the circle (radians) - num_points : int - Number of points to generate around the circle +from matplotlib.figure import Figure +from matplotlib.axes import Axes - Returns: - -------- - azimuth_points : ndarray - Azimuth angles of points on the circle - polar_points : ndarray - Polar angles of points on the circle - """ - # Convert center to Cartesian coordinates - x_center = np.cos(polar_center) * np.cos(azimuth_center) - y_center = np.cos(polar_center) * np.sin(azimuth_center) - z_center = np.sin(polar_center) - center = np.array([x_center, y_center, z_center]) - - # Create orthonormal basis at the center point - # North pole direction - north = np.array([0, 0, 1]) - - # Tangent vector 1 (east direction) - if np.abs(z_center) < 0.999: # Not at poles - east = np.cross(north, center) - east = east / np.linalg.norm(east) - else: # At poles, use arbitrary east direction - east = np.array([1, 0, 0]) - - # Tangent vector 2 (north direction in tangent plane) - north_tangent = np.cross(center, east) - north_tangent = north_tangent / np.linalg.norm(north_tangent) - - # Generate points around the circle - theta = np.linspace(0, 2 * np.pi, num_points) - - # Points on the small circle using rotation - # The small circle is at angular distance radius_angle from center - azimuth_points = np.zeros(num_points) - polar_points = np.zeros(num_points) - - for i, t in enumerate(theta): - # Local coordinates: move radius_angle away from center in direction t - local_direction = np.cos(t) * east + np.sin(t) * north_tangent - - # Rotate center point toward local_direction by radius_angle - # Using Rodrigues' rotation formula - axis = np.cross(center, local_direction) - axis_norm = np.linalg.norm(axis) - - if axis_norm > 1e-10: - axis = axis / axis_norm - # Rotate center around axis by radius_angle - cos_angle = np.cos(radius_angle) - sin_angle = np.sin(radius_angle) - - point = (cos_angle * center + - sin_angle * np.cross(axis, center) + - (1 - cos_angle) * np.dot(axis, center) * axis) - else: - # If no rotation needed (shouldn't happen) - point = center - - # Convert back to spherical coordinates - x, y, z = point - r = np.sqrt(x**2 + y**2 + z**2) - x, y, z = x/r, y/r, z/r # Normalize to unit sphere - - polar_points[i] = np.arcsin(np.clip(z, -1, 1)) - azimuth_points[i] = np.arctan2(y, x) - if azimuth_points[i] < 0: - azimuth_points[i] += 2 * np.pi - - return azimuth_points, polar_points - - -def mercator_projection(azimuth, polar): - """ - Apply Mercator projection to spherical coordinates. - - Parameters: - ----------- - azimuth : float or ndarray - Azimuth angle(s) in radians (0 to 2π) - polar : float or ndarray - Polar angle(s) in radians (-π/2 to π/2) +from foundationTypes.mathTypes.UnitSphericalSmallCircle import UnitSphericalSmallCircle - Returns: - -------- - x : float or ndarray - Mercator x coordinate (longitude) - y : float or ndarray - Mercator y coordinate (latitude projection) - """ - x = azimuth - # Mercator projection: y = ln(tan(π/4 + lat/2)) - # Clamp polar to avoid infinities at poles - polar_clamped = np.clip(polar, -np.pi/2 + 0.001, np.pi/2 - 0.001) - y = np.log(np.tan(np.pi/4 + polar_clamped/2)) - return x, y +from pyMathTools.transforms.sphericalTransforms import ( + mercator_transform, + gauss_kruger_transform, +) +from pyMathTools.generators.sphericalGenerators import ( + generate_spherical_small_circle_points, +) def plot_spherical_small_circles( circles: Union[List[UnitSphericalSmallCircle], List[List[float]]], - num_points=120, - figsize=(10, 6), - title="Small Circles on Sphere (Mercator Projection)", -): + num_points: int = 120, + figsize: Tuple[int, int] = (10, 6), + title: str = "Small Circles on Sphere (Mercator Projection)", + show_legend: bool = True, +) -> Tuple[Figure, Axes]: """ Plot small circles on a sphere using Mercator projection. @@ -140,6 +36,8 @@ def plot_spherical_small_circles( Figure size (default: (10, 6)) title : str, optional Plot title + show_legend : bool, optional + Whether to show the legend (default: True) Returns: -------- @@ -166,7 +64,7 @@ def plot_spherical_small_circles( ) # Apply Mercator projection - x_points, y_points = mercator_projection(azimuth_points, polar_points) + x_points, y_points = mercator_transform(azimuth_points, polar_points) # Convert radius to degrees for legend radius_deg = np.degrees(radius_angle) @@ -176,8 +74,8 @@ def plot_spherical_small_circles( ax.scatter(x_points, y_points, s=10, alpha=0.6, label=label) # Set up axes - ax.set_xlabel("Longitude (radians)", fontsize=12) - ax.set_ylabel("Mercator Y", fontsize=12) + ax.set_xlabel("Azmuth (radians)", fontsize=12) + ax.set_ylabel("Polar (radians)", fontsize=12) ax.set_xlim(0, 2 * np.pi) # Add grid @@ -189,8 +87,9 @@ def plot_spherical_small_circles( ax.set_title(title, fontsize=14) - # Add legend - ax.legend() + # Add legend if requested + if show_legend: + ax.legend() plt.tight_layout() @@ -199,13 +98,13 @@ def plot_spherical_small_circles( def plot_spherical_small_circles_advanced( circles: Union[List[UnitSphericalSmallCircle], List[List[float]]], - num_points=120, - figsize=(10, 6), - title="Small Circles on Sphere (Flattened)", - colors=None, - labels=None, - show_centers=True, -): + num_points: int = 120, + figsize: Tuple[int, int] = (10, 6), + title: str = "Small Circles on Sphere (Flattened)", + colors: Optional[List[str]] = None, + labels: Optional[List[str]] = None, + show_centers: bool = True, +) -> Tuple[Figure, Axes]: """ Advanced version with more customization options. @@ -249,7 +148,7 @@ def plot_spherical_small_circles_advanced( ) # Apply Mercator projection - x_points, y_points = mercator_projection(azimuth_points, polar_points) + x_points, y_points = mercator_transform(azimuth_points, polar_points) # Determine color and label color = colors[idx] if colors is not None and idx < len(colors) else None @@ -260,7 +159,7 @@ def plot_spherical_small_circles_advanced( # Plot center point if requested if show_centers: - x_center, y_center = mercator_projection(azimuth_center, polar_center) + x_center, y_center = mercator_transform(azimuth_center, polar_center) ax.scatter( x_center, y_center, @@ -271,8 +170,8 @@ def plot_spherical_small_circles_advanced( ) # Set up axes - ax.set_xlabel("Longitude (radians)", fontsize=12) - ax.set_ylabel("Mercator Y", fontsize=12) + ax.set_xlabel("Azmuth (radians)", fontsize=12) + ax.set_ylabel("Polar (radians)", fontsize=12) ax.set_xlim(0, 2 * np.pi) # Add grid @@ -292,12 +191,89 @@ def plot_spherical_small_circles_advanced( return fig, ax +def plot_spherical_small_circles_gauss_kruger( + circles: Union[List[UnitSphericalSmallCircle], List[List[float]]], + num_points: int = 120, + figsize: Tuple[int, int] = (10, 6), + title: str = "Small Circles on Sphere (Gauss-Krüger Projection)", + central_meridian: float = np.pi, + show_legend: bool = True, +) -> Tuple[Figure, Axes]: + """ + Plot small circles on a sphere using Gauss-Krüger (Transverse Mercator) projection. + + Parameters: + ----------- + circles : List[UnitSphericalSmallCircle] or array-like + List of UnitSphericalSmallCircle objects or array of [azimuth, polar, radius] lists + num_points : int, optional + Number of points to generate around each circle (default: 120) + figsize : tuple, optional + Figure size (default: (10, 6)) + title : str, optional + Plot title + central_meridian : float, optional + Central meridian in radians (default: π) + show_legend : bool, optional + Whether to show the legend (default: True) + + Returns: + -------- + fig, ax : matplotlib figure and axes objects + """ + fig, ax = plt.subplots(figsize=figsize) + + # Convert UnitSphericalSmallCircle objects to arrays if needed + if circles and isinstance(circles[0], UnitSphericalSmallCircle): + circles = [[c.azimuth, c.polar, c.radius_angle] for c in circles] + + circles = np.array(circles) + if circles.ndim == 1: + circles = circles.reshape(1, -1) + + for idx, circle in enumerate(circles): + azimuth_center, polar_center, radius_angle = circle + + # Generate true spherical small circle points + azimuth_points, polar_points = generate_spherical_small_circle_points( + azimuth_center, polar_center, radius_angle, num_points + ) + + # Apply Gauss-Krüger projection + easting, northing = gauss_kruger_transform( + azimuth_points, polar_points, central_meridian + ) + + # Convert radius to degrees for legend + radius_deg = np.degrees(radius_angle) + label = f"Circle {idx + 1}: {radius_deg:.1f}°" + + # Plot with swapped coordinates to view down at equator + # x-axis (azimuth-related) = northing + # y-axis (polar-related) = easting + ax.scatter(northing, easting, s=10, alpha=0.6, label=label) + + ax.set_xlabel("Azimuth-related (radians)", fontsize=12) + ax.set_ylabel("Polar-related (radians)", fontsize=12) + ax.grid(True, alpha=0.3) + ax.set_title(title, fontsize=14) + ax.set_aspect("equal", adjustable="box") + + # Add legend if requested + if show_legend: + ax.legend() + + plt.tight_layout() + + return fig, ax + + def plot_spherical_small_circles_polar( circles: Union[List[UnitSphericalSmallCircle], List[List[float]]], - num_points=120, - figsize=(8, 8), - title="Small Circles on Sphere (Polar View)", -): + num_points: int = 120, + figsize: Tuple[int, int] = (8, 8), + title: str = "Small Circles on Sphere (Polar View)", +) -> Tuple[Figure, Axes]: """ Plot small circles on a sphere using polar projection. Azimuth is shown as the angle, polar angle is mapped to radius. @@ -317,7 +293,7 @@ def plot_spherical_small_circles_polar( -------- fig, ax : matplotlib figure and axes objects """ - fig, ax = plt.subplots(figsize=figsize, subplot_kw={'projection': 'polar'}) + fig, ax = plt.subplots(figsize=figsize, subplot_kw={"projection": "polar"}) # Convert UnitSphericalSmallCircle objects to arrays if needed if circles and isinstance(circles[0], UnitSphericalSmallCircle): @@ -330,20 +306,12 @@ def plot_spherical_small_circles_polar( for idx, circle in enumerate(circles): azimuth_center, polar_center, radius_angle = circle - # Generate angles around the circle - theta = np.linspace(0, 2 * np.pi, num_points) - - # Generate circle points in spherical space - azimuth_points = azimuth_center + radius_angle * np.cos(theta) - polar_points = polar_center + radius_angle * np.sin(theta) - - # Handle wrapping for azimuth (0 to 2*pi) - azimuth_points = np.mod(azimuth_points, 2 * np.pi) - - # Clip polar angle to valid range (-pi/2 to pi/2) - polar_points = np.clip(polar_points, -np.pi / 2, np.pi / 2) + # Generate true spherical small circle points + azimuth_points, polar_points = generate_spherical_small_circle_points( + azimuth_center, polar_center, radius_angle, num_points + ) - # Map polar angle to radius (0 to 1) + # Map polar angle to radius (0 to 1) for polar plot # polar ranges from -pi/2 to pi/2, map to 0 to 1 radius_points = (polar_points + np.pi / 2) / np.pi @@ -355,7 +323,7 @@ def plot_spherical_small_circles_polar( ax.scatter(azimuth_points, radius_points, s=10, alpha=0.6, label=label) ax.set_title(title, fontsize=14, pad=20) - ax.legend(loc='upper right', bbox_to_anchor=(1.3, 1.1)) + ax.legend(loc="upper right", bbox_to_anchor=(1.3, 1.1)) ax.grid(True, alpha=0.3) plt.tight_layout() @@ -365,11 +333,12 @@ def plot_spherical_small_circles_polar( def plot_spherical_small_circles_3d( circles: Union[List[UnitSphericalSmallCircle], List[List[float]]], - num_points=120, - figsize=(10, 10), - title="Small Circles on Unit Sphere (3D)", - show_sphere=True, -): + num_points: int = 120, + figsize: Tuple[int, int] = (10, 10), + title: str = "Small Circles on Unit Sphere (3D)", + show_sphere: bool = True, + show_legend: bool = True, +) -> Tuple[Figure, Axes]: """ Plot small circles on a 3D unit sphere. @@ -385,13 +354,15 @@ def plot_spherical_small_circles_3d( Plot title show_sphere : bool, optional Whether to show the sphere surface (default: True) + show_legend : bool, optional + Whether to show the legend (default: True) Returns: -------- fig, ax : matplotlib figure and axes objects """ fig = plt.figure(figsize=figsize) - ax = fig.add_subplot(111, projection='3d') + ax = fig.add_subplot(111, projection="3d") # Convert UnitSphericalSmallCircle objects to arrays if needed if circles and isinstance(circles[0], UnitSphericalSmallCircle): @@ -408,22 +379,16 @@ def plot_spherical_small_circles_3d( x_sphere = np.outer(np.cos(v), np.cos(u)) y_sphere = np.outer(np.cos(v), np.sin(u)) z_sphere = np.outer(np.sin(v), np.ones(np.size(u))) - ax.plot_surface(x_sphere, y_sphere, z_sphere, alpha=0.1, color='lightblue') + ax.plot_surface(x_sphere, y_sphere, z_sphere, alpha=0.1, color="lightblue") # Plot each circle for idx, circle in enumerate(circles): azimuth_center, polar_center, radius_angle = circle - # Generate points around the small circle on the sphere - theta = np.linspace(0, 2 * np.pi, num_points) - - # Generate circle points in spherical coordinates - azimuth_points = azimuth_center + radius_angle * np.cos(theta) - polar_points = polar_center + radius_angle * np.sin(theta) - - # Handle wrapping and clipping - azimuth_points = np.mod(azimuth_points, 2 * np.pi) - polar_points = np.clip(polar_points, -np.pi / 2, np.pi / 2) + # Generate true spherical small circle points + azimuth_points, polar_points = generate_spherical_small_circle_points( + azimuth_center, polar_center, radius_angle, num_points + ) # Convert to Cartesian coordinates (unit sphere, r=1) x = np.cos(polar_points) * np.cos(azimuth_points) @@ -437,11 +402,14 @@ def plot_spherical_small_circles_3d( # Plot the circle ax.scatter(x, y, z, s=10, alpha=0.8, label=label) - ax.set_xlabel('X') - ax.set_ylabel('Y') - ax.set_zlabel('Z') + ax.set_xlabel("X") + ax.set_ylabel("Y") + ax.set_zlabel("Z") ax.set_title(title, fontsize=14) - ax.legend() + + # Add legend if requested + if show_legend: + ax.legend() # Set equal aspect ratio ax.set_box_aspect([1, 1, 1]) @@ -453,12 +421,12 @@ def plot_spherical_small_circles_3d( def plot_spherical_small_circles_multiplot( circles: Union[List[UnitSphericalSmallCircle], List[List[float]]], - num_points=120, - figsize=(18, 6), - title="Small Circles on Unit Sphere - Multiple Views", -): + num_points: int = 120, + figsize: Tuple[int, int] = (18, 6), + title: str = "Small Circles on Unit Sphere - Multiple Views", +) -> Figure: """ - Create a multiplot showing flattened, polar, and 3D views of spherical small circles. + Create a multiplot showing Mercator, Gauss-Krüger, and 3D views of spherical small circles. Parameters: ----------- @@ -487,47 +455,50 @@ def plot_spherical_small_circles_multiplot( if circles_array.ndim == 1: circles_array = circles_array.reshape(1, -1) - # 1. Flattened view + # 1. Mercator projection view ax1 = fig.add_subplot(131) for idx, circle in enumerate(circles_array): azimuth_center, polar_center, radius_angle = circle - theta = np.linspace(0, 2 * np.pi, num_points) - azimuth_points = azimuth_center + radius_angle * np.cos(theta) - polar_points = polar_center + radius_angle * np.sin(theta) - azimuth_points = np.mod(azimuth_points, 2 * np.pi) - polar_points = np.clip(polar_points, -np.pi / 2, np.pi / 2) + azimuth_points, polar_points = generate_spherical_small_circle_points( + azimuth_center, polar_center, radius_angle, num_points + ) + x_points, y_points = mercator_transform(azimuth_points, polar_points) radius_deg = np.degrees(radius_angle) label = f"Circle {idx + 1}: {radius_deg:.1f}°" - ax1.scatter(azimuth_points, polar_points, s=10, alpha=0.6, label=label) + ax1.scatter(x_points, y_points, s=10, alpha=0.6, label=label) - ax1.set_xlabel("Azimuthal Angle (radians)", fontsize=10) - ax1.set_ylabel("Polar Angle (radians)", fontsize=10) + ax1.set_xlabel("Azmuth (radians)", fontsize=10) + ax1.set_ylabel("Polar (radians)", fontsize=10) ax1.set_xlim(0, 2 * np.pi) - ax1.set_ylim(-np.pi / 2, np.pi / 2) ax1.grid(True, alpha=0.3) - ax1.set_title("Flattened View", fontsize=12) - ax1.legend(fontsize=8) + ax1.set_title("Mercator Projection", fontsize=12) + # No legend in Mercator view - # 2. Polar view - ax2 = fig.add_subplot(132, projection='polar') + # 2. Gauss-Krüger view + ax2 = fig.add_subplot(132) + central_meridian = np.pi # Central meridian at 180° for idx, circle in enumerate(circles_array): azimuth_center, polar_center, radius_angle = circle - theta = np.linspace(0, 2 * np.pi, num_points) - azimuth_points = azimuth_center + radius_angle * np.cos(theta) - polar_points = polar_center + radius_angle * np.sin(theta) - azimuth_points = np.mod(azimuth_points, 2 * np.pi) - polar_points = np.clip(polar_points, -np.pi / 2, np.pi / 2) - radius_points = (polar_points + np.pi / 2) / np.pi + azimuth_points, polar_points = generate_spherical_small_circle_points( + azimuth_center, polar_center, radius_angle, num_points + ) + easting, northing = gauss_kruger_transform( + azimuth_points, polar_points, central_meridian + ) radius_deg = np.degrees(radius_angle) label = f"Circle {idx + 1}: {radius_deg:.1f}°" - ax2.scatter(azimuth_points, radius_points, s=10, alpha=0.6, label=label) + # Swap coordinates: x=northing (azimuth-related), y=easting (polar-related) + ax2.scatter(northing, easting, s=10, alpha=0.6, label=label) - ax2.set_title("Polar View", fontsize=12, pad=20) + ax2.set_xlabel("Azimuth-related (radians)", fontsize=10) + ax2.set_ylabel("Polar-related (radians)", fontsize=10) + ax2.set_title("Gauss-Krüger Projection", fontsize=12) + ax2.set_aspect("equal", adjustable="box") ax2.grid(True, alpha=0.3) - ax2.legend(fontsize=8, loc='upper right', bbox_to_anchor=(1.2, 1.1)) + # No legend in Gauss-Krüger view # 3. 3D view - ax3 = fig.add_subplot(133, projection='3d') + ax3 = fig.add_subplot(133, projection="3d") # Plot sphere surface u = np.linspace(0, 2 * np.pi, 30) @@ -535,15 +506,13 @@ def plot_spherical_small_circles_multiplot( x_sphere = np.outer(np.cos(v), np.cos(u)) y_sphere = np.outer(np.cos(v), np.sin(u)) z_sphere = np.outer(np.sin(v), np.ones(np.size(u))) - ax3.plot_surface(x_sphere, y_sphere, z_sphere, alpha=0.1, color='lightblue') + ax3.plot_surface(x_sphere, y_sphere, z_sphere, alpha=0.1, color="lightblue") for idx, circle in enumerate(circles_array): azimuth_center, polar_center, radius_angle = circle - theta = np.linspace(0, 2 * np.pi, num_points) - azimuth_points = azimuth_center + radius_angle * np.cos(theta) - polar_points = polar_center + radius_angle * np.sin(theta) - azimuth_points = np.mod(azimuth_points, 2 * np.pi) - polar_points = np.clip(polar_points, -np.pi / 2, np.pi / 2) + azimuth_points, polar_points = generate_spherical_small_circle_points( + azimuth_center, polar_center, radius_angle, num_points + ) x = np.cos(polar_points) * np.cos(azimuth_points) y = np.cos(polar_points) * np.sin(azimuth_points) z = np.sin(polar_points) @@ -551,9 +520,9 @@ def plot_spherical_small_circles_multiplot( label = f"Circle {idx + 1}: {radius_deg:.1f}°" ax3.scatter(x, y, z, s=10, alpha=0.8, label=label) - ax3.set_xlabel('X', fontsize=10) - ax3.set_ylabel('Y', fontsize=10) - ax3.set_zlabel('Z', fontsize=10) + ax3.set_xlabel("X", fontsize=10) + ax3.set_ylabel("Y", fontsize=10) + ax3.set_zlabel("Z", fontsize=10) ax3.set_title("3D View", fontsize=12) ax3.legend(fontsize=8) ax3.set_box_aspect([1, 1, 1]) @@ -563,38 +532,73 @@ def plot_spherical_small_circles_multiplot( return fig -def demo(): +def demo() -> None: """Demonstration of multiplot view showing spherical small circles in three different projections.""" # Define 4 circles using UnitSphericalSmallCircle dataclass - circles = [ + circles: List[UnitSphericalSmallCircle] = [ UnitSphericalSmallCircle( azimuth=0, polar=0, - radius_angle=0.6, # Circle at origin (0°, equator) with large radius + radius_angle=np.deg2rad(45), ), UnitSphericalSmallCircle( - azimuth=np.pi, # 180° - polar=np.pi / 4, # 45° north - radius_angle=0.2, + azimuth=0, + polar=np.pi / 2, + radius_angle=np.deg2rad(45), ), UnitSphericalSmallCircle( - azimuth=3 * np.pi / 2, # 270° - polar=-np.pi / 6, # -30° south - radius_angle=0.25, + azimuth=0, + polar=np.pi, + radius_angle=np.deg2rad(45), + ), + UnitSphericalSmallCircle( + azimuth=0, + polar=np.deg2rad(-90), + radius_angle=np.deg2rad(45), + ), + UnitSphericalSmallCircle( + azimuth=np.deg2rad(90), + polar=0, + radius_angle=np.deg2rad(45), + ), + UnitSphericalSmallCircle( + azimuth=np.deg2rad(-90), + polar=0, + radius_angle=np.deg2rad(45), ), UnitSphericalSmallCircle( - azimuth=0, # 0° - polar=-np.pi / 2, # -90° south pole - radius_angle=np.pi / 2, + azimuth=np.deg2rad(30), + polar=np.deg2rad(30), + radius_angle=np.deg2rad(30), ), + # UnitSphericalSmallCircle( + # azimuth=np.deg2rad(10), + # polar=np.deg2rad(-80), + # radius_angle=np.deg2rad(45), + # ), + # UnitSphericalSmallCircle( + # azimuth=np.deg2rad(10), + # polar=np.deg2rad(80), + # radius_angle=np.deg2rad(45), + # ), + # UnitSphericalSmallCircle( + # azimuth=np.pi, # 180° + # polar=np.pi / 4, # 45° north + # radius_angle=0.2, + # ), + # UnitSphericalSmallCircle( + # azimuth=3 * np.pi / 2, # 270° + # polar=-np.pi / 6, # -30° south + # radius_angle=0.25, + # ), ] # Create the multiplot showing all three views fig = plot_spherical_small_circles_multiplot( circles, num_points=240, - title="Unit Spherical Small Circles - Flattened, Polar, and 3D Views" + title="Unit Spherical Small Circles - Mercator, Gauss-Krüger, and 3D Views", ) # Show the plot From 8a94967c7d2e6e9e34dc277d1670f36ee076116f Mon Sep 17 00:00:00 2001 From: kopecn Date: Mon, 1 Dec 2025 16:35:22 -0800 Subject: [PATCH 02/22] still have some issues. --- examples/plotQuatUnitCircles.py | 27 ++ pyproject.toml | 1 + src/pyMathTools/constructors/__init__.py | 0 .../constructUnitSphericalSmallCircle.py | 24 ++ .../generators/sphericalGenerators.py | 4 +- .../transforms/sphericalTransforms.py | 2 +- .../plotUnitSphericalSmallCircles.py | 239 +++++++++++------- 7 files changed, 198 insertions(+), 99 deletions(-) create mode 100644 examples/plotQuatUnitCircles.py create mode 100644 src/pyMathTools/constructors/__init__.py create mode 100644 src/pyMathTools/constructors/constructUnitSphericalSmallCircle.py diff --git a/examples/plotQuatUnitCircles.py b/examples/plotQuatUnitCircles.py new file mode 100644 index 0000000..1882bb1 --- /dev/null +++ b/examples/plotQuatUnitCircles.py @@ -0,0 +1,27 @@ +import numpy as np +import quaternion + +from pyMathTools.constructors.constructUnitSphericalSmallCircle import ( + quat_to_unitSphericalSmallCircle, +) + +from pyMathToolsPlotHelpers.plotUnitSphericalSmallCircles import ( + plot_spherical_small_circles_multiplot, +) + + +def main(): + + _ = plot_spherical_small_circles_multiplot( + circles=[ + quat_to_unitSphericalSmallCircle(q=quaternion.one, radius_angle=np.pi / 8), + quat_to_unitSphericalSmallCircle(q=quaternion.x, radius_angle=np.pi / 8), + quat_to_unitSphericalSmallCircle(q=quaternion.y), + quat_to_unitSphericalSmallCircle(q=quaternion.z), + ], + show_plot=True, + ) + + +if __name__ == "__main__": + main() diff --git a/pyproject.toml b/pyproject.toml index 09e2439..afb03c6 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -22,6 +22,7 @@ license = {text = "MIT"} dependencies = [ "pyFoundationTools @ git+https://github.com/kopecn/py-foundationTools.git@dev", "matplotlib", + "numpy-quaternion", ] [project.optional-dependencies] diff --git a/src/pyMathTools/constructors/__init__.py b/src/pyMathTools/constructors/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/src/pyMathTools/constructors/constructUnitSphericalSmallCircle.py b/src/pyMathTools/constructors/constructUnitSphericalSmallCircle.py new file mode 100644 index 0000000..6c6d167 --- /dev/null +++ b/src/pyMathTools/constructors/constructUnitSphericalSmallCircle.py @@ -0,0 +1,24 @@ +import numpy as np +from numpy import atan2, acos +import quaternion +from numpy import pi +from foundationTypes.mathTypes.UnitSphericalSmallCircle import UnitSphericalSmallCircle + + +def quat_to_unitSphericalSmallCircle( + q: np.quaternion, radius_angle: float = pi / 4 +) -> UnitSphericalSmallCircle: + # Extract vector part + + r = np.sqrt(q.x**2 + q.y**2 + q.z**2) + + if r == 0: + # Default to north pole if vector part is zero + azimuth = 0.0 + polar = pi / 2 + else: + polar = atan2(q.y, q.x) + azimuth = acos(q.z / r) + return UnitSphericalSmallCircle( + azimuth=polar, polar=azimuth, radius_angle=radius_angle + ) diff --git a/src/pyMathTools/generators/sphericalGenerators.py b/src/pyMathTools/generators/sphericalGenerators.py index 3ecbe52..5e00eb8 100644 --- a/src/pyMathTools/generators/sphericalGenerators.py +++ b/src/pyMathTools/generators/sphericalGenerators.py @@ -34,9 +34,9 @@ def generate_spherical_small_circle_points( Parameters: ----------- azimuth_center : float or NDArray - Azimuth angle of circle center (radians, 0 to 2π) + Azimuth angle of circle center (radians, -2π to 2π) polar_center : float or NDArray - Polar angle of circle center (radians, -π/2 to π/2) + Polar angle of circle center (radians, 0 to π) radius_angle : float or NDArray Angular radius of the circle (radians) num_points : int diff --git a/src/pyMathTools/transforms/sphericalTransforms.py b/src/pyMathTools/transforms/sphericalTransforms.py index 79b7b1b..d046bcd 100644 --- a/src/pyMathTools/transforms/sphericalTransforms.py +++ b/src/pyMathTools/transforms/sphericalTransforms.py @@ -41,7 +41,7 @@ def mercator_transform( def gauss_kruger_transform( azimuth: FloatOrNDArray, polar: FloatOrNDArray, - central_meridian: float = pi, + central_meridian: float = 0, ) -> tuple[FloatOrNDArray, FloatOrNDArray]: """ Apply Gauss-Krüger (Transverse Mercator) projection to spherical coordinates. diff --git a/src/pyMathToolsPlotHelpers/plotUnitSphericalSmallCircles.py b/src/pyMathToolsPlotHelpers/plotUnitSphericalSmallCircles.py index 47be1ce..4da19dd 100644 --- a/src/pyMathToolsPlotHelpers/plotUnitSphericalSmallCircles.py +++ b/src/pyMathToolsPlotHelpers/plotUnitSphericalSmallCircles.py @@ -21,6 +21,7 @@ def plot_spherical_small_circles( figsize: Tuple[int, int] = (10, 6), title: str = "Small Circles on Sphere (Mercator Projection)", show_legend: bool = True, + show_plot: bool = False, ) -> Tuple[Figure, Axes]: """ Plot small circles on a sphere using Mercator projection. @@ -93,6 +94,9 @@ def plot_spherical_small_circles( plt.tight_layout() + if show_plot: + plt.show() + return fig, ax @@ -104,6 +108,8 @@ def plot_spherical_small_circles_advanced( colors: Optional[List[str]] = None, labels: Optional[List[str]] = None, show_centers: bool = True, + show_plot: bool = False, + ax: Optional[Axes] = None, ) -> Tuple[Figure, Axes]: """ Advanced version with more customization options. @@ -115,7 +121,7 @@ def plot_spherical_small_circles_advanced( num_points : int, optional Number of points to generate around each circle figsize : tuple, optional - Figure size + Figure size (ignored if ax is provided) title : str, optional Plot title colors : list, optional @@ -124,12 +130,22 @@ def plot_spherical_small_circles_advanced( List of labels for each circle show_centers : bool, optional Whether to show center points (default: True) + show_plot : bool, optional + Whether to show the plot (default: False, ignored if ax is provided) + ax : matplotlib.axes.Axes, optional + Axes to plot on. If None, creates new figure and axes Returns: -------- fig, ax : matplotlib figure and axes objects """ - fig, ax = plt.subplots(figsize=figsize) + # Track whether we created the axes or it was provided + created_ax = ax is None + + if ax is None: + fig, ax = plt.subplots(figsize=figsize) + else: + fig = ax.get_figure() # Convert UnitSphericalSmallCircle objects to arrays if needed if circles and isinstance(circles[0], UnitSphericalSmallCircle): @@ -186,7 +202,11 @@ def plot_spherical_small_circles_advanced( if labels is not None: ax.legend() - plt.tight_layout() + # Only call tight_layout and show if we created the axes + if created_ax: + plt.tight_layout() + if show_plot: + plt.show() return fig, ax @@ -196,8 +216,10 @@ def plot_spherical_small_circles_gauss_kruger( num_points: int = 120, figsize: Tuple[int, int] = (10, 6), title: str = "Small Circles on Sphere (Gauss-Krüger Projection)", - central_meridian: float = np.pi, + central_meridian: float = 0, show_legend: bool = True, + show_plot: bool = False, + ax: Optional[Axes] = None, ) -> Tuple[Figure, Axes]: """ Plot small circles on a sphere using Gauss-Krüger (Transverse Mercator) projection. @@ -209,19 +231,29 @@ def plot_spherical_small_circles_gauss_kruger( num_points : int, optional Number of points to generate around each circle (default: 120) figsize : tuple, optional - Figure size (default: (10, 6)) + Figure size (default: (10, 6), ignored if ax is provided) title : str, optional Plot title central_meridian : float, optional - Central meridian in radians (default: π) + Central meridian in radians (default: 0) show_legend : bool, optional Whether to show the legend (default: True) + show_plot : bool, optional + Whether to show the plot (default: False, ignored if ax is provided) + ax : matplotlib.axes.Axes, optional + Axes to plot on. If None, creates new figure and axes Returns: -------- fig, ax : matplotlib figure and axes objects """ - fig, ax = plt.subplots(figsize=figsize) + # Track whether we created the axes or it was provided + created_ax = ax is None + + if ax is None: + fig, ax = plt.subplots(figsize=figsize) + else: + fig = ax.get_figure() # Convert UnitSphericalSmallCircle objects to arrays if needed if circles and isinstance(circles[0], UnitSphericalSmallCircle): @@ -251,7 +283,7 @@ def plot_spherical_small_circles_gauss_kruger( # Plot with swapped coordinates to view down at equator # x-axis (azimuth-related) = northing # y-axis (polar-related) = easting - ax.scatter(northing, easting, s=10, alpha=0.6, label=label) + ax.scatter(easting, northing, s=10, alpha=0.6, label=label) ax.set_xlabel("Azimuth-related (radians)", fontsize=12) ax.set_ylabel("Polar-related (radians)", fontsize=12) @@ -263,7 +295,11 @@ def plot_spherical_small_circles_gauss_kruger( if show_legend: ax.legend() - plt.tight_layout() + # Only call tight_layout and show if we created the axes + if created_ax: + plt.tight_layout() + if show_plot: + plt.show() return fig, ax @@ -273,6 +309,7 @@ def plot_spherical_small_circles_polar( num_points: int = 120, figsize: Tuple[int, int] = (8, 8), title: str = "Small Circles on Sphere (Polar View)", + show_plot: bool = False, ) -> Tuple[Figure, Axes]: """ Plot small circles on a sphere using polar projection. @@ -328,6 +365,9 @@ def plot_spherical_small_circles_polar( plt.tight_layout() + if show_plot: + plt.show() + return fig, ax @@ -338,6 +378,8 @@ def plot_spherical_small_circles_3d( title: str = "Small Circles on Unit Sphere (3D)", show_sphere: bool = True, show_legend: bool = True, + show_plot: bool = False, + ax: Optional[Axes] = None, ) -> Tuple[Figure, Axes]: """ Plot small circles on a 3D unit sphere. @@ -349,20 +391,30 @@ def plot_spherical_small_circles_3d( num_points : int, optional Number of points to generate around each circle (default: 120) figsize : tuple, optional - Figure size (default: (10, 10)) + Figure size (default: (10, 10), ignored if ax is provided) title : str, optional Plot title show_sphere : bool, optional Whether to show the sphere surface (default: True) show_legend : bool, optional Whether to show the legend (default: True) + show_plot : bool, optional + Whether to show the plot (default: False, ignored if ax is provided) + ax : matplotlib.axes.Axes, optional + 3D axes to plot on. If None, creates new figure and axes Returns: -------- fig, ax : matplotlib figure and axes objects """ - fig = plt.figure(figsize=figsize) - ax = fig.add_subplot(111, projection="3d") + # Track whether we created the axes or it was provided + created_ax = ax is None + + if ax is None: + fig = plt.figure(figsize=figsize) + ax = fig.add_subplot(111, projection="3d") + else: + fig = ax.get_figure() # Convert UnitSphericalSmallCircle objects to arrays if needed if circles and isinstance(circles[0], UnitSphericalSmallCircle): @@ -397,7 +449,11 @@ def plot_spherical_small_circles_3d( # Convert radius to degrees for legend radius_deg = np.degrees(radius_angle) - label = f"Circle {idx + 1}: {radius_deg:.1f}°" + azimuth_deg = np.degrees(azimuth_center) + polar_deg = np.degrees(polar_center) + label = ( + f"Circle {idx + 1}: {azimuth_deg:.1f}°, {polar_deg:.1f}°, {radius_deg:.1f}°" + ) # Plot the circle ax.scatter(x, y, z, s=10, alpha=0.8, label=label) @@ -407,14 +463,22 @@ def plot_spherical_small_circles_3d( ax.set_zlabel("Z") ax.set_title(title, fontsize=14) - # Add legend if requested + # Add legend if requested, positioned outside the plot area if show_legend: - ax.legend() + ax.legend(loc="upper left", bbox_to_anchor=(1.05, 1.0)) # Set equal aspect ratio ax.set_box_aspect([1, 1, 1]) - plt.tight_layout() + # Set initial view angle: physics convention with Z up, Y right, X down-left + # elev=20 looks from slightly above, azim=45 gives the proper orientation + ax.view_init(elev=20, azim=45) + + # Only call tight_layout and show if we created the axes + if created_ax: + plt.tight_layout() + if show_plot: + plt.show() return fig, ax @@ -424,6 +488,7 @@ def plot_spherical_small_circles_multiplot( num_points: int = 120, figsize: Tuple[int, int] = (18, 6), title: str = "Small Circles on Unit Sphere - Multiple Views", + show_plot: bool = False, ) -> Figure: """ Create a multiplot showing Mercator, Gauss-Krüger, and 3D views of spherical small circles. @@ -438,97 +503,71 @@ def plot_spherical_small_circles_multiplot( Figure size (default: (18, 6)) title : str, optional Overall title + show_plot : bool, optional + Whether to show the plot (default: False) Returns: -------- fig : matplotlib figure object """ + # Create figure with overall title fig = plt.figure(figsize=figsize) fig.suptitle(title, fontsize=16) - # Convert UnitSphericalSmallCircle objects to arrays if needed - circles_array = circles - if circles and isinstance(circles[0], UnitSphericalSmallCircle): - circles_array = [[c.azimuth, c.polar, c.radius_angle] for c in circles] - - circles_array = np.array(circles_array) - if circles_array.ndim == 1: - circles_array = circles_array.reshape(1, -1) - - # 1. Mercator projection view - ax1 = fig.add_subplot(131) - for idx, circle in enumerate(circles_array): - azimuth_center, polar_center, radius_angle = circle - azimuth_points, polar_points = generate_spherical_small_circle_points( - azimuth_center, polar_center, radius_angle, num_points - ) - x_points, y_points = mercator_transform(azimuth_points, polar_points) - radius_deg = np.degrees(radius_angle) - label = f"Circle {idx + 1}: {radius_deg:.1f}°" - ax1.scatter(x_points, y_points, s=10, alpha=0.6, label=label) - + # Create subplots + ax1 = fig.add_subplot(131) # Mercator projection + ax2 = fig.add_subplot(132) # Gauss-Krüger projection + ax3 = fig.add_subplot(133, projection="3d") # 3D view + + # 1. Mercator projection using advanced plot function + plot_spherical_small_circles_advanced( + circles=circles, + num_points=num_points, + title="Mercator Projection", + show_centers=False, + ax=ax1, + ) + # Customize for multiplot ax1.set_xlabel("Azmuth (radians)", fontsize=10) ax1.set_ylabel("Polar (radians)", fontsize=10) - ax1.set_xlim(0, 2 * np.pi) - ax1.grid(True, alpha=0.3) ax1.set_title("Mercator Projection", fontsize=12) - # No legend in Mercator view - - # 2. Gauss-Krüger view - ax2 = fig.add_subplot(132) - central_meridian = np.pi # Central meridian at 180° - for idx, circle in enumerate(circles_array): - azimuth_center, polar_center, radius_angle = circle - azimuth_points, polar_points = generate_spherical_small_circle_points( - azimuth_center, polar_center, radius_angle, num_points - ) - easting, northing = gauss_kruger_transform( - azimuth_points, polar_points, central_meridian - ) - radius_deg = np.degrees(radius_angle) - label = f"Circle {idx + 1}: {radius_deg:.1f}°" - # Swap coordinates: x=northing (azimuth-related), y=easting (polar-related) - ax2.scatter(northing, easting, s=10, alpha=0.6, label=label) + # 2. Gauss-Krüger projection + plot_spherical_small_circles_gauss_kruger( + circles=circles, + num_points=num_points, + title="Gauss-Krüger Projection", + central_meridian=0, + show_legend=False, + ax=ax2, + ) + # Customize for multiplot ax2.set_xlabel("Azimuth-related (radians)", fontsize=10) ax2.set_ylabel("Polar-related (radians)", fontsize=10) ax2.set_title("Gauss-Krüger Projection", fontsize=12) - ax2.set_aspect("equal", adjustable="box") - ax2.grid(True, alpha=0.3) - # No legend in Gauss-Krüger view # 3. 3D view - ax3 = fig.add_subplot(133, projection="3d") - - # Plot sphere surface - u = np.linspace(0, 2 * np.pi, 30) - v = np.linspace(-np.pi / 2, np.pi / 2, 20) - x_sphere = np.outer(np.cos(v), np.cos(u)) - y_sphere = np.outer(np.cos(v), np.sin(u)) - z_sphere = np.outer(np.sin(v), np.ones(np.size(u))) - ax3.plot_surface(x_sphere, y_sphere, z_sphere, alpha=0.1, color="lightblue") - - for idx, circle in enumerate(circles_array): - azimuth_center, polar_center, radius_angle = circle - azimuth_points, polar_points = generate_spherical_small_circle_points( - azimuth_center, polar_center, radius_angle, num_points - ) - x = np.cos(polar_points) * np.cos(azimuth_points) - y = np.cos(polar_points) * np.sin(azimuth_points) - z = np.sin(polar_points) - radius_deg = np.degrees(radius_angle) - label = f"Circle {idx + 1}: {radius_deg:.1f}°" - ax3.scatter(x, y, z, s=10, alpha=0.8, label=label) - + plot_spherical_small_circles_3d( + circles=circles, + num_points=num_points, + title="3D View", + show_sphere=True, + show_legend=True, + ax=ax3, + ) + # Customize for multiplot ax3.set_xlabel("X", fontsize=10) ax3.set_ylabel("Y", fontsize=10) ax3.set_zlabel("Z", fontsize=10) ax3.set_title("3D View", fontsize=12) - ax3.legend(fontsize=8) - ax3.set_box_aspect([1, 1, 1]) + # Position legend outside plot area with smaller font + ax3.legend(loc="upper left", bbox_to_anchor=(1.05, 1.0), fontsize=8) plt.tight_layout() + if show_plot: + plt.show() + return fig @@ -543,28 +582,28 @@ def demo() -> None: radius_angle=np.deg2rad(45), ), UnitSphericalSmallCircle( - azimuth=0, - polar=np.pi / 2, + azimuth=np.pi / 2, + polar=0, radius_angle=np.deg2rad(45), ), UnitSphericalSmallCircle( - azimuth=0, - polar=np.pi, + azimuth=np.pi, + polar=0, radius_angle=np.deg2rad(45), ), UnitSphericalSmallCircle( - azimuth=0, - polar=np.deg2rad(-90), + azimuth=np.deg2rad(-90), + polar=0, radius_angle=np.deg2rad(45), ), UnitSphericalSmallCircle( - azimuth=np.deg2rad(90), - polar=0, + azimuth=0, + polar=np.deg2rad(90), radius_angle=np.deg2rad(45), ), UnitSphericalSmallCircle( - azimuth=np.deg2rad(-90), - polar=0, + azimuth=0, + polar=np.deg2rad(-90), radius_angle=np.deg2rad(45), ), UnitSphericalSmallCircle( @@ -595,14 +634,22 @@ def demo() -> None: ] # Create the multiplot showing all three views - fig = plot_spherical_small_circles_multiplot( + _ = plot_spherical_small_circles_multiplot( circles, num_points=240, title="Unit Spherical Small Circles - Mercator, Gauss-Krüger, and 3D Views", + show_plot=True, + ) + + _ = plot_spherical_small_circles( + circles=circles, + show_plot=True, ) - # Show the plot - plt.show() + _ = plot_spherical_small_circles_polar( + circles=circles, + show_plot=True, + ) if __name__ == "__main__": From c5575595eefcf69fb6878b7d7ec71f6c1c669374 Mon Sep 17 00:00:00 2001 From: kopecn Date: Mon, 1 Dec 2025 18:16:08 -0800 Subject: [PATCH 03/22] fixed... --- .../generators/sphericalGenerators.py | 18 +++++++++++------- 1 file changed, 11 insertions(+), 7 deletions(-) diff --git a/src/pyMathTools/generators/sphericalGenerators.py b/src/pyMathTools/generators/sphericalGenerators.py index 5e00eb8..0b2b869 100644 --- a/src/pyMathTools/generators/sphericalGenerators.py +++ b/src/pyMathTools/generators/sphericalGenerators.py @@ -10,7 +10,7 @@ linspace, zeros, sqrt, - arcsin, + arccos, dot, ) from numpy.linalg import norm @@ -49,10 +49,12 @@ def generate_spherical_small_circle_points( polar_points : NDArray Polar angles of points on the circle """ - # Convert center to Cartesian coordinates - x_center = cos(polar_center) * cos(azimuth_center) - y_center = cos(polar_center) * sin(azimuth_center) - z_center = sin(polar_center) + # Convert center to Cartesian coordinates (physics convention) + # theta (polar) = 0 at north pole, pi at south pole + # phi (azimuth) = rotation about Z axis + x_center = sin(polar_center) * cos(azimuth_center) + y_center = sin(polar_center) * sin(azimuth_center) + z_center = cos(polar_center) center = array([x_center, y_center, z_center]) # Create orthonormal basis at the center point @@ -102,12 +104,14 @@ def generate_spherical_small_circle_points( # If no rotation needed (shouldn't happen) point = center - # Convert back to spherical coordinates + # Convert back to spherical coordinates (physics convention) x, y, z = point r = sqrt(x**2 + y**2 + z**2) x, y, z = x / r, y / r, z / r # Normalize to unit sphere - polar_points[i] = arcsin(clip(z, -1, 1)) + # theta (polar) = arccos(z), ranges from 0 (north pole) to pi (south pole) + polar_points[i] = arccos(clip(z, -1, 1)) + # phi (azimuth) = atan2(y, x) azimuth_points[i] = arctan2(y, x) if azimuth_points[i] < 0: azimuth_points[i] += 2 * pi From 3c3a5a3d7f6161ef751b10b8d7c6ba1314e52019 Mon Sep 17 00:00:00 2001 From: kopecn Date: Mon, 1 Dec 2025 18:16:35 -0800 Subject: [PATCH 04/22] keep the 3d fixes --- .../plotUnitSphericalSmallCircles.py | 47 ++++++++++--------- 1 file changed, 26 insertions(+), 21 deletions(-) diff --git a/src/pyMathToolsPlotHelpers/plotUnitSphericalSmallCircles.py b/src/pyMathToolsPlotHelpers/plotUnitSphericalSmallCircles.py index 4da19dd..eae500e 100644 --- a/src/pyMathToolsPlotHelpers/plotUnitSphericalSmallCircles.py +++ b/src/pyMathToolsPlotHelpers/plotUnitSphericalSmallCircles.py @@ -424,13 +424,13 @@ def plot_spherical_small_circles_3d( if circles.ndim == 1: circles = circles.reshape(1, -1) - # Optionally plot the sphere surface + # Optionally plot the sphere surface (physics convention) if show_sphere: - u = np.linspace(0, 2 * np.pi, 30) - v = np.linspace(-np.pi / 2, np.pi / 2, 20) - x_sphere = np.outer(np.cos(v), np.cos(u)) - y_sphere = np.outer(np.cos(v), np.sin(u)) - z_sphere = np.outer(np.sin(v), np.ones(np.size(u))) + phi = np.linspace(0, 2 * np.pi, 30) # azimuth + theta = np.linspace(0, np.pi, 20) # polar angle from north to south + x_sphere = np.outer(np.sin(theta), np.cos(phi)) + y_sphere = np.outer(np.sin(theta), np.sin(phi)) + z_sphere = np.outer(np.cos(theta), np.ones(np.size(phi))) ax.plot_surface(x_sphere, y_sphere, z_sphere, alpha=0.1, color="lightblue") # Plot each circle @@ -442,18 +442,18 @@ def plot_spherical_small_circles_3d( azimuth_center, polar_center, radius_angle, num_points ) - # Convert to Cartesian coordinates (unit sphere, r=1) - x = np.cos(polar_points) * np.cos(azimuth_points) - y = np.cos(polar_points) * np.sin(azimuth_points) - z = np.sin(polar_points) + # Convert to Cartesian coordinates (unit sphere, r=1, physics convention) + # theta (polar) from 0 (north pole) to pi (south pole) + # phi (azimuth) rotation about Z axis + x = np.sin(polar_points) * np.cos(azimuth_points) + y = np.sin(polar_points) * np.sin(azimuth_points) + z = np.cos(polar_points) # Convert radius to degrees for legend radius_deg = np.degrees(radius_angle) azimuth_deg = np.degrees(azimuth_center) polar_deg = np.degrees(polar_center) - label = ( - f"Circle {idx + 1}: {azimuth_deg:.1f}°, {polar_deg:.1f}°, {radius_deg:.1f}°" - ) + label = f"Circle {idx + 1}: A:{azimuth_deg:.1f}°, P:{polar_deg:.1f}°, R:{radius_deg:.1f}°" # Plot the circle ax.scatter(x, y, z, s=10, alpha=0.8, label=label) @@ -581,29 +581,34 @@ def demo() -> None: polar=0, radius_angle=np.deg2rad(45), ), - UnitSphericalSmallCircle( + UnitSphericalSmallCircle( # should overlap circle 1 azimuth=np.pi / 2, polar=0, radius_angle=np.deg2rad(45), ), UnitSphericalSmallCircle( - azimuth=np.pi, - polar=0, + azimuth=0, + polar=np.pi / 2, radius_angle=np.deg2rad(45), ), UnitSphericalSmallCircle( - azimuth=np.deg2rad(-90), - polar=0, + azimuth=0, + polar=np.pi, radius_angle=np.deg2rad(45), ), UnitSphericalSmallCircle( azimuth=0, - polar=np.deg2rad(90), + polar=np.pi, # Fixed: was -π/2 (invalid), now π (south pole) radius_angle=np.deg2rad(45), ), UnitSphericalSmallCircle( - azimuth=0, - polar=np.deg2rad(-90), + azimuth=-np.pi / 2, + polar=np.pi / 2, + radius_angle=np.deg2rad(45), + ), + UnitSphericalSmallCircle( + azimuth=np.pi / 2, + polar=np.pi / 2, radius_angle=np.deg2rad(45), ), UnitSphericalSmallCircle( From c0be50b09e5a221975a92f917e48ec125d6908d6 Mon Sep 17 00:00:00 2001 From: kopecn Date: Mon, 1 Dec 2025 18:20:29 -0800 Subject: [PATCH 05/22] this is correct... --- .../plotUnitSphericalSmallCircles.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/pyMathToolsPlotHelpers/plotUnitSphericalSmallCircles.py b/src/pyMathToolsPlotHelpers/plotUnitSphericalSmallCircles.py index eae500e..9ec2a3e 100644 --- a/src/pyMathToolsPlotHelpers/plotUnitSphericalSmallCircles.py +++ b/src/pyMathToolsPlotHelpers/plotUnitSphericalSmallCircles.py @@ -581,11 +581,6 @@ def demo() -> None: polar=0, radius_angle=np.deg2rad(45), ), - UnitSphericalSmallCircle( # should overlap circle 1 - azimuth=np.pi / 2, - polar=0, - radius_angle=np.deg2rad(45), - ), UnitSphericalSmallCircle( azimuth=0, polar=np.pi / 2, @@ -598,7 +593,7 @@ def demo() -> None: ), UnitSphericalSmallCircle( azimuth=0, - polar=np.pi, # Fixed: was -π/2 (invalid), now π (south pole) + polar=-np.pi / 2, radius_angle=np.deg2rad(45), ), UnitSphericalSmallCircle( @@ -616,6 +611,11 @@ def demo() -> None: polar=np.deg2rad(30), radius_angle=np.deg2rad(30), ), + UnitSphericalSmallCircle( # should overlap circle 1 + azimuth=np.pi / 2, + polar=0, + radius_angle=np.deg2rad(45), + ), # UnitSphericalSmallCircle( # azimuth=np.deg2rad(10), # polar=np.deg2rad(-80), From 907a0c04de7bd0902a8caded261d91ebf8f24ed3 Mon Sep 17 00:00:00 2001 From: kopecn Date: Mon, 1 Dec 2025 18:58:44 -0800 Subject: [PATCH 06/22] fixes for plate-carree projection --- .../transforms/sphericalTransforms.py | 45 +++++++----- .../plotUnitSphericalSmallCircles.py | 71 ++++++++++--------- 2 files changed, 65 insertions(+), 51 deletions(-) diff --git a/src/pyMathTools/transforms/sphericalTransforms.py b/src/pyMathTools/transforms/sphericalTransforms.py index d046bcd..8b642fe 100644 --- a/src/pyMathTools/transforms/sphericalTransforms.py +++ b/src/pyMathTools/transforms/sphericalTransforms.py @@ -10,31 +10,36 @@ from pyMathTools.hints import FloatOrNDArray -def mercator_transform( +def plate_carree_transform( azimuth: FloatOrNDArray, polar: FloatOrNDArray ) -> tuple[FloatOrNDArray, FloatOrNDArray]: """ - Apply Mercator projection to spherical coordinates. + Apply plate carrée (equirectangular) projection to spherical coordinates (physics convention). + + This is a simple cylindrical projection that preserves angular measurements, + making it suitable for robot positioning where actual angles in radians are needed. Parameters: ----------- azimuth : float or ndarray - Azimuth angle(s) in radians (0 to 2π) + Azimuth angle(s) in radians (0 to 2π), rotation about Z axis polar : float or ndarray - Polar angle(s) in radians (-π/2 to π/2) + Polar angle(s) in radians (0 to π), angle from north pole Returns: -------- - x : float or ndarray - Mercator x coordinate (longitude) + Azimuth in radians (0 to 2π) - y : float or ndarray - Polar (radians) coordinate (latitude projection) + Latitude in radians (-π/2 to π/2), where: + - π/2 = north pole + - 0 = equator + - -π/2 = south pole """ x = azimuth - # Mercator projection: y = ln(tan(π/4 + lat/2)) - # Clamp polar to avoid infinities at poles - polar_clamped = clip(polar, -pi / 2 + 0.001, pi / 2 - 0.001) - y = log(tan(pi / 4 + polar_clamped / 2)) + # Convert physics polar angle (theta, 0 to π) to latitude (-π/2 to π/2) + # lat = π/2 - theta + y = pi / 2 - polar return x, y @@ -44,16 +49,16 @@ def gauss_kruger_transform( central_meridian: float = 0, ) -> tuple[FloatOrNDArray, FloatOrNDArray]: """ - Apply Gauss-Krüger (Transverse Mercator) projection to spherical coordinates. + Apply Gauss-Krüger (Transverse Mercator) projection to spherical coordinates (physics convention). Parameters: ----------- azimuth : float or ndarray - Azimuth angle(s) in radians (0 to 2π) + Azimuth angle(s) in radians (0 to 2π), rotation about Z axis polar : float or ndarray - Polar angle(s) in radians (-π/2 to π/2) + Polar angle(s) in radians (0 to π), angle from north pole central_meridian : float, optional - Central meridian in radians (default: π) + Central meridian in radians (default: 0) Returns: -------- @@ -68,14 +73,18 @@ def gauss_kruger_transform( # Normalize to [-π, π] lambda_rel = arctan2(sin(lambda_rel), cos(lambda_rel)) - # Clamp polar to avoid numerical issues - polar_clamped = clip(polar, -pi / 2 + 0.001, pi / 2 - 0.001) + # Convert physics polar angle (theta, 0 to π) to latitude (-π/2 to π/2) + # lat = π/2 - theta + latitude = pi / 2 - polar + + # Clamp latitude to avoid numerical issues + lat_clamped = clip(latitude, -pi / 2 + 0.001, pi / 2 - 0.001) # Transverse Mercator formulas for a sphere # Easting - easting = arctanh(cos(polar_clamped) * sin(lambda_rel)) + easting = arctanh(cos(lat_clamped) * sin(lambda_rel)) # Northing - northing = arctan2(sin(polar_clamped), cos(polar_clamped) * cos(lambda_rel)) + northing = arctan2(sin(lat_clamped), cos(lat_clamped) * cos(lambda_rel)) return easting, northing diff --git a/src/pyMathToolsPlotHelpers/plotUnitSphericalSmallCircles.py b/src/pyMathToolsPlotHelpers/plotUnitSphericalSmallCircles.py index 9ec2a3e..0a38a69 100644 --- a/src/pyMathToolsPlotHelpers/plotUnitSphericalSmallCircles.py +++ b/src/pyMathToolsPlotHelpers/plotUnitSphericalSmallCircles.py @@ -7,7 +7,7 @@ from foundationTypes.mathTypes.UnitSphericalSmallCircle import UnitSphericalSmallCircle from pyMathTools.transforms.sphericalTransforms import ( - mercator_transform, + plate_carree_transform, gauss_kruger_transform, ) from pyMathTools.generators.sphericalGenerators import ( @@ -19,12 +19,12 @@ def plot_spherical_small_circles( circles: Union[List[UnitSphericalSmallCircle], List[List[float]]], num_points: int = 120, figsize: Tuple[int, int] = (10, 6), - title: str = "Small Circles on Sphere (Mercator Projection)", + title: str = "Small Circles on Sphere (Plate Carrée Projection)", show_legend: bool = True, show_plot: bool = False, ) -> Tuple[Figure, Axes]: """ - Plot small circles on a sphere using Mercator projection. + Plot small circles on a sphere using plate carrée (equirectangular) projection. Parameters: ----------- @@ -64,19 +64,21 @@ def plot_spherical_small_circles( azimuth_center, polar_center, radius_angle, num_points ) - # Apply Mercator projection - x_points, y_points = mercator_transform(azimuth_points, polar_points) + # Apply Plate Carrée projection + x_points, y_points = plate_carree_transform(azimuth_points, polar_points) # Convert radius to degrees for legend radius_deg = np.degrees(radius_angle) - label = f"Circle {idx + 1}: {radius_deg:.1f}°" + azimuth_deg = np.degrees(azimuth_center) + polar_deg = np.degrees(polar_center) + label = f"Circle {idx + 1}: A:{azimuth_deg:.1f}°, P:{polar_deg:.1f}°, R:{radius_deg:.1f}°" # Plot the circle ax.scatter(x_points, y_points, s=10, alpha=0.6, label=label) # Set up axes - ax.set_xlabel("Azmuth (radians)", fontsize=12) - ax.set_ylabel("Polar (radians)", fontsize=12) + ax.set_xlabel("Azimuth (radians)", fontsize=12) + ax.set_ylabel("Latitude (radians)", fontsize=12) ax.set_xlim(0, 2 * np.pi) # Add grid @@ -163,8 +165,8 @@ def plot_spherical_small_circles_advanced( azimuth_center, polar_center, radius_angle, num_points ) - # Apply Mercator projection - x_points, y_points = mercator_transform(azimuth_points, polar_points) + # Apply Plate Carrée projection + x_points, y_points = plate_carree_transform(azimuth_points, polar_points) # Determine color and label color = colors[idx] if colors is not None and idx < len(colors) else None @@ -175,7 +177,7 @@ def plot_spherical_small_circles_advanced( # Plot center point if requested if show_centers: - x_center, y_center = mercator_transform(azimuth_center, polar_center) + x_center, y_center = plate_carree_transform(azimuth_center, polar_center) ax.scatter( x_center, y_center, @@ -186,8 +188,8 @@ def plot_spherical_small_circles_advanced( ) # Set up axes - ax.set_xlabel("Azmuth (radians)", fontsize=12) - ax.set_ylabel("Polar (radians)", fontsize=12) + ax.set_xlabel("Azimuth (radians)", fontsize=12) + ax.set_ylabel("Latitude (radians)", fontsize=12) ax.set_xlim(0, 2 * np.pi) # Add grid @@ -278,15 +280,16 @@ def plot_spherical_small_circles_gauss_kruger( # Convert radius to degrees for legend radius_deg = np.degrees(radius_angle) - label = f"Circle {idx + 1}: {radius_deg:.1f}°" + azimuth_deg = np.degrees(azimuth_center) + polar_deg = np.degrees(polar_center) + label = f"Circle {idx + 1}: A:{azimuth_deg:.1f}°, P:{polar_deg:.1f}°, R:{radius_deg:.1f}°" - # Plot with swapped coordinates to view down at equator - # x-axis (azimuth-related) = northing - # y-axis (polar-related) = easting + # Plot in Gauss-Krüger projection coordinates + # x-axis = easting, y-axis = northing ax.scatter(easting, northing, s=10, alpha=0.6, label=label) - ax.set_xlabel("Azimuth-related (radians)", fontsize=12) - ax.set_ylabel("Polar-related (radians)", fontsize=12) + ax.set_xlabel("Easting", fontsize=12) + ax.set_ylabel("Northing", fontsize=12) ax.grid(True, alpha=0.3) ax.set_title(title, fontsize=14) ax.set_aspect("equal", adjustable="box") @@ -354,7 +357,9 @@ def plot_spherical_small_circles_polar( # Convert radius to degrees for legend radius_deg = np.degrees(radius_angle) - label = f"Circle {idx + 1}: {radius_deg:.1f}°" + azimuth_deg = np.degrees(azimuth_center) + polar_deg = np.degrees(polar_center) + label = f"Circle {idx + 1}: A:{azimuth_deg:.1f}°, P:{polar_deg:.1f}°, R:{radius_deg:.1f}°" # Plot in polar coordinates ax.scatter(azimuth_points, radius_points, s=10, alpha=0.6, label=label) @@ -491,7 +496,7 @@ def plot_spherical_small_circles_multiplot( show_plot: bool = False, ) -> Figure: """ - Create a multiplot showing Mercator, Gauss-Krüger, and 3D views of spherical small circles. + Create a multiplot showing Plate Carrée, Gauss-Krüger, and 3D views of spherical small circles. Parameters: ----------- @@ -515,22 +520,22 @@ def plot_spherical_small_circles_multiplot( fig.suptitle(title, fontsize=16) # Create subplots - ax1 = fig.add_subplot(131) # Mercator projection + ax1 = fig.add_subplot(131) # Plate Carrée projection ax2 = fig.add_subplot(132) # Gauss-Krüger projection ax3 = fig.add_subplot(133, projection="3d") # 3D view - # 1. Mercator projection using advanced plot function + # 1. Plate Carrée projection using advanced plot function plot_spherical_small_circles_advanced( circles=circles, num_points=num_points, - title="Mercator Projection", + title="Plate Carrée Projection", show_centers=False, ax=ax1, ) # Customize for multiplot - ax1.set_xlabel("Azmuth (radians)", fontsize=10) - ax1.set_ylabel("Polar (radians)", fontsize=10) - ax1.set_title("Mercator Projection", fontsize=12) + ax1.set_xlabel("Azimuth (radians)", fontsize=10) + ax1.set_ylabel("Latitude (radians)", fontsize=10) + ax1.set_title("Plate Carrée Projection", fontsize=12) # 2. Gauss-Krüger projection plot_spherical_small_circles_gauss_kruger( @@ -611,11 +616,11 @@ def demo() -> None: polar=np.deg2rad(30), radius_angle=np.deg2rad(30), ), - UnitSphericalSmallCircle( # should overlap circle 1 - azimuth=np.pi / 2, - polar=0, - radius_angle=np.deg2rad(45), - ), + # UnitSphericalSmallCircle( # should overlap circle 1 + # azimuth=np.pi / 2, + # polar=0, + # radius_angle=np.deg2rad(45), + # ), # UnitSphericalSmallCircle( # azimuth=np.deg2rad(10), # polar=np.deg2rad(-80), @@ -642,7 +647,7 @@ def demo() -> None: _ = plot_spherical_small_circles_multiplot( circles, num_points=240, - title="Unit Spherical Small Circles - Mercator, Gauss-Krüger, and 3D Views", + title="Unit Spherical Small Circles - Plate Carrée, Gauss-Krüger, and 3D Views", show_plot=True, ) From a5f26ac187e420bb566adc2f0ab8362ad2ea81d1 Mon Sep 17 00:00:00 2001 From: kopecn Date: Mon, 1 Dec 2025 19:43:47 -0800 Subject: [PATCH 07/22] humph... --- .../transforms/sphericalTransforms.py | 49 +---- .../plotUnitSphericalSmallCircles.py | 186 +++++------------- 2 files changed, 47 insertions(+), 188 deletions(-) diff --git a/src/pyMathTools/transforms/sphericalTransforms.py b/src/pyMathTools/transforms/sphericalTransforms.py index 8b642fe..a4fe0fa 100644 --- a/src/pyMathTools/transforms/sphericalTransforms.py +++ b/src/pyMathTools/transforms/sphericalTransforms.py @@ -6,7 +6,7 @@ geometry, mapping, and scientific computing. """ -from numpy import log, clip, pi, tan, arctan2, cos, sin, arctanh +from numpy import pi from pyMathTools.hints import FloatOrNDArray @@ -41,50 +41,3 @@ def plate_carree_transform( # lat = π/2 - theta y = pi / 2 - polar return x, y - - -def gauss_kruger_transform( - azimuth: FloatOrNDArray, - polar: FloatOrNDArray, - central_meridian: float = 0, -) -> tuple[FloatOrNDArray, FloatOrNDArray]: - """ - Apply Gauss-Krüger (Transverse Mercator) projection to spherical coordinates (physics convention). - - Parameters: - ----------- - azimuth : float or ndarray - Azimuth angle(s) in radians (0 to 2π), rotation about Z axis - polar : float or ndarray - Polar angle(s) in radians (0 to π), angle from north pole - central_meridian : float, optional - Central meridian in radians (default: 0) - - Returns: - -------- - easting : float or ndarray - Gauss-Krüger easting coordinate - northing : float or ndarray - Gauss-Krüger northing coordinate - """ - # Adjust azimuth relative to central meridian - lambda_rel = azimuth - central_meridian - - # Normalize to [-π, π] - lambda_rel = arctan2(sin(lambda_rel), cos(lambda_rel)) - - # Convert physics polar angle (theta, 0 to π) to latitude (-π/2 to π/2) - # lat = π/2 - theta - latitude = pi / 2 - polar - - # Clamp latitude to avoid numerical issues - lat_clamped = clip(latitude, -pi / 2 + 0.001, pi / 2 - 0.001) - - # Transverse Mercator formulas for a sphere - # Easting - easting = arctanh(cos(lat_clamped) * sin(lambda_rel)) - - # Northing - northing = arctan2(sin(lat_clamped), cos(lat_clamped) * cos(lambda_rel)) - - return easting, northing diff --git a/src/pyMathToolsPlotHelpers/plotUnitSphericalSmallCircles.py b/src/pyMathToolsPlotHelpers/plotUnitSphericalSmallCircles.py index 0a38a69..4199584 100644 --- a/src/pyMathToolsPlotHelpers/plotUnitSphericalSmallCircles.py +++ b/src/pyMathToolsPlotHelpers/plotUnitSphericalSmallCircles.py @@ -6,10 +6,7 @@ from foundationTypes.mathTypes.UnitSphericalSmallCircle import UnitSphericalSmallCircle -from pyMathTools.transforms.sphericalTransforms import ( - plate_carree_transform, - gauss_kruger_transform, -) +from pyMathTools.transforms.sphericalTransforms import plate_carree_transform from pyMathTools.generators.sphericalGenerators import ( generate_spherical_small_circle_points, ) @@ -74,19 +71,25 @@ def plot_spherical_small_circles( label = f"Circle {idx + 1}: A:{azimuth_deg:.1f}°, P:{polar_deg:.1f}°, R:{radius_deg:.1f}°" # Plot the circle - ax.scatter(x_points, y_points, s=10, alpha=0.6, label=label) + ax.scatter( + np.degrees(x_points), + np.degrees(y_points), + s=10, + alpha=0.6, + label=label, + ) # Set up axes - ax.set_xlabel("Azimuth (radians)", fontsize=12) - ax.set_ylabel("Latitude (radians)", fontsize=12) - ax.set_xlim(0, 2 * np.pi) + ax.set_xlabel("Azimuth (deg)", fontsize=12) + ax.set_ylabel("Latitude (deg)", fontsize=12) + ax.set_xlim(0, 180) # Add grid ax.grid(True, alpha=0.3) - # Set x-axis ticks to show multiples of pi - ax.set_xticks([0, np.pi / 2, np.pi, 3 * np.pi / 2, 2 * np.pi]) - ax.set_xticklabels([r"$0$", r"$\pi/2$", r"$\pi$", r"$3\pi/2$", r"$2\pi$"]) + # Set x-axis ticks + ax.set_xticks([0, 90, 180, 270, 360]) + ax.set_xticklabels(["0°", "90°", "180°", "270°", "360°"]) ax.set_title(title, fontsize=14) @@ -173,14 +176,21 @@ def plot_spherical_small_circles_advanced( label = labels[idx] if labels is not None and idx < len(labels) else None # Plot the circle - ax.scatter(x_points, y_points, s=10, alpha=0.6, c=color, label=label) + ax.scatter( + np.degrees(x_points), + np.degrees(y_points), + s=10, + alpha=0.6, + c=color, + label=label, + ) # Plot center point if requested if show_centers: x_center, y_center = plate_carree_transform(azimuth_center, polar_center) ax.scatter( - x_center, - y_center, + np.degrees(x_center), + np.degrees(y_center), s=100, marker="x", c=color if color else "red", @@ -188,16 +198,16 @@ def plot_spherical_small_circles_advanced( ) # Set up axes - ax.set_xlabel("Azimuth (radians)", fontsize=12) - ax.set_ylabel("Latitude (radians)", fontsize=12) - ax.set_xlim(0, 2 * np.pi) + ax.set_xlabel("Azimuth (deg)", fontsize=12) + ax.set_ylabel("Latitude (deg)", fontsize=12) + ax.set_xlim(0, 180) # Add grid ax.grid(True, alpha=0.3) # Set ticks - ax.set_xticks([0, np.pi / 2, np.pi, 3 * np.pi / 2, 2 * np.pi]) - ax.set_xticklabels([r"$0$", r"$\pi/2$", r"$\pi$", r"$3\pi/2$", r"$2\pi$"]) + ax.set_xticks([0, 90, 180, 270, 360]) + ax.set_xticklabels(["0°", "90°", "180°", "270°", "360°"]) ax.set_title(title, fontsize=14) @@ -213,100 +223,6 @@ def plot_spherical_small_circles_advanced( return fig, ax -def plot_spherical_small_circles_gauss_kruger( - circles: Union[List[UnitSphericalSmallCircle], List[List[float]]], - num_points: int = 120, - figsize: Tuple[int, int] = (10, 6), - title: str = "Small Circles on Sphere (Gauss-Krüger Projection)", - central_meridian: float = 0, - show_legend: bool = True, - show_plot: bool = False, - ax: Optional[Axes] = None, -) -> Tuple[Figure, Axes]: - """ - Plot small circles on a sphere using Gauss-Krüger (Transverse Mercator) projection. - - Parameters: - ----------- - circles : List[UnitSphericalSmallCircle] or array-like - List of UnitSphericalSmallCircle objects or array of [azimuth, polar, radius] lists - num_points : int, optional - Number of points to generate around each circle (default: 120) - figsize : tuple, optional - Figure size (default: (10, 6), ignored if ax is provided) - title : str, optional - Plot title - central_meridian : float, optional - Central meridian in radians (default: 0) - show_legend : bool, optional - Whether to show the legend (default: True) - show_plot : bool, optional - Whether to show the plot (default: False, ignored if ax is provided) - ax : matplotlib.axes.Axes, optional - Axes to plot on. If None, creates new figure and axes - - Returns: - -------- - fig, ax : matplotlib figure and axes objects - """ - # Track whether we created the axes or it was provided - created_ax = ax is None - - if ax is None: - fig, ax = plt.subplots(figsize=figsize) - else: - fig = ax.get_figure() - - # Convert UnitSphericalSmallCircle objects to arrays if needed - if circles and isinstance(circles[0], UnitSphericalSmallCircle): - circles = [[c.azimuth, c.polar, c.radius_angle] for c in circles] - - circles = np.array(circles) - if circles.ndim == 1: - circles = circles.reshape(1, -1) - - for idx, circle in enumerate(circles): - azimuth_center, polar_center, radius_angle = circle - - # Generate true spherical small circle points - azimuth_points, polar_points = generate_spherical_small_circle_points( - azimuth_center, polar_center, radius_angle, num_points - ) - - # Apply Gauss-Krüger projection - easting, northing = gauss_kruger_transform( - azimuth_points, polar_points, central_meridian - ) - - # Convert radius to degrees for legend - radius_deg = np.degrees(radius_angle) - azimuth_deg = np.degrees(azimuth_center) - polar_deg = np.degrees(polar_center) - label = f"Circle {idx + 1}: A:{azimuth_deg:.1f}°, P:{polar_deg:.1f}°, R:{radius_deg:.1f}°" - - # Plot in Gauss-Krüger projection coordinates - # x-axis = easting, y-axis = northing - ax.scatter(easting, northing, s=10, alpha=0.6, label=label) - - ax.set_xlabel("Easting", fontsize=12) - ax.set_ylabel("Northing", fontsize=12) - ax.grid(True, alpha=0.3) - ax.set_title(title, fontsize=14) - ax.set_aspect("equal", adjustable="box") - - # Add legend if requested - if show_legend: - ax.legend() - - # Only call tight_layout and show if we created the axes - if created_ax: - plt.tight_layout() - if show_plot: - plt.show() - - return fig, ax - - def plot_spherical_small_circles_polar( circles: Union[List[UnitSphericalSmallCircle], List[List[float]]], num_points: int = 120, @@ -520,9 +436,8 @@ def plot_spherical_small_circles_multiplot( fig.suptitle(title, fontsize=16) # Create subplots - ax1 = fig.add_subplot(131) # Plate Carrée projection - ax2 = fig.add_subplot(132) # Gauss-Krüger projection - ax3 = fig.add_subplot(133, projection="3d") # 3D view + ax1 = fig.add_subplot(121) # Plate Carrée projection + ax2 = fig.add_subplot(122, projection="3d") # 3D view # 1. Plate Carrée projection using advanced plot function plot_spherical_small_circles_advanced( @@ -537,36 +452,22 @@ def plot_spherical_small_circles_multiplot( ax1.set_ylabel("Latitude (radians)", fontsize=10) ax1.set_title("Plate Carrée Projection", fontsize=12) - # 2. Gauss-Krüger projection - plot_spherical_small_circles_gauss_kruger( - circles=circles, - num_points=num_points, - title="Gauss-Krüger Projection", - central_meridian=0, - show_legend=False, - ax=ax2, - ) - # Customize for multiplot - ax2.set_xlabel("Azimuth-related (radians)", fontsize=10) - ax2.set_ylabel("Polar-related (radians)", fontsize=10) - ax2.set_title("Gauss-Krüger Projection", fontsize=12) - - # 3. 3D view + # 2. 3D view plot_spherical_small_circles_3d( circles=circles, num_points=num_points, title="3D View", show_sphere=True, show_legend=True, - ax=ax3, + ax=ax2, ) # Customize for multiplot - ax3.set_xlabel("X", fontsize=10) - ax3.set_ylabel("Y", fontsize=10) - ax3.set_zlabel("Z", fontsize=10) - ax3.set_title("3D View", fontsize=12) + ax2.set_xlabel("X", fontsize=10) + ax2.set_ylabel("Y", fontsize=10) + ax2.set_zlabel("Z", fontsize=10) + ax2.set_title("3D View", fontsize=12) # Position legend outside plot area with smaller font - ax3.legend(loc="upper left", bbox_to_anchor=(1.05, 1.0), fontsize=8) + ax2.legend(loc="upper left", bbox_to_anchor=(1.05, 1.0), fontsize=8) plt.tight_layout() @@ -612,9 +513,14 @@ def demo() -> None: radius_angle=np.deg2rad(45), ), UnitSphericalSmallCircle( - azimuth=np.deg2rad(30), - polar=np.deg2rad(30), - radius_angle=np.deg2rad(30), + azimuth=-np.deg2rad(45), + polar=np.deg2rad(45), + radius_angle=np.deg2rad(45), + ), + UnitSphericalSmallCircle( + azimuth=np.deg2rad(135), + polar=np.deg2rad(45), + radius_angle=np.deg2rad(45), ), # UnitSphericalSmallCircle( # should overlap circle 1 # azimuth=np.pi / 2, From 29ce6153e9e0f89b71e2285b60c0e9fd10a4ce1c Mon Sep 17 00:00:00 2001 From: kopecn Date: Mon, 1 Dec 2025 19:52:37 -0800 Subject: [PATCH 08/22] reducing code to use the UnitSphericalSmallCircle --- .../generators/sphericalGenerators.py | 20 +++-- .../plotUnitSphericalSmallCircles.py | 73 +++++-------------- 2 files changed, 27 insertions(+), 66 deletions(-) diff --git a/src/pyMathTools/generators/sphericalGenerators.py b/src/pyMathTools/generators/sphericalGenerators.py index 0b2b869..01abcb3 100644 --- a/src/pyMathTools/generators/sphericalGenerators.py +++ b/src/pyMathTools/generators/sphericalGenerators.py @@ -14,15 +14,13 @@ dot, ) from numpy.linalg import norm +from foundationTypes.mathTypes.UnitSphericalSmallCircle import UnitSphericalSmallCircle - -from pyMathTools.hints import FloatOrNDArray, FloatNDArray +from pyMathTools.hints import FloatNDArray def generate_spherical_small_circle_points( - azimuth_center: FloatOrNDArray, - polar_center: FloatOrNDArray, - radius_angle: FloatOrNDArray, + circle: UnitSphericalSmallCircle, num_points: int = 120, ) -> Tuple[FloatNDArray, FloatNDArray]: """ @@ -33,7 +31,7 @@ def generate_spherical_small_circle_points( Parameters: ----------- - azimuth_center : float or NDArray + circle.azimuth : float or NDArray Azimuth angle of circle center (radians, -2π to 2π) polar_center : float or NDArray Polar angle of circle center (radians, 0 to π) @@ -52,9 +50,9 @@ def generate_spherical_small_circle_points( # Convert center to Cartesian coordinates (physics convention) # theta (polar) = 0 at north pole, pi at south pole # phi (azimuth) = rotation about Z axis - x_center = sin(polar_center) * cos(azimuth_center) - y_center = sin(polar_center) * sin(azimuth_center) - z_center = cos(polar_center) + x_center = sin(circle.polar) * cos(circle.azimuth) + y_center = sin(circle.polar) * sin(circle.azimuth) + z_center = cos(circle.polar) center = array([x_center, y_center, z_center]) # Create orthonormal basis at the center point @@ -92,8 +90,8 @@ def generate_spherical_small_circle_points( if axis_norm > 1e-10: axis = axis / axis_norm # Rotate center around axis by radius_angle - cos_angle = cos(radius_angle) - sin_angle = sin(radius_angle) + cos_angle = cos(circle.radius_angle) + sin_angle = sin(circle.radius_angle) point = ( cos_angle * center diff --git a/src/pyMathToolsPlotHelpers/plotUnitSphericalSmallCircles.py b/src/pyMathToolsPlotHelpers/plotUnitSphericalSmallCircles.py index 4199584..a461d2f 100644 --- a/src/pyMathToolsPlotHelpers/plotUnitSphericalSmallCircles.py +++ b/src/pyMathToolsPlotHelpers/plotUnitSphericalSmallCircles.py @@ -13,7 +13,7 @@ def plot_spherical_small_circles( - circles: Union[List[UnitSphericalSmallCircle], List[List[float]]], + circles: List[UnitSphericalSmallCircle], num_points: int = 120, figsize: Tuple[int, int] = (10, 6), title: str = "Small Circles on Sphere (Plate Carrée Projection)", @@ -43,31 +43,21 @@ def plot_spherical_small_circles( """ fig, ax = plt.subplots(figsize=figsize) - # Convert UnitSphericalSmallCircle objects to arrays if needed - if circles and isinstance(circles[0], UnitSphericalSmallCircle): - circles = [[c.azimuth, c.polar, c.radius_angle] for c in circles] - - # Convert to numpy array - circles = np.array(circles) - if circles.ndim == 1: - circles = circles.reshape(1, -1) - # Generate points for each circle for idx, circle in enumerate(circles): - azimuth_center, polar_center, radius_angle = circle # Generate true spherical small circle points azimuth_points, polar_points = generate_spherical_small_circle_points( - azimuth_center, polar_center, radius_angle, num_points + circle, num_points ) # Apply Plate Carrée projection x_points, y_points = plate_carree_transform(azimuth_points, polar_points) # Convert radius to degrees for legend - radius_deg = np.degrees(radius_angle) - azimuth_deg = np.degrees(azimuth_center) - polar_deg = np.degrees(polar_center) + radius_deg = np.degrees(circle.radius_angle) + azimuth_deg = np.degrees(circle.azimuth) + polar_deg = np.degrees(circle.polar) label = f"Circle {idx + 1}: A:{azimuth_deg:.1f}°, P:{polar_deg:.1f}°, R:{radius_deg:.1f}°" # Plot the circle @@ -106,7 +96,7 @@ def plot_spherical_small_circles( def plot_spherical_small_circles_advanced( - circles: Union[List[UnitSphericalSmallCircle], List[List[float]]], + circles: List[UnitSphericalSmallCircle], num_points: int = 120, figsize: Tuple[int, int] = (10, 6), title: str = "Small Circles on Sphere (Flattened)", @@ -152,20 +142,11 @@ def plot_spherical_small_circles_advanced( else: fig = ax.get_figure() - # Convert UnitSphericalSmallCircle objects to arrays if needed - if circles and isinstance(circles[0], UnitSphericalSmallCircle): - circles = [[c.azimuth, c.polar, c.radius_angle] for c in circles] - - circles = np.array(circles) - if circles.ndim == 1: - circles = circles.reshape(1, -1) - for idx, circle in enumerate(circles): - azimuth_center, polar_center, radius_angle = circle # Generate true spherical small circle points azimuth_points, polar_points = generate_spherical_small_circle_points( - azimuth_center, polar_center, radius_angle, num_points + circle, num_points ) # Apply Plate Carrée projection @@ -187,7 +168,7 @@ def plot_spherical_small_circles_advanced( # Plot center point if requested if show_centers: - x_center, y_center = plate_carree_transform(azimuth_center, polar_center) + x_center, y_center = plate_carree_transform(circle.azimuth, circle.polar) ax.scatter( np.degrees(x_center), np.degrees(y_center), @@ -224,7 +205,7 @@ def plot_spherical_small_circles_advanced( def plot_spherical_small_circles_polar( - circles: Union[List[UnitSphericalSmallCircle], List[List[float]]], + circles: List[UnitSphericalSmallCircle], num_points: int = 120, figsize: Tuple[int, int] = (8, 8), title: str = "Small Circles on Sphere (Polar View)", @@ -251,20 +232,11 @@ def plot_spherical_small_circles_polar( """ fig, ax = plt.subplots(figsize=figsize, subplot_kw={"projection": "polar"}) - # Convert UnitSphericalSmallCircle objects to arrays if needed - if circles and isinstance(circles[0], UnitSphericalSmallCircle): - circles = [[c.azimuth, c.polar, c.radius_angle] for c in circles] - - circles = np.array(circles) - if circles.ndim == 1: - circles = circles.reshape(1, -1) - for idx, circle in enumerate(circles): - azimuth_center, polar_center, radius_angle = circle # Generate true spherical small circle points azimuth_points, polar_points = generate_spherical_small_circle_points( - azimuth_center, polar_center, radius_angle, num_points + circle, num_points ) # Map polar angle to radius (0 to 1) for polar plot @@ -272,9 +244,9 @@ def plot_spherical_small_circles_polar( radius_points = (polar_points + np.pi / 2) / np.pi # Convert radius to degrees for legend - radius_deg = np.degrees(radius_angle) - azimuth_deg = np.degrees(azimuth_center) - polar_deg = np.degrees(polar_center) + radius_deg = np.degrees(circle.radius_angle) + azimuth_deg = np.degrees(circle.azimuth) + polar_deg = np.degrees(circle.polar) label = f"Circle {idx + 1}: A:{azimuth_deg:.1f}°, P:{polar_deg:.1f}°, R:{radius_deg:.1f}°" # Plot in polar coordinates @@ -293,7 +265,7 @@ def plot_spherical_small_circles_polar( def plot_spherical_small_circles_3d( - circles: Union[List[UnitSphericalSmallCircle], List[List[float]]], + circles: List[UnitSphericalSmallCircle], num_points: int = 120, figsize: Tuple[int, int] = (10, 10), title: str = "Small Circles on Unit Sphere (3D)", @@ -337,14 +309,6 @@ def plot_spherical_small_circles_3d( else: fig = ax.get_figure() - # Convert UnitSphericalSmallCircle objects to arrays if needed - if circles and isinstance(circles[0], UnitSphericalSmallCircle): - circles = [[c.azimuth, c.polar, c.radius_angle] for c in circles] - - circles = np.array(circles) - if circles.ndim == 1: - circles = circles.reshape(1, -1) - # Optionally plot the sphere surface (physics convention) if show_sphere: phi = np.linspace(0, 2 * np.pi, 30) # azimuth @@ -356,11 +320,10 @@ def plot_spherical_small_circles_3d( # Plot each circle for idx, circle in enumerate(circles): - azimuth_center, polar_center, radius_angle = circle # Generate true spherical small circle points azimuth_points, polar_points = generate_spherical_small_circle_points( - azimuth_center, polar_center, radius_angle, num_points + circle, num_points ) # Convert to Cartesian coordinates (unit sphere, r=1, physics convention) @@ -371,9 +334,9 @@ def plot_spherical_small_circles_3d( z = np.cos(polar_points) # Convert radius to degrees for legend - radius_deg = np.degrees(radius_angle) - azimuth_deg = np.degrees(azimuth_center) - polar_deg = np.degrees(polar_center) + radius_deg = np.degrees(circle.radius_angle) + azimuth_deg = np.degrees(circle.azimuth) + polar_deg = np.degrees(circle.polar) label = f"Circle {idx + 1}: A:{azimuth_deg:.1f}°, P:{polar_deg:.1f}°, R:{radius_deg:.1f}°" # Plot the circle From fd109dd1e265557b79f2174d193ed983329e373f Mon Sep 17 00:00:00 2001 From: kopecn Date: Mon, 1 Dec 2025 19:59:49 -0800 Subject: [PATCH 09/22] humph... --- examples/plotQuatUnitCircles.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/examples/plotQuatUnitCircles.py b/examples/plotQuatUnitCircles.py index 1882bb1..0addbef 100644 --- a/examples/plotQuatUnitCircles.py +++ b/examples/plotQuatUnitCircles.py @@ -1,4 +1,4 @@ -import numpy as np +from numpy import pi import quaternion from pyMathTools.constructors.constructUnitSphericalSmallCircle import ( @@ -14,8 +14,8 @@ def main(): _ = plot_spherical_small_circles_multiplot( circles=[ - quat_to_unitSphericalSmallCircle(q=quaternion.one, radius_angle=np.pi / 8), - quat_to_unitSphericalSmallCircle(q=quaternion.x, radius_angle=np.pi / 8), + quat_to_unitSphericalSmallCircle(q=quaternion.one, radius_angle=pi / 8), + quat_to_unitSphericalSmallCircle(q=quaternion.x, radius_angle=pi / 8), quat_to_unitSphericalSmallCircle(q=quaternion.y), quat_to_unitSphericalSmallCircle(q=quaternion.z), ], From e85d3c224972d3746496930b7245bc3c6620addf Mon Sep 17 00:00:00 2001 From: kopecn Date: Mon, 1 Dec 2025 20:35:00 -0800 Subject: [PATCH 10/22] cont... --- examples/plotQuatUnitCircles.py | 36 ++++++++ .../generators/sphericalGenerators.py | 6 +- .../plotUnitSphericalSmallCircles.py | 82 ++++++++++++------- 3 files changed, 88 insertions(+), 36 deletions(-) diff --git a/examples/plotQuatUnitCircles.py b/examples/plotQuatUnitCircles.py index 0addbef..e5905c1 100644 --- a/examples/plotQuatUnitCircles.py +++ b/examples/plotQuatUnitCircles.py @@ -1,3 +1,30 @@ +""" +Example: Visualizing Quaternion-Derived Small Circles on the Unit Sphere + +This script demonstrates how to convert quaternions to spherical small circles +and visualize them using multiple projections. It serves as a validation and +walkthrough for the quaternion-to-spherical conversion and plotting utilities. + +Steps: +------ +1. Import required modules and functions. +2. Convert several canonical quaternions (identity and axis-aligned) to + UnitSphericalSmallCircle objects. +3. Plot these circles using the multiplot function, which shows Mercator, + Gauss-Krüger, and 3D views. + +Usage: +------ +Run this script directly to see the resulting plots and validate the conversion. + +Expected Output: +---------------- +- The identity quaternion (`quaternion.one`) should map to the north pole. +- `quaternion.x`, `quaternion.y`, and `quaternion.z` should map to circles + centered on the X, Y, and Z axes, respectively. + +""" + from numpy import pi import quaternion @@ -11,7 +38,16 @@ def main(): + """ + Convert canonical quaternions to spherical small circles and plot them. + + - quaternion.one: Identity quaternion (should map to north pole) + - quaternion.x: X-axis quaternion + - quaternion.y: Y-axis quaternion + - quaternion.z: Z-axis quaternion + All circles use a radius angle of pi/8 except y and z, which use the default. + """ _ = plot_spherical_small_circles_multiplot( circles=[ quat_to_unitSphericalSmallCircle(q=quaternion.one, radius_angle=pi / 8), diff --git a/src/pyMathTools/generators/sphericalGenerators.py b/src/pyMathTools/generators/sphericalGenerators.py index 01abcb3..b50e67a 100644 --- a/src/pyMathTools/generators/sphericalGenerators.py +++ b/src/pyMathTools/generators/sphericalGenerators.py @@ -71,7 +71,7 @@ def generate_spherical_small_circle_points( north_tangent = north_tangent / norm(north_tangent) # Generate points around the circle - theta = linspace(0, 2 * pi, num_points) + theta = linspace(-pi, pi, num_points) # Points on the small circle using rotation # The small circle is at angular distance radius_angle from center @@ -107,11 +107,7 @@ def generate_spherical_small_circle_points( r = sqrt(x**2 + y**2 + z**2) x, y, z = x / r, y / r, z / r # Normalize to unit sphere - # theta (polar) = arccos(z), ranges from 0 (north pole) to pi (south pole) polar_points[i] = arccos(clip(z, -1, 1)) - # phi (azimuth) = atan2(y, x) azimuth_points[i] = arctan2(y, x) - if azimuth_points[i] < 0: - azimuth_points[i] += 2 * pi return azimuth_points, polar_points diff --git a/src/pyMathToolsPlotHelpers/plotUnitSphericalSmallCircles.py b/src/pyMathToolsPlotHelpers/plotUnitSphericalSmallCircles.py index a461d2f..436bd36 100644 --- a/src/pyMathToolsPlotHelpers/plotUnitSphericalSmallCircles.py +++ b/src/pyMathToolsPlotHelpers/plotUnitSphericalSmallCircles.py @@ -72,14 +72,14 @@ def plot_spherical_small_circles( # Set up axes ax.set_xlabel("Azimuth (deg)", fontsize=12) ax.set_ylabel("Latitude (deg)", fontsize=12) - ax.set_xlim(0, 180) + ax.set_xlim(-180, 180) # Add grid ax.grid(True, alpha=0.3) # Set x-axis ticks - ax.set_xticks([0, 90, 180, 270, 360]) - ax.set_xticklabels(["0°", "90°", "180°", "270°", "360°"]) + ax.set_xticks([-180, -90, 0, 90, 180]) + ax.set_xticklabels(["-180°", "-90°", "0°", "90°", "180°"]) ax.set_title(title, fontsize=14) @@ -181,14 +181,14 @@ def plot_spherical_small_circles_advanced( # Set up axes ax.set_xlabel("Azimuth (deg)", fontsize=12) ax.set_ylabel("Latitude (deg)", fontsize=12) - ax.set_xlim(0, 180) + ax.set_xlim(-180, 180) # Add grid ax.grid(True, alpha=0.3) # Set ticks - ax.set_xticks([0, 90, 180, 270, 360]) - ax.set_xticklabels(["0°", "90°", "180°", "270°", "360°"]) + ax.set_xticks([-180, -90, 0, 90, 180]) + ax.set_xticklabels(["-180°", "-90°", "0°", "90°", "180°"]) ax.set_title(title, fontsize=14) @@ -209,6 +209,8 @@ def plot_spherical_small_circles_polar( num_points: int = 120, figsize: Tuple[int, int] = (8, 8), title: str = "Small Circles on Sphere (Polar View)", + labels: Optional[List[str]] = None, + ax: Optional[Axes] = None, show_plot: bool = False, ) -> Tuple[Figure, Axes]: """ @@ -230,7 +232,13 @@ def plot_spherical_small_circles_polar( -------- fig, ax : matplotlib figure and axes objects """ - fig, ax = plt.subplots(figsize=figsize, subplot_kw={"projection": "polar"}) + # Track whether we created the axes or it was provided + created_ax = ax is None + + if ax is None: + fig, ax = plt.subplots(figsize=figsize, subplot_kw={"projection": "polar"}) + else: + fig = ax.get_figure() for idx, circle in enumerate(circles): @@ -253,13 +261,17 @@ def plot_spherical_small_circles_polar( ax.scatter(azimuth_points, radius_points, s=10, alpha=0.6, label=label) ax.set_title(title, fontsize=14, pad=20) - ax.legend(loc="upper right", bbox_to_anchor=(1.3, 1.1)) - ax.grid(True, alpha=0.3) - plt.tight_layout() + if labels is not None: + ax.legend(loc="upper right", bbox_to_anchor=(1.3, 1.1)) - if show_plot: - plt.show() + ax.grid(True, alpha=0.3) + + # Only call tight_layout and show if we created the axes + if created_ax: + plt.tight_layout() + if show_plot: + plt.show() return fig, ax @@ -370,12 +382,12 @@ def plot_spherical_small_circles_3d( def plot_spherical_small_circles_multiplot( circles: Union[List[UnitSphericalSmallCircle], List[List[float]]], num_points: int = 120, - figsize: Tuple[int, int] = (18, 6), + figsize: Tuple[int, int] = (24, 6), title: str = "Small Circles on Unit Sphere - Multiple Views", show_plot: bool = False, ) -> Figure: """ - Create a multiplot showing Plate Carrée, Gauss-Krüger, and 3D views of spherical small circles. + Create a multiplot showing Plate Carrée, Top View, and 3D views of spherical small circles. Parameters: ----------- @@ -384,7 +396,7 @@ def plot_spherical_small_circles_multiplot( num_points : int, optional Number of points to generate around each circle (default: 120) figsize : tuple, optional - Figure size (default: (18, 6)) + Figure size (default: (24, 6)) title : str, optional Overall title show_plot : bool, optional @@ -399,8 +411,9 @@ def plot_spherical_small_circles_multiplot( fig.suptitle(title, fontsize=16) # Create subplots - ax1 = fig.add_subplot(121) # Plate Carrée projection - ax2 = fig.add_subplot(122, projection="3d") # 3D view + ax1 = fig.add_subplot(131) # Plate Carrée projection + ax2 = fig.add_subplot(132, projection="polar") # Top view + ax3 = fig.add_subplot(133, projection="3d") # 3D view # 1. Plate Carrée projection using advanced plot function plot_spherical_small_circles_advanced( @@ -411,26 +424,38 @@ def plot_spherical_small_circles_multiplot( ax=ax1, ) # Customize for multiplot - ax1.set_xlabel("Azimuth (radians)", fontsize=10) - ax1.set_ylabel("Latitude (radians)", fontsize=10) + ax1.set_xlabel("Azimuth (deg)", fontsize=10) + ax1.set_ylabel("Latitude (deg)", fontsize=10) ax1.set_title("Plate Carrée Projection", fontsize=12) - # 2. 3D view + # 2. Top view (looking down z-axis) + plot_spherical_small_circles_polar( + circles=circles, + num_points=num_points, + title="Top View (Down Z-Axis)", + ax=ax2, + ) + # Customize for multiplot + ax2.set_xlabel("Azimuth (deg)", fontsize=10) + ax2.set_ylabel("Latitude (deg)", fontsize=10) + ax2.set_title("Top View (Down Z-Axis)", fontsize=12) + + # 3. 3D view plot_spherical_small_circles_3d( circles=circles, num_points=num_points, title="3D View", show_sphere=True, show_legend=True, - ax=ax2, + ax=ax3, ) # Customize for multiplot - ax2.set_xlabel("X", fontsize=10) - ax2.set_ylabel("Y", fontsize=10) - ax2.set_zlabel("Z", fontsize=10) - ax2.set_title("3D View", fontsize=12) + ax3.set_xlabel("X", fontsize=10) + ax3.set_ylabel("Y", fontsize=10) + ax3.set_zlabel("Z", fontsize=10) + ax3.set_title("3D View", fontsize=12) # Position legend outside plot area with smaller font - ax2.legend(loc="upper left", bbox_to_anchor=(1.05, 1.0), fontsize=8) + ax3.legend(loc="upper left", bbox_to_anchor=(1.05, 1.0), fontsize=8) plt.tight_layout() @@ -525,11 +550,6 @@ def demo() -> None: show_plot=True, ) - _ = plot_spherical_small_circles_polar( - circles=circles, - show_plot=True, - ) - if __name__ == "__main__": demo() From ed29d25de4c13957bad883bfd13af30f656feee4 Mon Sep 17 00:00:00 2001 From: kopecn Date: Mon, 1 Dec 2025 20:42:50 -0800 Subject: [PATCH 11/22] neat... --- .../plotUnitSphericalSmallCircles.py | 32 +++++++++++++++++-- 1 file changed, 30 insertions(+), 2 deletions(-) diff --git a/src/pyMathToolsPlotHelpers/plotUnitSphericalSmallCircles.py b/src/pyMathToolsPlotHelpers/plotUnitSphericalSmallCircles.py index 436bd36..b9cb670 100644 --- a/src/pyMathToolsPlotHelpers/plotUnitSphericalSmallCircles.py +++ b/src/pyMathToolsPlotHelpers/plotUnitSphericalSmallCircles.py @@ -251,14 +251,42 @@ def plot_spherical_small_circles_polar( # polar ranges from -pi/2 to pi/2, map to 0 to 1 radius_points = (polar_points + np.pi / 2) / np.pi + # Handle azimuth wraparound by detecting large jumps in the ORIGINAL order + # Insert NaN values at discontinuities to break the line + azimuth_plot = [azimuth_points[0]] + radius_plot = [radius_points[0]] + + for i in range(1, len(azimuth_points)): + # If there's a jump greater than π, we've crossed the -π to π discontinuity + if abs(azimuth_points[i] - azimuth_points[i-1]) > np.pi: + # Insert NaN to break the line + azimuth_plot.append(np.nan) + radius_plot.append(np.nan) + + azimuth_plot.append(azimuth_points[i]) + radius_plot.append(radius_points[i]) + + # Close the circle by connecting back to the first point + # Check if we need to handle wraparound at the end + if abs(azimuth_points[0] - azimuth_points[-1]) > np.pi: + # Large jump at the end - don't connect + pass + else: + # Normal case - connect back to start + azimuth_plot.append(azimuth_points[0]) + radius_plot.append(radius_points[0]) + + azimuth_plot = np.array(azimuth_plot) + radius_plot = np.array(radius_plot) + # Convert radius to degrees for legend radius_deg = np.degrees(circle.radius_angle) azimuth_deg = np.degrees(circle.azimuth) polar_deg = np.degrees(circle.polar) label = f"Circle {idx + 1}: A:{azimuth_deg:.1f}°, P:{polar_deg:.1f}°, R:{radius_deg:.1f}°" - # Plot in polar coordinates - ax.scatter(azimuth_points, radius_points, s=10, alpha=0.6, label=label) + # Plot in polar coordinates using plot instead of scatter for continuous lines + ax.plot(azimuth_plot, radius_plot, linewidth=2, alpha=0.7, label=label) ax.set_title(title, fontsize=14, pad=20) From 5f6fdbdd3ca6e9b79acee37849ca96bbb6b80296 Mon Sep 17 00:00:00 2001 From: kopecn Date: Tue, 2 Dec 2025 00:32:52 -0800 Subject: [PATCH 12/22] neat... --- .../generators/sphericalGenerators.py | 123 +++- .../plotUnitSphericalSmallCircles.py | 539 ++++++++++++------ 2 files changed, 466 insertions(+), 196 deletions(-) diff --git a/src/pyMathTools/generators/sphericalGenerators.py b/src/pyMathTools/generators/sphericalGenerators.py index b50e67a..099534c 100644 --- a/src/pyMathTools/generators/sphericalGenerators.py +++ b/src/pyMathTools/generators/sphericalGenerators.py @@ -15,6 +15,7 @@ ) from numpy.linalg import norm from foundationTypes.mathTypes.UnitSphericalSmallCircle import UnitSphericalSmallCircle +from foundationTypes.mathTypes.UnitSphericalArc import UnitSphericalArc from pyMathTools.hints import FloatNDArray @@ -33,9 +34,10 @@ def generate_spherical_small_circle_points( ----------- circle.azimuth : float or NDArray Azimuth angle of circle center (radians, -2π to 2π) - polar_center : float or NDArray - Polar angle of circle center (radians, 0 to π) - radius_angle : float or NDArray + circle.polar : float or NDArray + Polar angle of circle center (radians, physics convention: 0 to π) + where 0 = north pole, π/2 = equator, π = south pole + circle.radius_angle : float or NDArray Angular radius of the circle (radians) num_points : int Number of points to generate around the circle @@ -43,12 +45,13 @@ def generate_spherical_small_circle_points( Returns: -------- azimuth_points : NDArray - Azimuth angles of points on the circle + Azimuth angles of points on the circle (radians) polar_points : NDArray - Polar angles of points on the circle + Polar angles of points on the circle (physics convention: 0 to π, colatitude) """ + # Use physics convention directly: polar ∈ [0, π], where 0 = north pole, π/2 = equator # Convert center to Cartesian coordinates (physics convention) - # theta (polar) = 0 at north pole, pi at south pole + # theta (polar/colatitude) = 0 at north pole, π at south pole # phi (azimuth) = rotation about Z axis x_center = sin(circle.polar) * cos(circle.azimuth) y_center = sin(circle.polar) * sin(circle.azimuth) @@ -107,7 +110,111 @@ def generate_spherical_small_circle_points( r = sqrt(x**2 + y**2 + z**2) x, y, z = x / r, y / r, z / r # Normalize to unit sphere - polar_points[i] = arccos(clip(z, -1, 1)) - azimuth_points[i] = arctan2(y, x) + # Output in physics convention for compatibility with downstream code + polar_points[i] = arccos(clip(z, -1, 1)) # colatitude (0 to π) + azimuth_points[i] = arctan2(y, x) # azimuth + + return azimuth_points, polar_points + + +def generate_spherical_arc_points( + arc: UnitSphericalArc, + num_points: int = 120, +) -> Tuple[FloatNDArray, FloatNDArray]: + """ + Generate points along an arc on a unit sphere. + + The arc starts at the point (azimuth, polar) and extends along a great circle + in the direction determined by the orient parameter for the specified arc_length. + + The orientation follows the right-hand rule: point your right thumb along the radial + vector (from origin through the start point), and your fingers curl in the positive + direction of rotation (positive orient). + + Parameters: + ----------- + arc : UnitSphericalArc + Arc specification with: + - azimuth: starting point longitude (radians) + - polar: starting point colatitude in physics convention (0 to π, 0 = north pole) + - orient: rotation angle about the radial vector (radians, -π to π) + Following right-hand rule: thumb along radial, fingers curl in +orient direction + - arc_length: extension along great circle (radians, can be negative) + num_points : int + Number of points to generate along the arc + + Returns: + -------- + azimuth_points : NDArray + Azimuth angles of points along the arc (radians) + polar_points : NDArray + Polar angles of points along the arc (physics convention: 0 to π, colatitude) + """ + # Use physics convention directly: polar ∈ [0, π], where 0 = north pole, π/2 = equator + # Convert starting point to Cartesian coordinates (physics convention) + x_start = sin(arc.polar) * cos(arc.azimuth) + y_start = sin(arc.polar) * sin(arc.azimuth) + z_start = cos(arc.polar) + start_point = array([x_start, y_start, z_start]) + + # Create tangent direction based on the orient parameter + # The orient defines the rotation angle about the radial vector (right-hand rule) + # We need to create an orthonormal basis in the tangent plane and rotate by orient + + north = array([0, 0, 1]) + + # Create a reference direction in the tangent plane + # Use the north direction or an arbitrary perpendicular if at poles + if abs(z_start) < 0.999: # Not at poles + # Reference direction: perpendicular to both north and radial + ref_tangent = cross(north, start_point) + ref_tangent = ref_tangent / norm(ref_tangent) + else: # At poles, use arbitrary reference + ref_tangent = array([1, 0, 0]) + + # Second tangent direction (perpendicular to both radial and reference) + second_tangent = cross(start_point, ref_tangent) + second_tangent = second_tangent / norm(second_tangent) + + # Rotate the reference tangent by the orient angle about the radial vector + # Using the rotation formula in the tangent plane + initial_direction = cos(arc.orient) * ref_tangent + sin(arc.orient) * second_tangent + initial_direction = initial_direction / norm(initial_direction) + + # The rotation axis for the arc is perpendicular to the initial direction and radial + rotation_axis = cross(start_point, initial_direction) + rotation_axis = rotation_axis / norm(rotation_axis) + + # Generate points along the arc + # Arc length parameter: 0 to arc_length + t = linspace(0, arc.arc_length, num_points) + + azimuth_points = zeros(num_points) + polar_points = zeros(num_points) + + for i, arc_dist in enumerate(t): + # Move along great circle in the direction determined by orient + # Using Rodrigues' rotation formula + # Rotation axis is perpendicular to initial direction (follows from right-hand rule) + axis = rotation_axis + + # Rotate start_point around axis by arc_dist + cos_angle = cos(arc_dist) + sin_angle = sin(arc_dist) + + point = ( + cos_angle * start_point + + sin_angle * cross(axis, start_point) + + (1 - cos_angle) * dot(axis, start_point) * axis + ) + + # Convert back to spherical coordinates (physics convention) + x, y, z = point + r = sqrt(x**2 + y**2 + z**2) + x, y, z = x / r, y / r, z / r # Normalize to unit sphere + + # Output in physics convention for compatibility with downstream code + polar_points[i] = arccos(clip(z, -1, 1)) # colatitude (0 to π) + azimuth_points[i] = arctan2(y, x) # azimuth return azimuth_points, polar_points diff --git a/src/pyMathToolsPlotHelpers/plotUnitSphericalSmallCircles.py b/src/pyMathToolsPlotHelpers/plotUnitSphericalSmallCircles.py index b9cb670..da4535e 100644 --- a/src/pyMathToolsPlotHelpers/plotUnitSphericalSmallCircles.py +++ b/src/pyMathToolsPlotHelpers/plotUnitSphericalSmallCircles.py @@ -1,19 +1,27 @@ -from typing import List, Union, Tuple, Optional +from typing import List, Tuple, Optional import numpy as np import matplotlib.pyplot as plt from matplotlib.figure import Figure from matplotlib.axes import Axes from foundationTypes.mathTypes.UnitSphericalSmallCircle import UnitSphericalSmallCircle +from foundationTypes.mathTypes.UnitSphericalArc import UnitSphericalArc from pyMathTools.transforms.sphericalTransforms import plate_carree_transform from pyMathTools.generators.sphericalGenerators import ( generate_spherical_small_circle_points, + generate_spherical_arc_points, ) +deg45 = np.deg2rad(45) +deg90 = np.deg2rad(90) +deg180 = np.deg2rad(180) +deg22_5 = np.pi / 8 -def plot_spherical_small_circles( - circles: List[UnitSphericalSmallCircle], + +def plot_unit_spherical( + circles: List[UnitSphericalSmallCircle] | None = None, + arcs: List[UnitSphericalArc] | None = None, num_points: int = 120, figsize: Tuple[int, int] = (10, 6), title: str = "Small Circles on Sphere (Plate Carrée Projection)", @@ -21,21 +29,27 @@ def plot_spherical_small_circles( show_plot: bool = False, ) -> Tuple[Figure, Axes]: """ - Plot small circles on a sphere using plate carrée (equirectangular) projection. + Plot small circles and arcs on a sphere using plate carrée (equirectangular) projection. Parameters: ----------- - circles : List[UnitSphericalSmallCircle] or array-like - List of UnitSphericalSmallCircle objects or array of [azimuth, polar, radius] lists. + circles : List[UnitSphericalSmallCircle], optional + List of UnitSphericalSmallCircle objects. Each circle represents a small circle on the unit sphere + arcs : List[UnitSphericalArc], optional + List of UnitSphericalArc objects. + Each arc starts at (azimuth, polar) and extends along a great circle + for the specified arc_length (positive or negative) num_points : int, optional - Number of points to generate around each circle (default: 120) + Number of points to generate for each circle/arc (default: 120) figsize : tuple, optional Figure size (default: (10, 6)) title : str, optional Plot title show_legend : bool, optional Whether to show the legend (default: True) + show_plot : bool, optional + Whether to display the plot (default: False) Returns: -------- @@ -44,30 +58,60 @@ def plot_spherical_small_circles( fig, ax = plt.subplots(figsize=figsize) # Generate points for each circle - for idx, circle in enumerate(circles): - - # Generate true spherical small circle points - azimuth_points, polar_points = generate_spherical_small_circle_points( - circle, num_points - ) - - # Apply Plate Carrée projection - x_points, y_points = plate_carree_transform(azimuth_points, polar_points) - - # Convert radius to degrees for legend - radius_deg = np.degrees(circle.radius_angle) - azimuth_deg = np.degrees(circle.azimuth) - polar_deg = np.degrees(circle.polar) - label = f"Circle {idx + 1}: A:{azimuth_deg:.1f}°, P:{polar_deg:.1f}°, R:{radius_deg:.1f}°" - - # Plot the circle - ax.scatter( - np.degrees(x_points), - np.degrees(y_points), - s=10, - alpha=0.6, - label=label, - ) + if circles is not None: + for idx, circle in enumerate(circles): + + # Generate true spherical small circle points + azimuth_points, polar_points = generate_spherical_small_circle_points( + circle, num_points + ) + + # Apply Plate Carrée projection + x_points, y_points = plate_carree_transform(azimuth_points, polar_points) + + # Convert radius to degrees for legend + radius_deg = np.degrees(circle.radius_angle) + azimuth_deg = np.degrees(circle.azimuth) + polar_deg = np.degrees(circle.polar) + label = f"Circle {idx + 1}: A:{azimuth_deg:.1f}°, P:{polar_deg:.1f}°, R:{radius_deg:.1f}°" + + # Plot the circle + ax.scatter( + np.degrees(x_points), + np.degrees(y_points), + s=10, + alpha=0.6, + label=label, + ) + + # Generate points for each arc + if arcs is not None: + for idx, arc in enumerate(arcs): + + # Generate spherical arc points + azimuth_points, polar_points = generate_spherical_arc_points( + arc, num_points + ) + + # Apply Plate Carrée projection + x_points, y_points = plate_carree_transform(azimuth_points, polar_points) + + # Convert to degrees for legend + arc_length_deg = np.degrees(arc.arc_length) + azimuth_deg = np.degrees(arc.azimuth) + polar_deg = np.degrees(arc.polar) + label = f"Arc {idx + 1}: A:{azimuth_deg:.1f}°, P:{polar_deg:.1f}°, L:{arc_length_deg:.1f}°" + + # Plot the arc as a line + ax.plot( + np.degrees(x_points), + np.degrees(y_points), + linewidth=2, + alpha=0.8, + label=label, + marker="o", + markersize=3, + ) # Set up axes ax.set_xlabel("Azimuth (deg)", fontsize=12) @@ -95,8 +139,9 @@ def plot_spherical_small_circles( return fig, ax -def plot_spherical_small_circles_advanced( - circles: List[UnitSphericalSmallCircle], +def plot_unit_spherical_advanced( + circles: List[UnitSphericalSmallCircle] | None = None, + arcs: List[UnitSphericalArc] | None = None, num_points: int = 120, figsize: Tuple[int, int] = (10, 6), title: str = "Small Circles on Sphere (Flattened)", @@ -113,6 +158,10 @@ def plot_spherical_small_circles_advanced( ----------- circles : List[UnitSphericalSmallCircle] or array-like List of UnitSphericalSmallCircle objects or array of [azimuth, polar, radius] lists + arcs : List[UnitSphericalArc], optional + List of UnitSphericalArc objects. + Each arc starts at (azimuth, polar) and extends along a great circle + for the specified arc_length (positive or negative) num_points : int, optional Number of points to generate around each circle figsize : tuple, optional @@ -142,40 +191,71 @@ def plot_spherical_small_circles_advanced( else: fig = ax.get_figure() - for idx, circle in enumerate(circles): - - # Generate true spherical small circle points - azimuth_points, polar_points = generate_spherical_small_circle_points( - circle, num_points - ) - - # Apply Plate Carrée projection - x_points, y_points = plate_carree_transform(azimuth_points, polar_points) - - # Determine color and label - color = colors[idx] if colors is not None and idx < len(colors) else None - label = labels[idx] if labels is not None and idx < len(labels) else None - - # Plot the circle - ax.scatter( - np.degrees(x_points), - np.degrees(y_points), - s=10, - alpha=0.6, - c=color, - label=label, - ) - - # Plot center point if requested - if show_centers: - x_center, y_center = plate_carree_transform(circle.azimuth, circle.polar) + if circles is not None: + for idx, circle in enumerate(circles): + + # Generate true spherical small circle points + azimuth_points, polar_points = generate_spherical_small_circle_points( + circle, num_points + ) + + # Apply Plate Carrée projection + x_points, y_points = plate_carree_transform(azimuth_points, polar_points) + + # Determine color and label + color = colors[idx] if colors is not None and idx < len(colors) else None + label = labels[idx] if labels is not None and idx < len(labels) else None + + # Plot the circle ax.scatter( - np.degrees(x_center), - np.degrees(y_center), - s=100, - marker="x", - c=color if color else "red", - linewidths=2, + np.degrees(x_points), + np.degrees(y_points), + s=10, + alpha=0.6, + c=color, + label=label, + ) + + # Plot center point if requested + if show_centers: + x_center, y_center = plate_carree_transform( + circle.azimuth, circle.polar + ) + ax.scatter( + np.degrees(x_center), + np.degrees(y_center), + s=100, + marker="x", + c=color if color else "red", + linewidths=2, + ) + + # Generate points for each arc + if arcs is not None: + for idx, arc in enumerate(arcs): + + # Generate spherical arc points + azimuth_points, polar_points = generate_spherical_arc_points( + arc, num_points + ) + + # Apply Plate Carrée projection + x_points, y_points = plate_carree_transform(azimuth_points, polar_points) + + # Determine color and label + color = colors[idx] if colors is not None and idx < len(colors) else None + label = labels[idx] if labels is not None and idx < len(labels) else None + + # Plot the arc as a line + ax.plot( + np.degrees(x_points), + np.degrees(y_points), + linewidth=2, + alpha=0.8, + c=color, + label=label, + marker="o", + markersize=3, ) # Set up axes @@ -204,8 +284,9 @@ def plot_spherical_small_circles_advanced( return fig, ax -def plot_spherical_small_circles_polar( - circles: List[UnitSphericalSmallCircle], +def plot_unit_spherical_polar( + circles: List[UnitSphericalSmallCircle] | None = None, + arcs: List[UnitSphericalArc] | None = None, num_points: int = 120, figsize: Tuple[int, int] = (8, 8), title: str = "Small Circles on Sphere (Polar View)", @@ -214,13 +295,17 @@ def plot_spherical_small_circles_polar( show_plot: bool = False, ) -> Tuple[Figure, Axes]: """ - Plot small circles on a sphere using polar projection. + Plot small circles and arcs on a sphere using polar projection. Azimuth is shown as the angle, polar angle is mapped to radius. Parameters: ----------- circles : List[UnitSphericalSmallCircle] or array-like List of UnitSphericalSmallCircle objects or array of [azimuth, polar, radius] lists + arcs : List[UnitSphericalArc], optional + List of UnitSphericalArc objects. + Each arc starts at (azimuth, polar) and extends along a great circle + for the specified arc_length (positive or negative) num_points : int, optional Number of points to generate around each circle (default: 120) figsize : tuple, optional @@ -240,53 +325,100 @@ def plot_spherical_small_circles_polar( else: fig = ax.get_figure() - for idx, circle in enumerate(circles): - - # Generate true spherical small circle points - azimuth_points, polar_points = generate_spherical_small_circle_points( - circle, num_points - ) - - # Map polar angle to radius (0 to 1) for polar plot - # polar ranges from -pi/2 to pi/2, map to 0 to 1 - radius_points = (polar_points + np.pi / 2) / np.pi - - # Handle azimuth wraparound by detecting large jumps in the ORIGINAL order - # Insert NaN values at discontinuities to break the line - azimuth_plot = [azimuth_points[0]] - radius_plot = [radius_points[0]] - - for i in range(1, len(azimuth_points)): - # If there's a jump greater than π, we've crossed the -π to π discontinuity - if abs(azimuth_points[i] - azimuth_points[i-1]) > np.pi: - # Insert NaN to break the line - azimuth_plot.append(np.nan) - radius_plot.append(np.nan) - - azimuth_plot.append(azimuth_points[i]) - radius_plot.append(radius_points[i]) - - # Close the circle by connecting back to the first point - # Check if we need to handle wraparound at the end - if abs(azimuth_points[0] - azimuth_points[-1]) > np.pi: - # Large jump at the end - don't connect - pass - else: - # Normal case - connect back to start - azimuth_plot.append(azimuth_points[0]) - radius_plot.append(radius_points[0]) - - azimuth_plot = np.array(azimuth_plot) - radius_plot = np.array(radius_plot) - - # Convert radius to degrees for legend - radius_deg = np.degrees(circle.radius_angle) - azimuth_deg = np.degrees(circle.azimuth) - polar_deg = np.degrees(circle.polar) - label = f"Circle {idx + 1}: A:{azimuth_deg:.1f}°, P:{polar_deg:.1f}°, R:{radius_deg:.1f}°" - - # Plot in polar coordinates using plot instead of scatter for continuous lines - ax.plot(azimuth_plot, radius_plot, linewidth=2, alpha=0.7, label=label) + if circles is not None: + for idx, circle in enumerate(circles): + + # Generate true spherical small circle points + azimuth_points, polar_points = generate_spherical_small_circle_points( + circle, num_points + ) + + # Map polar angle to radius (0 to 1) for polar plot + # polar ranges from -pi/2 to pi/2, map to 0 to 1 + radius_points = (polar_points + np.pi / 2) / np.pi + + # Handle azimuth wraparound by detecting large jumps in the ORIGINAL order + # Insert NaN values at discontinuities to break the line + azimuth_plot = [azimuth_points[0]] + radius_plot = [radius_points[0]] + + for i in range(1, len(azimuth_points)): + # If there's a jump greater than π, we've crossed the -π to π discontinuity + if abs(azimuth_points[i] - azimuth_points[i - 1]) > np.pi: + # Insert NaN to break the line + azimuth_plot.append(np.nan) + radius_plot.append(np.nan) + + azimuth_plot.append(azimuth_points[i]) + radius_plot.append(radius_points[i]) + + # Close the circle by connecting back to the first point + # Check if we need to handle wraparound at the end + if abs(azimuth_points[0] - azimuth_points[-1]) > np.pi: + # Large jump at the end - don't connect + pass + else: + # Normal case - connect back to start + azimuth_plot.append(azimuth_points[0]) + radius_plot.append(radius_points[0]) + + azimuth_plot = np.array(azimuth_plot) + radius_plot = np.array(radius_plot) + + # Convert radius to degrees for legend + radius_deg = np.degrees(circle.radius_angle) + azimuth_deg = np.degrees(circle.azimuth) + polar_deg = np.degrees(circle.polar) + label = f"Circle {idx + 1}: A:{azimuth_deg:.1f}°, P:{polar_deg:.1f}°, R:{radius_deg:.1f}°" + + # Plot in polar coordinates using plot instead of scatter for continuous lines + ax.plot(azimuth_plot, radius_plot, linewidth=2, alpha=0.7, label=label) + + # Generate points for each arc + if arcs is not None: + for idx, arc in enumerate(arcs): + + # Generate spherical arc points + azimuth_points, polar_points = generate_spherical_arc_points( + arc, num_points + ) + + # Map polar angle to radius (0 to 1) for polar plot + radius_points = (polar_points + np.pi / 2) / np.pi + + # Handle azimuth wraparound + azimuth_plot = [azimuth_points[0]] + radius_plot = [radius_points[0]] + + for i in range(1, len(azimuth_points)): + # If there's a jump greater than π, we've crossed the -π to π discontinuity + if abs(azimuth_points[i] - azimuth_points[i - 1]) > np.pi: + # Insert NaN to break the line + azimuth_plot.append(np.nan) + radius_plot.append(np.nan) + + azimuth_plot.append(azimuth_points[i]) + radius_plot.append(radius_points[i]) + + azimuth_plot = np.array(azimuth_plot) + radius_plot = np.array(radius_plot) + + # Convert to degrees for legend + arc_length_deg = np.degrees(arc.arc_length) + azimuth_deg = np.degrees(arc.azimuth) + polar_deg = np.degrees(arc.polar) + label = f"Arc {idx + 1}: A:{azimuth_deg:.1f}°, P:{polar_deg:.1f}°, L:{arc_length_deg:.1f}°" + + # Plot the arc + ax.plot( + azimuth_plot, + radius_plot, + linewidth=2, + alpha=0.8, + label=label, + marker="o", + markersize=3, + ) ax.set_title(title, fontsize=14, pad=20) @@ -304,8 +436,9 @@ def plot_spherical_small_circles_polar( return fig, ax -def plot_spherical_small_circles_3d( - circles: List[UnitSphericalSmallCircle], +def plot_unit_spherical_3d( + circles: List[UnitSphericalSmallCircle] | None = None, + arcs: List[UnitSphericalArc] | None = None, num_points: int = 120, figsize: Tuple[int, int] = (10, 10), title: str = "Small Circles on Unit Sphere (3D)", @@ -315,12 +448,16 @@ def plot_spherical_small_circles_3d( ax: Optional[Axes] = None, ) -> Tuple[Figure, Axes]: """ - Plot small circles on a 3D unit sphere. + Plot small circles and arcs on a 3D unit sphere. Parameters: ----------- circles : List[UnitSphericalSmallCircle] or array-like List of UnitSphericalSmallCircle objects or array of [azimuth, polar, radius] lists + arcs : List[UnitSphericalArc], optional + List of UnitSphericalArc objects. + Each arc starts at (azimuth, polar) and extends along a great circle + for the specified arc_length (positive or negative) num_points : int, optional Number of points to generate around each circle (default: 120) figsize : tuple, optional @@ -359,28 +496,61 @@ def plot_spherical_small_circles_3d( ax.plot_surface(x_sphere, y_sphere, z_sphere, alpha=0.1, color="lightblue") # Plot each circle - for idx, circle in enumerate(circles): - - # Generate true spherical small circle points - azimuth_points, polar_points = generate_spherical_small_circle_points( - circle, num_points - ) + if circles is not None: + for idx, circle in enumerate(circles): - # Convert to Cartesian coordinates (unit sphere, r=1, physics convention) - # theta (polar) from 0 (north pole) to pi (south pole) - # phi (azimuth) rotation about Z axis - x = np.sin(polar_points) * np.cos(azimuth_points) - y = np.sin(polar_points) * np.sin(azimuth_points) - z = np.cos(polar_points) + # Generate true spherical small circle points + azimuth_points, polar_points = generate_spherical_small_circle_points( + circle, num_points + ) - # Convert radius to degrees for legend - radius_deg = np.degrees(circle.radius_angle) - azimuth_deg = np.degrees(circle.azimuth) - polar_deg = np.degrees(circle.polar) - label = f"Circle {idx + 1}: A:{azimuth_deg:.1f}°, P:{polar_deg:.1f}°, R:{radius_deg:.1f}°" + # Convert to Cartesian coordinates (unit sphere, r=1, physics convention) + # theta (polar) from 0 (north pole) to pi (south pole) + # phi (azimuth) rotation about Z axis + x = np.sin(polar_points) * np.cos(azimuth_points) + y = np.sin(polar_points) * np.sin(azimuth_points) + z = np.cos(polar_points) + + # Convert radius to degrees for legend + radius_deg = np.degrees(circle.radius_angle) + azimuth_deg = np.degrees(circle.azimuth) + polar_deg = np.degrees(circle.polar) + label = f"Circle {idx + 1}: A:{azimuth_deg:.1f}°, P:{polar_deg:.1f}°, R:{radius_deg:.1f}°" + + # Plot the circle + ax.scatter(x, y, z, s=10, alpha=0.8, label=label) + + # Plot each arc + if arcs is not None: + for idx, arc in enumerate(arcs): + + # Generate spherical arc points + azimuth_points, polar_points = generate_spherical_arc_points( + arc, num_points + ) - # Plot the circle - ax.scatter(x, y, z, s=10, alpha=0.8, label=label) + # Convert to Cartesian coordinates (unit sphere, r=1, physics convention) + x = np.sin(polar_points) * np.cos(azimuth_points) + y = np.sin(polar_points) * np.sin(azimuth_points) + z = np.cos(polar_points) + + # Convert to degrees for legend + arc_length_deg = np.degrees(arc.arc_length) + azimuth_deg = np.degrees(arc.azimuth) + polar_deg = np.degrees(arc.polar) + label = f"Arc {idx + 1}: A:{azimuth_deg:.1f}°, P:{polar_deg:.1f}°, L:{arc_length_deg:.1f}°" + + # Plot the arc as a line + ax.plot( + x, + y, + z, + linewidth=2, + alpha=0.8, + label=label, + marker="o", + markersize=3, + ) ax.set_xlabel("X") ax.set_ylabel("Y") @@ -407,8 +577,9 @@ def plot_spherical_small_circles_3d( return fig, ax -def plot_spherical_small_circles_multiplot( - circles: Union[List[UnitSphericalSmallCircle], List[List[float]]], +def plot_unit_spherical_multiplot( + circles: List[UnitSphericalSmallCircle] | None = None, + arcs: List[UnitSphericalArc] | None = None, num_points: int = 120, figsize: Tuple[int, int] = (24, 6), title: str = "Small Circles on Unit Sphere - Multiple Views", @@ -444,8 +615,9 @@ def plot_spherical_small_circles_multiplot( ax3 = fig.add_subplot(133, projection="3d") # 3D view # 1. Plate Carrée projection using advanced plot function - plot_spherical_small_circles_advanced( + plot_unit_spherical_advanced( circles=circles, + arcs=arcs, num_points=num_points, title="Plate Carrée Projection", show_centers=False, @@ -457,8 +629,9 @@ def plot_spherical_small_circles_multiplot( ax1.set_title("Plate Carrée Projection", fontsize=12) # 2. Top view (looking down z-axis) - plot_spherical_small_circles_polar( + plot_unit_spherical_polar( circles=circles, + arcs=arcs, num_points=num_points, title="Top View (Down Z-Axis)", ax=ax2, @@ -469,8 +642,9 @@ def plot_spherical_small_circles_multiplot( ax2.set_title("Top View (Down Z-Axis)", fontsize=12) # 3. 3D view - plot_spherical_small_circles_3d( + plot_unit_spherical_3d( circles=circles, + arcs=arcs, num_points=num_points, title="3D View", show_sphere=True, @@ -501,80 +675,69 @@ def demo() -> None: UnitSphericalSmallCircle( azimuth=0, polar=0, - radius_angle=np.deg2rad(45), + radius_angle=deg45, ), UnitSphericalSmallCircle( azimuth=0, polar=np.pi / 2, - radius_angle=np.deg2rad(45), + radius_angle=deg45, ), UnitSphericalSmallCircle( azimuth=0, polar=np.pi, - radius_angle=np.deg2rad(45), + radius_angle=deg45, ), UnitSphericalSmallCircle( azimuth=0, polar=-np.pi / 2, - radius_angle=np.deg2rad(45), + radius_angle=deg45, ), UnitSphericalSmallCircle( azimuth=-np.pi / 2, polar=np.pi / 2, - radius_angle=np.deg2rad(45), + radius_angle=deg45, ), UnitSphericalSmallCircle( azimuth=np.pi / 2, polar=np.pi / 2, - radius_angle=np.deg2rad(45), + radius_angle=deg45, ), UnitSphericalSmallCircle( - azimuth=-np.deg2rad(45), - polar=np.deg2rad(45), - radius_angle=np.deg2rad(45), + azimuth=-deg45, + polar=deg45, + radius_angle=deg45, ), UnitSphericalSmallCircle( azimuth=np.deg2rad(135), - polar=np.deg2rad(45), - radius_angle=np.deg2rad(45), + polar=deg45, + radius_angle=deg45, ), - # UnitSphericalSmallCircle( # should overlap circle 1 - # azimuth=np.pi / 2, - # polar=0, - # radius_angle=np.deg2rad(45), - # ), - # UnitSphericalSmallCircle( - # azimuth=np.deg2rad(10), - # polar=np.deg2rad(-80), - # radius_angle=np.deg2rad(45), - # ), - # UnitSphericalSmallCircle( - # azimuth=np.deg2rad(10), - # polar=np.deg2rad(80), - # radius_angle=np.deg2rad(45), - # ), - # UnitSphericalSmallCircle( - # azimuth=np.pi, # 180° - # polar=np.pi / 4, # 45° north - # radius_angle=0.2, - # ), - # UnitSphericalSmallCircle( - # azimuth=3 * np.pi / 2, # 270° - # polar=-np.pi / 6, # -30° south - # radius_angle=0.25, - # ), ] - # Create the multiplot showing all three views - _ = plot_spherical_small_circles_multiplot( - circles, - num_points=240, - title="Unit Spherical Small Circles - Plate Carrée, Gauss-Krüger, and 3D Views", - show_plot=True, - ) + # # Create the multiplot showing all three views + # _ = plot_unit_spherical_multiplot( + # circles=circles, + # num_points=240, + # title="Unit Spherical Small Circles - Plate Carrée, Gauss-Krüger, and 3D Views", + # show_plot=True, + # ) + + # _ = plot_unit_spherical( + # circles=circles, + # show_plot=True, + # ) + + arcs: List[UnitSphericalArc] = [ + UnitSphericalArc(orient=0, azimuth=0, polar=0, arc_length=deg22_5), + UnitSphericalArc(orient=deg45, azimuth=0, polar=0, arc_length=deg22_5), + UnitSphericalArc(orient=deg90, azimuth=0, polar=0, arc_length=deg22_5), + UnitSphericalArc(orient=0, azimuth=deg45, polar=deg45, arc_length=deg22_5), + UnitSphericalArc(orient=deg45, azimuth=deg45, polar=deg45, arc_length=deg22_5), + UnitSphericalArc(orient=deg90, azimuth=deg45, polar=deg45, arc_length=deg22_5), + ] - _ = plot_spherical_small_circles( - circles=circles, + _ = plot_unit_spherical_multiplot( + arcs=arcs, show_plot=True, ) From 4ed3a85a62c6adc8aecb13816a999c7d1225c037 Mon Sep 17 00:00:00 2001 From: kopecn Date: Tue, 2 Dec 2025 11:56:02 -0800 Subject: [PATCH 13/22] neat... --- examples/plotQuatUnitCircles.py | 8 +++--- .../generators/sphericalGenerators.py | 28 +++++++++++-------- .../plotUnitSphericalSmallCircles.py | 23 +++++++-------- 3 files changed, 31 insertions(+), 28 deletions(-) diff --git a/examples/plotQuatUnitCircles.py b/examples/plotQuatUnitCircles.py index e5905c1..5ac5ea2 100644 --- a/examples/plotQuatUnitCircles.py +++ b/examples/plotQuatUnitCircles.py @@ -10,8 +10,8 @@ 1. Import required modules and functions. 2. Convert several canonical quaternions (identity and axis-aligned) to UnitSphericalSmallCircle objects. -3. Plot these circles using the multiplot function, which shows Mercator, - Gauss-Krüger, and 3D views. +3. Plot these circles using the multiplot function, which shows flat projection, + polar, and 3D views. Usage: ------ @@ -33,7 +33,7 @@ ) from pyMathToolsPlotHelpers.plotUnitSphericalSmallCircles import ( - plot_spherical_small_circles_multiplot, + plot_unit_spherical_multiplot, ) @@ -48,7 +48,7 @@ def main(): All circles use a radius angle of pi/8 except y and z, which use the default. """ - _ = plot_spherical_small_circles_multiplot( + _ = plot_unit_spherical_multiplot( circles=[ quat_to_unitSphericalSmallCircle(q=quaternion.one, radius_angle=pi / 8), quat_to_unitSphericalSmallCircle(q=quaternion.x, radius_angle=pi / 8), diff --git a/src/pyMathTools/generators/sphericalGenerators.py b/src/pyMathTools/generators/sphericalGenerators.py index 099534c..b01be06 100644 --- a/src/pyMathTools/generators/sphericalGenerators.py +++ b/src/pyMathTools/generators/sphericalGenerators.py @@ -160,24 +160,30 @@ def generate_spherical_arc_points( # Create tangent direction based on the orient parameter # The orient defines the rotation angle about the radial vector (right-hand rule) # We need to create an orthonormal basis in the tangent plane and rotate by orient + # Convention: orient=0 points towards the north pole in the tangent plane north = array([0, 0, 1]) - # Create a reference direction in the tangent plane - # Use the north direction or an arbitrary perpendicular if at poles + # Create an orthonormal basis in the tangent plane + # Following right-hand rule: thumb along radial, fingers curl in positive orient direction if abs(z_start) < 0.999: # Not at poles - # Reference direction: perpendicular to both north and radial - ref_tangent = cross(north, start_point) - ref_tangent = ref_tangent / norm(ref_tangent) + # East direction: perpendicular to both north and radial + east = cross(north, start_point) + east = east / norm(east) + + # North direction in tangent plane: completes right-handed system + # cross(north_tangent, east) = start_point (outward radial) + north_tangent = cross(east, start_point) + north_tangent = north_tangent / norm(north_tangent) + + ref_tangent = north_tangent + second_tangent = east else: # At poles, use arbitrary reference ref_tangent = array([1, 0, 0]) + second_tangent = array([0, 1, 0]) - # Second tangent direction (perpendicular to both radial and reference) - second_tangent = cross(start_point, ref_tangent) - second_tangent = second_tangent / norm(second_tangent) - - # Rotate the reference tangent by the orient angle about the radial vector - # Using the rotation formula in the tangent plane + # Rotate by the orient angle about the radial vector (right-hand rule) + # orient=0 points north, positive orient rotates towards east initial_direction = cos(arc.orient) * ref_tangent + sin(arc.orient) * second_tangent initial_direction = initial_direction / norm(initial_direction) diff --git a/src/pyMathToolsPlotHelpers/plotUnitSphericalSmallCircles.py b/src/pyMathToolsPlotHelpers/plotUnitSphericalSmallCircles.py index da4535e..8dedb97 100644 --- a/src/pyMathToolsPlotHelpers/plotUnitSphericalSmallCircles.py +++ b/src/pyMathToolsPlotHelpers/plotUnitSphericalSmallCircles.py @@ -714,29 +714,26 @@ def demo() -> None: ), ] - # # Create the multiplot showing all three views - # _ = plot_unit_spherical_multiplot( - # circles=circles, - # num_points=240, - # title="Unit Spherical Small Circles - Plate Carrée, Gauss-Krüger, and 3D Views", - # show_plot=True, - # ) - - # _ = plot_unit_spherical( - # circles=circles, - # show_plot=True, - # ) - arcs: List[UnitSphericalArc] = [ UnitSphericalArc(orient=0, azimuth=0, polar=0, arc_length=deg22_5), UnitSphericalArc(orient=deg45, azimuth=0, polar=0, arc_length=deg22_5), UnitSphericalArc(orient=deg90, azimuth=0, polar=0, arc_length=deg22_5), + UnitSphericalArc(orient=0, azimuth=0, polar=deg45, arc_length=deg22_5), UnitSphericalArc(orient=0, azimuth=deg45, polar=deg45, arc_length=deg22_5), + UnitSphericalArc(orient=0, azimuth=deg90, polar=deg45, arc_length=deg22_5), + UnitSphericalArc(orient=deg45, azimuth=0, polar=deg45, arc_length=deg22_5), UnitSphericalArc(orient=deg45, azimuth=deg45, polar=deg45, arc_length=deg22_5), + UnitSphericalArc(orient=deg45, azimuth=deg90, polar=deg45, arc_length=deg22_5), + UnitSphericalArc(orient=deg90, azimuth=0, polar=deg45, arc_length=deg22_5), UnitSphericalArc(orient=deg90, azimuth=deg45, polar=deg45, arc_length=deg22_5), + UnitSphericalArc(orient=deg90, azimuth=deg90, polar=deg45, arc_length=deg22_5), + UnitSphericalArc(orient=0, azimuth=0, polar=deg45, arc_length=-deg22_5), + UnitSphericalArc(orient=0, azimuth=deg45, polar=deg45, arc_length=-deg22_5), + UnitSphericalArc(orient=0, azimuth=deg90, polar=deg45, arc_length=-deg22_5), ] _ = plot_unit_spherical_multiplot( + circles=circles, arcs=arcs, show_plot=True, ) From 2860e2ce6942fec5ec0b43e299a014ba8e218ecd Mon Sep 17 00:00:00 2001 From: kopecn Date: Tue, 2 Dec 2025 14:06:56 -0800 Subject: [PATCH 14/22] neat... --- .claude/settings.local.json | 9 + .../constructUnitSphericalSmallCircle.py | 9 +- src/pyMathTools/hints.py | 21 +- src/pyMathTools/spatial/Quaternion.py | 577 ++++++++++++ src/pyMathTools/spatial/__init__.py | 0 tests/test_quaternion.py | 846 ++++++++++++++++++ tests/testpyMathTools.py | 11 - 7 files changed, 1456 insertions(+), 17 deletions(-) create mode 100644 .claude/settings.local.json create mode 100644 src/pyMathTools/spatial/Quaternion.py create mode 100644 src/pyMathTools/spatial/__init__.py create mode 100644 tests/test_quaternion.py delete mode 100644 tests/testpyMathTools.py diff --git a/.claude/settings.local.json b/.claude/settings.local.json new file mode 100644 index 0000000..f315077 --- /dev/null +++ b/.claude/settings.local.json @@ -0,0 +1,9 @@ +{ + "permissions": { + "allow": [ + "Bash(python:*)" + ], + "deny": [], + "ask": [] + } +} diff --git a/src/pyMathTools/constructors/constructUnitSphericalSmallCircle.py b/src/pyMathTools/constructors/constructUnitSphericalSmallCircle.py index 6c6d167..5cbffab 100644 --- a/src/pyMathTools/constructors/constructUnitSphericalSmallCircle.py +++ b/src/pyMathTools/constructors/constructUnitSphericalSmallCircle.py @@ -1,12 +1,13 @@ import numpy as np -from numpy import atan2, acos -import quaternion -from numpy import pi +from numpy import atan2, acos, pi + from foundationTypes.mathTypes.UnitSphericalSmallCircle import UnitSphericalSmallCircle +from pyMathTools.spatial.Quaternion import Quaternion + def quat_to_unitSphericalSmallCircle( - q: np.quaternion, radius_angle: float = pi / 4 + q: Quaternion, radius_angle: float = pi / 4 ) -> UnitSphericalSmallCircle: # Extract vector part diff --git a/src/pyMathTools/hints.py b/src/pyMathTools/hints.py index 1033936..2305056 100644 --- a/src/pyMathTools/hints.py +++ b/src/pyMathTools/hints.py @@ -5,10 +5,27 @@ enabling strict and clear typing in mathematical code throughout the package. """ -from typing import Union +from typing import Union, TYPE_CHECKING +import numpy as np from numpy import float64 from numpy.typing import NDArray +import quaternion -FloatOrNDArray = Union[float, NDArray[float64]] +if TYPE_CHECKING: + from pyMathTools.spatial.quaternion import Quaternion +# Basic numeric types +FloatOrNDArray = Union[float, NDArray[float64]] FloatNDArray = NDArray[float64] + +# Quaternion types +QuaternionLike = Union["Quaternion", np.quaternion, NDArray[np.quaternion]] +FloatOrQuaternion = Union["Quaternion", float, np.quaternion] +FloatArray3 = NDArray[float64] # Shape (..., 3) +FloatArray4 = NDArray[float64] # Shape (..., 4) +RotationMatrix = NDArray[float64] # Shape (..., 3, 3) + + +if __name__ == "__main__": + q = quaternion.one + assert isinstance(q, np.quaternion) diff --git a/src/pyMathTools/spatial/Quaternion.py b/src/pyMathTools/spatial/Quaternion.py new file mode 100644 index 0000000..0294127 --- /dev/null +++ b/src/pyMathTools/spatial/Quaternion.py @@ -0,0 +1,577 @@ +""" +Fully-typed OOP wrapper around numpy-quaternions library +""" + +from __future__ import annotations +from dataclasses import dataclass +from typing import TypeVar, Type, Any +import numpy as np +from foundationTypes.mathTypes.QuaternionType import QuaternionType + +# Import all necessary functions from quaternion module +from quaternion import ( + as_float_array, + as_quat_array, + from_float_array, + as_rotation_matrix, + from_rotation_matrix, + as_rotation_vector, + from_rotation_vector, + as_euler_angles, + from_euler_angles, + as_vector_part, + from_vector_part, + rotate_vectors, + isclose as quat_isclose, + allclose as quat_allclose, +) +from quaternion.quaternion_time_series import slerp as quat_slerp + +# Import custom type hints +from pyMathTools.hints import ( + FloatArray3, + FloatArray4, + RotationMatrix, + FloatOrQuaternion, +) + +T = TypeVar("T", bound="Quaternion") + + +def from_float(x: Any) -> float: + assert isinstance(x, (float, int)) and not isinstance(x, bool) + return float(x) + + +@dataclass +class Quaternion(QuaternionType): + """ + wrapper class for numpy-quaternion for allowing + """ + + __q: np.quaternion + + # MARK: - Constructors + + @staticmethod + def fromNumpyQuaternion(q: np.quaternion) -> "Quaternion": + return Quaternion(q) + + @classmethod + def from_components( + cls: Type[T], w: float = 1, x: float = 0, y: float = 0, z: float = 0 + ) -> T: + """Create a quaternion instance from individual w, x, y, z components. + + Args: + w: The scalar (real) component + x: The x component of the vector part + y: The y component of the vector part + z: The z component of the vector part + + Returns: + A concrete QuaternionType instance of the calling class type + """ + return cls(np.quaternion(w, x, y, z)) + + @classmethod + def from_float_array(cls: Type[T], array: FloatArray4) -> T: + """Create a quaternion from a 4-element float array [w, x, y, z]. + + Args: + array: 4-element array in [w, x, y, z] order + + Returns: + A concrete QuaternionType instance of the calling class type + + Raises: + ValueError: If array doesn't have exactly 4 elements + """ + q = from_float_array(array) + return cls(q) + + @staticmethod + def as_quat_array(array: np.ndarray) -> np.ndarray: + """View a float array as an array of quaternions. + + The input array must have a final dimension whose size is 4. + Each set of 4 components will be interpreted as [w, x, y, z]. + + This is a fast operation (no data copy) when the array is C-contiguous. + + Args: + array: Float array with last dimension of size 4 + Shape can be (..., 4) where ... represents any number of dimensions + Examples: + - [w, x, y, z] -> single quaternion + - [[w1,x1,y1,z1], [w2,x2,y2,z2]] -> array of 2 quaternions + + Returns: + NumPy array of quaternions with shape (...) + + Examples: + >>> # Single quaternion + >>> Quaternion.as_quat_array([1, 0, 0, 0]) + quaternion(1, 0, 0, 0) + + >>> # Array of quaternions + >>> Quaternion.as_quat_array([[1,0,0,0], [0.707,0,0,0.707]]) + array([quaternion(1, 0, 0, 0), quaternion(0.707, 0, 0, 0.707)]) + """ + return as_quat_array(array) + + # MARK: - Properties + + @property + def w(self) -> float: + return self.__q.w + + @property + def x(self) -> float: + return self.__q.x + + @property + def y(self) -> float: + return self.__q.y + + @property + def z(self) -> float: + return self.__q.z + + @property + def q(self) -> np.quaternion: + return self.__q + + # MARK: - Arithmetic Operators + + def __add__(self, other: FloatOrQuaternion) -> Quaternion: + """Add two quaternions or a quaternion and a scalar.""" + if isinstance(other, Quaternion): + return Quaternion(self.__q + other.q) + elif isinstance(other, (float, int)): + return Quaternion(self.__q + other) + elif isinstance(other, np.quaternion): + return Quaternion(self.__q + other) + return NotImplemented + + def __radd__(self, other: FloatOrQuaternion) -> Quaternion: + """Right addition.""" + return self.__add__(other) + + def __sub__(self, other: FloatOrQuaternion) -> Quaternion: + """Subtract two quaternions or a scalar from a quaternion.""" + if isinstance(other, Quaternion): + return Quaternion(self.__q - other.q) + elif isinstance(other, (float, int)): + return Quaternion(self.__q - other) + elif isinstance(other, np.quaternion): + return Quaternion(self.__q - other) + return NotImplemented + + def __rsub__(self, other: FloatOrQuaternion) -> Quaternion: + """Right subtraction.""" + if isinstance(other, (float, int)): + return Quaternion(other - self.__q) + elif isinstance(other, np.quaternion): + return Quaternion(other - self.__q) + return NotImplemented + + def __mul__(self, other: FloatOrQuaternion) -> Quaternion: + """Multiply two quaternions or a quaternion by a scalar.""" + if isinstance(other, Quaternion): + return Quaternion(self.__q * other.q) + elif isinstance(other, (float, int)): + return Quaternion(self.__q * other) + elif isinstance(other, np.quaternion): + return Quaternion(self.__q * other) + return NotImplemented + + def __rmul__(self, other: FloatOrQuaternion) -> Quaternion: + """Right multiplication.""" + return self.__mul__(other) + + def __truediv__(self, other: FloatOrQuaternion) -> Quaternion: + """Divide a quaternion by another quaternion or scalar.""" + if isinstance(other, Quaternion): + return Quaternion(self.__q / other.q) + elif isinstance(other, (float, int)): + return Quaternion(self.__q / other) + elif isinstance(other, np.quaternion): + return Quaternion(self.__q / other) + return NotImplemented + + def __rtruediv__(self, other: FloatOrQuaternion) -> Quaternion: + """Right division.""" + if isinstance(other, (float, int)): + return Quaternion(other / self.__q) + elif isinstance(other, np.quaternion): + return Quaternion(other / self.__q) + return NotImplemented + + def __pow__(self, exponent: float) -> Quaternion: + """Raise quaternion to a power.""" + return Quaternion(self.__q**exponent) + + def __neg__(self) -> Quaternion: + """Negate the quaternion.""" + return Quaternion(-self.__q) + + def __pos__(self) -> Quaternion: + """Unary positive.""" + return Quaternion(+self.__q) + + def __abs__(self) -> float: + """Return the absolute value (norm) of the quaternion.""" + return float(abs(self.__q)) + + # MARK: - Comparison Operators + + def __eq__(self, other: object) -> bool: + """Check equality with another quaternion.""" + if isinstance(other, Quaternion): + return bool(self.__q == other.q) + elif isinstance(other, np.quaternion): + return bool(self.__q == other) + return False + + def __ne__(self, other: object) -> bool: + """Check inequality with another quaternion.""" + return not self.__eq__(other) + + def __repr__(self) -> str: + """Return a detailed string representation.""" + return f"Quaternion(w={self.w}, x={self.x}, y={self.y}, z={self.z})" + + def __str__(self) -> str: + """Return a user-friendly string representation.""" + return f"quaternion({self.w}, {self.x}, {self.y}, {self.z})" + + # MARK: - Unary Operations + + def conjugate(self) -> Quaternion: + """Return the conjugate of this quaternion. + + The conjugate of q = w + xi + yj + zk is q* = w - xi - yj - zk + """ + return Quaternion(self.__q.conjugate()) + + def inverse(self) -> Quaternion: + """Return the inverse of this quaternion. + + The inverse is defined as q^-1 = q* / |q|^2 + """ + return Quaternion(1.0 / self.__q) + + def normalized(self) -> Quaternion: + """Return a normalized (unit) version of this quaternion.""" + return Quaternion(self.__q / abs(self.__q)) + + # MARK: - Utility Properties + + @property + def norm(self) -> float: + """Return the norm (magnitude) of the quaternion.""" + return float(abs(self.__q)) + + @property + def is_unit(self) -> bool: + """Check if this is a unit quaternion (norm approximately 1).""" + return bool(np.isclose(abs(self.__q), 1.0)) + + @property + def angle(self) -> float: + """Return the rotation angle in radians (for unit quaternions).""" + return float(2.0 * np.arccos(np.clip(self.w, -1.0, 1.0))) + + @property + def axis(self) -> FloatArray3: + """Return the rotation axis as a unit vector (for unit quaternions). + + Returns the zero vector for the identity quaternion. + """ + vect = np.array([self.x, self.y, self.z]) + norm = np.linalg.norm(vect) + if norm < 1e-10: + return np.array([0.0, 0.0, 1.0]) # Default axis for identity + return vect / norm + + # MARK: - Array Conversion Methods + + def as_float_array(self) -> FloatArray4: + """View the quaternion as a float array [w, x, y, z]. + + This function is fast because no data is copied; the returned + quantity is just a "view" of the original. + """ + return as_float_array(self.__q) + + def to_components(self) -> tuple[float, float, float, float]: + """Return quaternion components as a tuple (w, x, y, z).""" + return (self.w, self.x, self.y, self.z) + + # MARK: - Rotation Matrix Conversions + + def to_rotation_matrix(self) -> RotationMatrix: + """Convert this quaternion to a 3x3 rotation matrix. + + For any quaternion q, this returns a matrix m such that, for every + vector v, we have: + m @ v.vec == q * v * q.conjugate() + + Returns: + 3x3 rotation matrix as a numpy array + + Raises: + ZeroDivisionError: If this quaternion has zero norm + """ + return as_rotation_matrix(self.__q) + + @classmethod + def from_rotation_matrix(cls: Type[T], matrix: RotationMatrix) -> T: + """Create a quaternion from a 3x3 rotation matrix. + + Args: + matrix: 3x3 rotation matrix. Should be orthogonal, but the + algorithm handles non-orthogonal matrices by finding + the closest orthogonal matrix. + + Returns: + Unit quaternion representing the rotation + + Raises: + ValueError: If matrix is not 3x3 + """ + q = from_rotation_matrix(matrix) + return cls(q) + + # MARK: - Rotation Vector (Axis-Angle) Conversions + + def to_rotation_vector(self) -> FloatArray3: + """Convert this quaternion to axis-angle representation. + + Returns a vector whose direction is the rotation axis and whose + magnitude is the rotation angle in radians. + + Returns: + 3-element array representing axis-angle rotation + """ + return as_rotation_vector(self.__q) + + @classmethod + def from_rotation_vector(cls: Type[T], rotation_vector: FloatArray3) -> T: + """Create a quaternion from axis-angle representation. + + Args: + rotation_vector: 3-element vector where the direction is the + rotation axis and the magnitude is the angle + in radians. + + Returns: + Unit quaternion representing the rotation + """ + q = from_rotation_vector(rotation_vector) + return cls(q) + + # MARK: - Euler Angle Conversions + + def to_euler_angles(self) -> FloatArray3: + """Convert this quaternion to Euler angles (alpha, beta, gamma). + + WARNING: Euler angles are problematic and should be avoided when + possible. Use quaternions or rotation matrices instead. + + Assumes the Euler angles correspond to the quaternion R via: + R = exp(alpha*z/2) * exp(beta*y/2) * exp(gamma*z/2) + + Returns: + Array of (alpha, beta, gamma) in radians + """ + return as_euler_angles(self.__q) + + @classmethod + def from_euler_angles( + cls: Type[T], + alpha: float, + beta: float, + gamma: float, + ) -> T: + """Create a quaternion from Euler angles. + + WARNING: Euler angles are problematic and should be avoided when + possible. Use quaternions or rotation matrices instead. + + Assumes the Euler angles correspond to the quaternion R via: + R = exp(alpha*z/2) * exp(beta*y/2) * exp(gamma*z/2) + + Args: + alpha: First Euler angle (radians) + beta: Second Euler angle (radians) + gamma: Third Euler angle (radians) + + Returns: + Quaternion representing the rotation + """ + q = from_euler_angles(alpha, beta, gamma) + return cls(q) + + # MARK: - Vector Part Operations + + @classmethod + def from_vector_part(cls: Type[T], vector: FloatArray3) -> T: + """Create a quaternion from a 3D vector (pure quaternion). + + This creates a quaternion with w=0 and vector part equal to the input. + + Args: + vector: 3-element array representing the vector part + + Returns: + Pure quaternion (w=0) + """ + q = from_vector_part(vector) + return cls(q) + + def to_vector_part(self) -> FloatArray3: + """Extract the vector (imaginary) part as a 3D array. + + Returns: + 3-element array [x, y, z] + """ + return as_vector_part(self.__q) + + def rotate_vector(self, vector: FloatArray3) -> FloatArray3: + """Rotate a 3D vector by this quaternion. + + This computes: q * v * q.conjugate() + where v is treated as a pure quaternion (w=0). + + Args: + vector: 3-element vector to rotate + + Returns: + Rotated 3-element vector + """ + return rotate_vectors(self.__q, vector) + + # MARK: - Interpolation Methods + + def slerp(self, other: Quaternion, t: float) -> Quaternion: + """Spherical linear interpolation between this and another quaternion. + + Uses the numpy-quaternions implementation for robust and efficient SLERP. + + Args: + other: Target quaternion + t: Interpolation parameter in [0, 1], where 0 returns self + and 1 returns other + + Returns: + Interpolated quaternion + """ + # Use the package's slerp with t1=0.0, t2=1.0 for normalized interpolation + result = quat_slerp(self.__q, other.q, 0.0, 1.0, t) + return Quaternion(result) + + # MARK: - Comparison Utilities + + def isclose( + self, + other: Quaternion, + rtol: float = 1e-9, + atol: float = 1e-11, + ) -> bool: + """Check if this quaternion is close to another within a tolerance. + + This method accounts for the double cover property of quaternions: + q and -q represent the same rotation, so both are checked. + + Args: + other: Quaternion to compare with + rtol: Relative tolerance (default: 1e-9) + atol: Absolute tolerance (default: 1e-11) + + Returns: + True if quaternions are close within tolerance (accounting for double cover) + """ + # Check both q and -q due to double cover + return bool( + quat_isclose(self.__q, other.q, rtol=rtol, atol=atol) + or quat_isclose(self.__q, -other.q, rtol=rtol, atol=atol) + ) + + @staticmethod + def allclose( + q1: Quaternion, + q2: Quaternion, + rtol: float = 1e-9, + atol: float = 1e-11, + ) -> bool: + """Check if two quaternions are close within a tolerance. + + This method accounts for the double cover property of quaternions: + q and -q represent the same rotation, so both are checked. + + Args: + q1: First quaternion + q2: Second quaternion + rtol: Relative tolerance (default: 1e-9) + atol: Absolute tolerance (default: 1e-11) + + Returns: + True if quaternions are close within tolerance (accounting for double cover) + """ + # Check both q2 and -q2 due to double cover + return bool( + quat_allclose(q1.q, q2.q, rtol=rtol, atol=atol) + or quat_allclose(q1.q, -q2.q, rtol=rtol, atol=atol) + ) + + # MARK: - Special Constructors + + @classmethod + def identity(cls: Type[T]) -> T: + """Create an identity quaternion (no rotation). + + Returns: + Quaternion(1, 0, 0, 0) + """ + return cls.from_components(w=1.0, x=0.0, y=0.0, z=0.0) + + @classmethod + def from_axis_angle(cls: Type[T], axis: FloatArray3, angle: float) -> T: + """Create a quaternion from a rotation axis and angle. + + Args: + axis: 3D rotation axis (will be normalized) + angle: Rotation angle in radians + + Returns: + Unit quaternion representing the rotation + """ + axis_normalized = axis / np.linalg.norm(axis) + half_angle = angle / 2.0 + sin_half = np.sin(half_angle) + return cls.from_components( + w=np.cos(half_angle), + x=axis_normalized[0] * sin_half, + y=axis_normalized[1] * sin_half, + z=axis_normalized[2] * sin_half, + ) + + @classmethod + def from_dict(cls, obj: Any) -> "Quaternion": + """Create a quaternion instance from a dictionary representation. + + Args: + obj: Dictionary containing 'w', 'x', 'y', 'z' keys + + Returns: + A concrete Quaternion instance of the calling class type + """ + w = from_float(obj.get("w")) + x = from_float(obj.get("x")) + y = from_float(obj.get("y")) + z = from_float(obj.get("z")) + # Note: This calls the concrete implementation's from_components + # cls ensures the correct subclass type is returned + return cls.from_components(w, x, y, z) diff --git a/src/pyMathTools/spatial/__init__.py b/src/pyMathTools/spatial/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/tests/test_quaternion.py b/tests/test_quaternion.py new file mode 100644 index 0000000..384379d --- /dev/null +++ b/tests/test_quaternion.py @@ -0,0 +1,846 @@ +""" +Comprehensive unit tests for the Quaternion class using unittest. + +Tests cover: +- Construction methods +- Properties and attributes +- Arithmetic operations +- Comparison operations +- Unary operations +- Rotation conversions +- Vector operations +- Interpolation (SLERP) +""" + +import unittest +import numpy as np +from numpy.testing import assert_allclose, assert_array_almost_equal + +from pyMathTools.spatial.Quaternion import Quaternion + + +class TestQuaternionConstruction(unittest.TestCase): + """Test various construction methods for Quaternion.""" + + def test_from_components_identity(self): + """Test creating identity quaternion from components.""" + q = Quaternion.from_components(1, 0, 0, 0) + self.assertEqual(q.w, 1.0) + self.assertEqual(q.x, 0.0) + self.assertEqual(q.y, 0.0) + self.assertEqual(q.z, 0.0) + + def test_from_components_general(self): + """Test creating general quaternion from components.""" + q = Quaternion.from_components(0.7071, 0.0, 0.0, 0.7071) + assert_allclose(q.w, 0.7071, rtol=1e-4) + assert_allclose(q.x, 0.0, atol=1e-10) + assert_allclose(q.y, 0.0, atol=1e-10) + assert_allclose(q.z, 0.7071, rtol=1e-4) + + def test_from_components_default(self): + """Test default values in from_components.""" + q = Quaternion.from_components() + self.assertEqual(q.w, 1.0) + self.assertEqual(q.x, 0.0) + self.assertEqual(q.y, 0.0) + self.assertEqual(q.z, 0.0) + + def test_identity(self): + """Test identity quaternion constructor.""" + q = Quaternion.identity() + self.assertEqual(q.w, 1.0) + self.assertEqual(q.x, 0.0) + self.assertEqual(q.y, 0.0) + self.assertEqual(q.z, 0.0) + self.assertTrue(q.is_unit) + + def test_from_axis_angle_z_90(self): + """Test construction from axis-angle (90° around Z).""" + axis = np.array([0, 0, 1]) + angle = np.pi / 2 + q = Quaternion.from_axis_angle(axis, angle) + + assert_allclose(q.w, np.cos(np.pi / 4), rtol=1e-10) + assert_allclose(q.x, 0.0, atol=1e-10) + assert_allclose(q.y, 0.0, atol=1e-10) + assert_allclose(q.z, np.sin(np.pi / 4), rtol=1e-10) + + def test_from_axis_angle_x_180(self): + """Test construction from axis-angle (180° around X).""" + axis = np.array([1, 0, 0]) + angle = np.pi + q = Quaternion.from_axis_angle(axis, angle) + + assert_allclose(q.w, 0.0, atol=1e-10) + assert_allclose(q.x, 1.0, rtol=1e-10) + assert_allclose(q.y, 0.0, atol=1e-10) + assert_allclose(q.z, 0.0, atol=1e-10) + + def test_from_axis_angle_unnormalized(self): + """Test that axis gets normalized automatically.""" + axis = np.array([2, 0, 0]) # Not normalized + angle = np.pi / 2 + q = Quaternion.from_axis_angle(axis, angle) + + # Should be same as normalized axis + self.assertTrue(q.is_unit) + + def test_from_vector_part(self): + """Test creating pure quaternion from vector.""" + vec = np.array([1.0, 2.0, 3.0]) + q = Quaternion.from_vector_part(vec) + + self.assertEqual(q.w, 0.0) + self.assertEqual(q.x, 1.0) + self.assertEqual(q.y, 2.0) + self.assertEqual(q.z, 3.0) + + +class TestQuaternionProperties(unittest.TestCase): + """Test quaternion properties and attributes.""" + + def test_norm_identity(self): + """Test norm of identity quaternion.""" + q = Quaternion.identity() + assert_allclose(q.norm, 1.0) + + def test_norm_general(self): + """Test norm of general quaternion.""" + q = Quaternion.from_components(1, 2, 3, 4) + expected_norm = np.sqrt(1 + 4 + 9 + 16) + assert_allclose(q.norm, expected_norm) + + def test_is_unit_true(self): + """Test is_unit for unit quaternion.""" + q = Quaternion.from_axis_angle(np.array([0, 0, 1]), np.pi / 2) + self.assertTrue(q.is_unit) + + def test_is_unit_false(self): + """Test is_unit for non-unit quaternion.""" + q = Quaternion.from_components(1, 1, 1, 1) + self.assertFalse(q.is_unit) + + def test_angle_identity(self): + """Test angle of identity quaternion.""" + q = Quaternion.identity() + assert_allclose(q.angle, 0.0, atol=1e-10) + + def test_angle_90_degrees(self): + """Test angle extraction for 90° rotation.""" + q = Quaternion.from_axis_angle(np.array([1, 0, 0]), np.pi / 2) + assert_allclose(q.angle, np.pi / 2, rtol=1e-10) + + def test_angle_180_degrees(self): + """Test angle extraction for 180° rotation.""" + q = Quaternion.from_axis_angle(np.array([0, 1, 0]), np.pi) + assert_allclose(q.angle, np.pi, rtol=1e-10) + + def test_axis_z(self): + """Test axis extraction for Z-axis rotation.""" + axis = np.array([0, 0, 1]) + q = Quaternion.from_axis_angle(axis, np.pi / 2) + assert_array_almost_equal(q.axis, axis, decimal=10) + + def test_axis_general(self): + """Test axis extraction for general rotation.""" + axis = np.array([1, 1, 1]) / np.sqrt(3) + q = Quaternion.from_axis_angle(axis, np.pi / 4) + assert_array_almost_equal(q.axis, axis, decimal=10) + + def test_axis_identity(self): + """Test axis for identity quaternion (defaults to Z).""" + q = Quaternion.identity() + # Identity should return default axis + self.assertEqual(q.axis.shape, (3,)) + + def test_to_components(self): + """Test components tuple extraction.""" + q = Quaternion.from_components(1, 2, 3, 4) + w, x, y, z = q.to_components() + self.assertEqual(w, 1.0) + self.assertEqual(x, 2.0) + self.assertEqual(y, 3.0) + self.assertEqual(z, 4.0) + + +class TestQuaternionArithmetic(unittest.TestCase): + """Test arithmetic operations on quaternions.""" + + def test_add_quaternions(self): + """Test quaternion addition.""" + q1 = Quaternion.from_components(1, 2, 3, 4) + q2 = Quaternion.from_components(5, 6, 7, 8) + q3 = q1 + q2 + + self.assertEqual(q3.w, 6.0) + self.assertEqual(q3.x, 8.0) + self.assertEqual(q3.y, 10.0) + self.assertEqual(q3.z, 12.0) + + def test_add_scalar(self): + """Test adding scalar to quaternion.""" + q1 = Quaternion.from_components(1, 2, 3, 4) + q2 = q1 + 5 + + self.assertEqual(q2.w, 6.0) + self.assertEqual(q2.x, 2.0) + self.assertEqual(q2.y, 3.0) + self.assertEqual(q2.z, 4.0) + + def test_radd_scalar(self): + """Test right addition with scalar.""" + q1 = Quaternion.from_components(1, 2, 3, 4) + q2 = 5 + q1 + + self.assertEqual(q2.w, 6.0) + self.assertEqual(q2.x, 2.0) + self.assertEqual(q2.y, 3.0) + self.assertEqual(q2.z, 4.0) + + def test_subtract_quaternions(self): + """Test quaternion subtraction.""" + q1 = Quaternion.from_components(5, 6, 7, 8) + q2 = Quaternion.from_components(1, 2, 3, 4) + q3 = q1 - q2 + + self.assertEqual(q3.w, 4.0) + self.assertEqual(q3.x, 4.0) + self.assertEqual(q3.y, 4.0) + self.assertEqual(q3.z, 4.0) + + def test_multiply_quaternions(self): + """Test quaternion multiplication (Hamilton product).""" + # 90° rotation around Z followed by 90° rotation around X + q1 = Quaternion.from_axis_angle(np.array([0, 0, 1]), np.pi / 2) + q2 = Quaternion.from_axis_angle(np.array([1, 0, 0]), np.pi / 2) + q3 = q1 * q2 + + # Result should be a valid rotation + self.assertTrue(q3.is_unit) + + def test_multiply_scalar(self): + """Test multiplying quaternion by scalar.""" + q1 = Quaternion.from_components(1, 2, 3, 4) + q2 = q1 * 2 + + self.assertEqual(q2.w, 2.0) + self.assertEqual(q2.x, 4.0) + self.assertEqual(q2.y, 6.0) + self.assertEqual(q2.z, 8.0) + + def test_divide_quaternions(self): + """Test quaternion division.""" + q1 = Quaternion.from_axis_angle(np.array([0, 0, 1]), np.pi / 2) + q2 = Quaternion.from_axis_angle(np.array([0, 0, 1]), np.pi / 4) + q3 = q1 / q2 + + # Should give approximately pi/4 rotation around Z + self.assertTrue(q3.is_unit) + assert_allclose(q3.angle, np.pi / 4, rtol=1e-10) + + def test_divide_scalar(self): + """Test dividing quaternion by scalar.""" + q1 = Quaternion.from_components(2, 4, 6, 8) + q2 = q1 / 2 + + self.assertEqual(q2.w, 1.0) + self.assertEqual(q2.x, 2.0) + self.assertEqual(q2.y, 3.0) + self.assertEqual(q2.z, 4.0) + + def test_power(self): + """Test quaternion power operation.""" + q1 = Quaternion.from_axis_angle(np.array([0, 0, 1]), np.pi / 2) + q2 = q1 ** 2 + + # (90° rotation)^2 should be 180° rotation + assert_allclose(q2.angle, np.pi, rtol=1e-10) + + def test_negate(self): + """Test quaternion negation.""" + q1 = Quaternion.from_components(1, 2, 3, 4) + q2 = -q1 + + self.assertEqual(q2.w, -1.0) + self.assertEqual(q2.x, -2.0) + self.assertEqual(q2.y, -3.0) + self.assertEqual(q2.z, -4.0) + + def test_positive(self): + """Test unary positive.""" + q1 = Quaternion.from_components(1, 2, 3, 4) + q2 = +q1 + + self.assertEqual(q2.w, 1.0) + self.assertEqual(q2.x, 2.0) + self.assertEqual(q2.y, 3.0) + self.assertEqual(q2.z, 4.0) + + def test_abs(self): + """Test absolute value (norm).""" + q = Quaternion.from_components(1, 2, 3, 4) + expected_norm = np.sqrt(30) + assert_allclose(abs(q), expected_norm) + + +class TestQuaternionComparison(unittest.TestCase): + """Test comparison operations.""" + + def test_equality_same(self): + """Test equality for identical quaternions.""" + q1 = Quaternion.from_components(1, 2, 3, 4) + q2 = Quaternion.from_components(1, 2, 3, 4) + self.assertEqual(q1, q2) + + def test_equality_different(self): + """Test equality for different quaternions.""" + q1 = Quaternion.from_components(1, 2, 3, 4) + q2 = Quaternion.from_components(1, 2, 3, 5) + self.assertNotEqual(q1, q2) + + def test_inequality(self): + """Test inequality operator.""" + q1 = Quaternion.from_components(1, 2, 3, 4) + q2 = Quaternion.from_components(1, 2, 3, 5) + self.assertNotEqual(q1, q2) + + def test_isclose_true(self): + """Test isclose for nearly equal quaternions.""" + q1 = Quaternion.from_components(1.0, 2.0, 3.0, 4.0) + q2 = Quaternion.from_components(1.0 + 1e-10, 2.0, 3.0, 4.0) + self.assertTrue(q1.isclose(q2)) + + def test_isclose_false(self): + """Test isclose for different quaternions.""" + q1 = Quaternion.from_components(1.0, 2.0, 3.0, 4.0) + q2 = Quaternion.from_components(1.1, 2.0, 3.0, 4.0) + self.assertFalse(q1.isclose(q2)) + + def test_allclose_true(self): + """Test static allclose method.""" + q1 = Quaternion.from_components(1.0, 2.0, 3.0, 4.0) + q2 = Quaternion.from_components(1.0 + 1e-10, 2.0, 3.0, 4.0) + self.assertTrue(Quaternion.allclose(q1, q2)) + + def test_repr(self): + """Test __repr__ output.""" + q = Quaternion.from_components(1, 2, 3, 4) + repr_str = repr(q) + self.assertIn("Quaternion", repr_str) + self.assertIn("w=1", repr_str) + self.assertIn("x=2", repr_str) + + def test_str(self): + """Test __str__ output.""" + q = Quaternion.from_components(1, 2, 3, 4) + str_repr = str(q) + self.assertIn("quaternion", str_repr) + + +class TestQuaternionUnaryOperations(unittest.TestCase): + """Test unary operations on quaternions.""" + + def test_conjugate(self): + """Test quaternion conjugate.""" + q = Quaternion.from_components(1, 2, 3, 4) + q_conj = q.conjugate() + + self.assertEqual(q_conj.w, 1.0) + self.assertEqual(q_conj.x, -2.0) + self.assertEqual(q_conj.y, -3.0) + self.assertEqual(q_conj.z, -4.0) + + def test_conjugate_identity(self): + """Test conjugate of identity.""" + q = Quaternion.identity() + q_conj = q.conjugate() + self.assertEqual(q, q_conj) + + def test_inverse_unit(self): + """Test inverse of unit quaternion.""" + q = Quaternion.from_axis_angle(np.array([0, 0, 1]), np.pi / 2) + q_inv = q.inverse() + + # For unit quaternions, inverse equals conjugate + self.assertTrue(q_inv.isclose(q.conjugate())) + + def test_inverse_times_quaternion(self): + """Test that q * q^-1 = identity.""" + q = Quaternion.from_components(1, 2, 3, 4) + q_inv = q.inverse() + result = q * q_inv + + # Should be close to identity + identity = Quaternion.identity() + self.assertTrue(result.isclose(identity, atol=1e-10)) + + def test_normalized_general(self): + """Test normalization of general quaternion.""" + q = Quaternion.from_components(1, 2, 3, 4) + q_norm = q.normalized() + + self.assertTrue(q_norm.is_unit) + assert_allclose(q_norm.norm, 1.0) + + def test_normalized_already_unit(self): + """Test normalizing already unit quaternion.""" + q = Quaternion.identity() + q_norm = q.normalized() + + self.assertTrue(q.isclose(q_norm)) + + +class TestRotationMatrixConversions(unittest.TestCase): + """Test conversions to/from rotation matrices.""" + + def test_to_rotation_matrix_identity(self): + """Test identity quaternion to rotation matrix.""" + q = Quaternion.identity() + matrix = q.to_rotation_matrix() + + expected = np.eye(3) + assert_array_almost_equal(matrix, expected, decimal=10) + + def test_to_rotation_matrix_90z(self): + """Test 90° Z rotation to matrix.""" + q = Quaternion.from_axis_angle(np.array([0, 0, 1]), np.pi / 2) + matrix = q.to_rotation_matrix() + + # 90° Z rotation matrix + expected = np.array([ + [0, -1, 0], + [1, 0, 0], + [0, 0, 1] + ]) + assert_array_almost_equal(matrix, expected, decimal=10) + + def test_to_rotation_matrix_180x(self): + """Test 180° X rotation to matrix.""" + q = Quaternion.from_axis_angle(np.array([1, 0, 0]), np.pi) + matrix = q.to_rotation_matrix() + + # 180° X rotation matrix + expected = np.array([ + [1, 0, 0], + [0, -1, 0], + [0, 0, -1] + ]) + assert_array_almost_equal(matrix, expected, decimal=10) + + def test_from_rotation_matrix_identity(self): + """Test creating quaternion from identity matrix.""" + matrix = np.eye(3) + q = Quaternion.from_rotation_matrix(matrix) + + self.assertTrue(q.isclose(Quaternion.identity())) + + def test_from_rotation_matrix_90z(self): + """Test creating quaternion from 90° Z rotation matrix.""" + matrix = np.array([ + [0, -1, 0], + [1, 0, 0], + [0, 0, 1] + ]) + q = Quaternion.from_rotation_matrix(matrix) + + expected = Quaternion.from_axis_angle(np.array([0, 0, 1]), np.pi / 2) + self.assertTrue(q.isclose(expected)) + + def test_roundtrip_matrix_conversion(self): + """Test quaternion -> matrix -> quaternion roundtrip.""" + q_original = Quaternion.from_axis_angle(np.array([1, 1, 1]) / np.sqrt(3), np.pi / 3) + matrix = q_original.to_rotation_matrix() + q_recovered = Quaternion.from_rotation_matrix(matrix) + + # Should be close (up to sign ambiguity) + self.assertTrue(q_recovered.isclose(q_original) or q_recovered.isclose(-q_original)) + + +class TestRotationVectorConversions(unittest.TestCase): + """Test conversions to/from rotation vectors (axis-angle).""" + + def test_to_rotation_vector_identity(self): + """Test identity quaternion to rotation vector.""" + q = Quaternion.identity() + vec = q.to_rotation_vector() + + assert_array_almost_equal(vec, np.zeros(3), decimal=10) + + def test_to_rotation_vector_90z(self): + """Test 90° Z rotation to rotation vector.""" + q = Quaternion.from_axis_angle(np.array([0, 0, 1]), np.pi / 2) + vec = q.to_rotation_vector() + + expected = np.array([0, 0, np.pi / 2]) + assert_array_almost_equal(vec, expected, decimal=10) + + def test_from_rotation_vector_identity(self): + """Test creating quaternion from zero rotation vector.""" + vec = np.zeros(3) + q = Quaternion.from_rotation_vector(vec) + + self.assertTrue(q.isclose(Quaternion.identity())) + + def test_from_rotation_vector_general(self): + """Test creating quaternion from rotation vector.""" + vec = np.array([0, np.pi / 2, 0]) # 90° around Y + q = Quaternion.from_rotation_vector(vec) + + expected = Quaternion.from_axis_angle(np.array([0, 1, 0]), np.pi / 2) + self.assertTrue(q.isclose(expected)) + + def test_roundtrip_rotation_vector(self): + """Test quaternion -> rotation vector -> quaternion roundtrip.""" + q_original = Quaternion.from_axis_angle(np.array([1, 0, 0]), np.pi / 4) + vec = q_original.to_rotation_vector() + q_recovered = Quaternion.from_rotation_vector(vec) + + self.assertTrue(q_recovered.isclose(q_original)) + + +class TestEulerAngleConversions(unittest.TestCase): + """Test conversions to/from Euler angles.""" + + def test_to_euler_angles_identity(self): + """Test identity quaternion to Euler angles.""" + q = Quaternion.identity() + euler = q.to_euler_angles() + + # Identity should give zero or small angles + self.assertEqual(euler.shape, (3,)) + + def test_from_euler_angles_zeros(self): + """Test creating quaternion from zero Euler angles.""" + q = Quaternion.from_euler_angles(0, 0, 0) + + self.assertTrue(q.isclose(Quaternion.identity())) + + def test_roundtrip_euler_angles(self): + """Test quaternion -> Euler -> quaternion roundtrip.""" + q_original = Quaternion.from_axis_angle(np.array([0, 0, 1]), np.pi / 2) + euler = q_original.to_euler_angles() + q_recovered = Quaternion.from_euler_angles(euler[0], euler[1], euler[2]) + + # Should be close (Euler angles have ambiguities) + self.assertTrue(q_recovered.isclose(q_original) or q_recovered.isclose(-q_original)) + + +class TestVectorOperations(unittest.TestCase): + """Test vector-related operations.""" + + def test_to_vector_part(self): + """Test extracting vector part.""" + q = Quaternion.from_components(1, 2, 3, 4) + vec = q.to_vector_part() + + assert_array_almost_equal(vec, np.array([2, 3, 4])) + + def test_rotate_vector_identity(self): + """Test rotating vector with identity quaternion.""" + q = Quaternion.identity() + vec = np.array([1, 2, 3]) + rotated = q.rotate_vector(vec) + + assert_array_almost_equal(rotated, vec, decimal=10) + + def test_rotate_vector_90z(self): + """Test rotating vector 90° around Z.""" + q = Quaternion.from_axis_angle(np.array([0, 0, 1]), np.pi / 2) + vec = np.array([1, 0, 0]) + rotated = q.rotate_vector(vec) + + expected = np.array([0, 1, 0]) + assert_array_almost_equal(rotated, expected, decimal=10) + + def test_rotate_vector_180x(self): + """Test rotating vector 180° around X.""" + q = Quaternion.from_axis_angle(np.array([1, 0, 0]), np.pi) + vec = np.array([0, 1, 0]) + rotated = q.rotate_vector(vec) + + expected = np.array([0, -1, 0]) + assert_array_almost_equal(rotated, expected, decimal=10) + + def test_rotate_vector_general(self): + """Test general vector rotation.""" + q = Quaternion.from_axis_angle(np.array([1, 1, 1]) / np.sqrt(3), 2 * np.pi / 3) + vec = np.array([1, 0, 0]) + rotated = q.rotate_vector(vec) + + # 120° rotation around (1,1,1) should permute coordinates + expected = np.array([0, 1, 0]) + assert_array_almost_equal(rotated, expected, decimal=8) + + +class TestSLERP(unittest.TestCase): + """Test spherical linear interpolation.""" + + def test_slerp_endpoints(self): + """Test SLERP at endpoints.""" + q1 = Quaternion.identity() + q2 = Quaternion.from_axis_angle(np.array([0, 0, 1]), np.pi / 2) + + # At t=0, should get q1 + q_start = q1.slerp(q2, 0.0) + self.assertTrue(q_start.isclose(q1)) + + # At t=1, should get q2 + q_end = q1.slerp(q2, 1.0) + self.assertTrue(q_end.isclose(q2)) + + def test_slerp_midpoint(self): + """Test SLERP at midpoint.""" + q1 = Quaternion.identity() + q2 = Quaternion.from_axis_angle(np.array([0, 0, 1]), np.pi / 2) + + q_mid = q1.slerp(q2, 0.5) + + # At t=0.5, should be 45° rotation + expected = Quaternion.from_axis_angle(np.array([0, 0, 1]), np.pi / 4) + self.assertTrue(q_mid.isclose(expected)) + + def test_slerp_unit_preservation(self): + """Test that SLERP preserves unit norm.""" + q1 = Quaternion.from_axis_angle(np.array([1, 0, 0]), np.pi / 4) + q2 = Quaternion.from_axis_angle(np.array([0, 1, 0]), np.pi / 3) + + for t in np.linspace(0, 1, 10): + q_interp = q1.slerp(q2, t) + self.assertTrue(q_interp.is_unit) + + def test_slerp_smooth_interpolation(self): + """Test that SLERP produces smooth interpolation.""" + q1 = Quaternion.identity() + q2 = Quaternion.from_axis_angle(np.array([0, 0, 1]), np.pi) + + angles = [] + for t in np.linspace(0, 1, 11): + q_interp = q1.slerp(q2, t) + angles.append(q_interp.angle) + + # Angles should increase monotonically + for i in range(len(angles) - 1): + self.assertLessEqual(angles[i], angles[i + 1]) + + def test_slerp_opposite_quaternions(self): + """Test SLERP between opposite quaternions.""" + q1 = Quaternion.from_axis_angle(np.array([0, 0, 1]), 0) + q2 = Quaternion.from_axis_angle(np.array([0, 0, 1]), np.pi) + + q_mid = q1.slerp(q2, 0.5) + + # Should interpolate correctly + assert_allclose(q_mid.angle, np.pi / 2, rtol=1e-10) + + +class TestArrayConversions(unittest.TestCase): + """Test array conversion methods.""" + + def test_as_float_array(self): + """Test conversion to float array.""" + q = Quaternion.from_components(1, 2, 3, 4) + arr = q.as_float_array() + + expected = np.array([1, 2, 3, 4]) + assert_array_almost_equal(arr, expected) + + def test_as_float_array_shape(self): + """Test shape of float array.""" + q = Quaternion.from_components(1, 2, 3, 4) + arr = q.as_float_array() + + self.assertEqual(arr.shape, (4,)) + + +class TestEdgeCases(unittest.TestCase): + """Test edge cases and special scenarios.""" + + def test_very_small_rotation(self): + """Test very small rotation angle.""" + q = Quaternion.from_axis_angle(np.array([1, 0, 0]), 1e-10) + + self.assertTrue(q.isclose(Quaternion.identity())) + + def test_very_large_angle(self): + """Test angle larger than 2π.""" + # Should wrap around + q = Quaternion.from_axis_angle(np.array([0, 0, 1]), 3 * np.pi) + + # 3π should be equivalent to π rotation + expected = Quaternion.from_axis_angle(np.array([0, 0, 1]), np.pi) + self.assertTrue(q.isclose(expected) or q.isclose(-expected)) + + def test_zero_norm_handling(self): + """Test that operations handle near-zero quaternions appropriately.""" + # This would be an invalid quaternion, but we test defensive handling + # Note: actual zero quaternion would cause division by zero in some operations + q = Quaternion.from_components(1e-15, 1e-15, 1e-15, 1e-15) + + # Should still be able to normalize + q_norm = q.normalized() + self.assertTrue(q_norm.is_unit) + + +class TestQuaternionSerialization(unittest.TestCase): + """Test serialization/deserialization methods and DataModelHelper inheritance.""" + + def test_to_dict_identity(self): + """Test to_dict for identity quaternion.""" + q = Quaternion.identity() + d = q.to_dict() + + self.assertIsInstance(d, dict) + self.assertIn("w", d) + self.assertIn("x", d) + self.assertIn("y", d) + self.assertIn("z", d) + self.assertEqual(d["w"], 1.0) + self.assertEqual(d["x"], 0.0) + self.assertEqual(d["y"], 0.0) + self.assertEqual(d["z"], 0.0) + + def test_to_dict_general(self): + """Test to_dict for general quaternion.""" + q = Quaternion.from_components(1.0, 2.0, 3.0, 4.0) + d = q.to_dict() + + self.assertEqual(d["w"], 1.0) + self.assertEqual(d["x"], 2.0) + self.assertEqual(d["y"], 3.0) + self.assertEqual(d["z"], 4.0) + + def test_to_dict_normalized(self): + """Test to_dict preserves values for normalized quaternion.""" + q = Quaternion.from_axis_angle(np.array([1, 0, 0]), np.pi / 2) + d = q.to_dict() + + # Check that all components are present and are floats + self.assertIsInstance(d["w"], float) + self.assertIsInstance(d["x"], float) + self.assertIsInstance(d["y"], float) + self.assertIsInstance(d["z"], float) + + # Verify components match + self.assertAlmostEqual(d["w"], q.w) + self.assertAlmostEqual(d["x"], q.x) + self.assertAlmostEqual(d["y"], q.y) + self.assertAlmostEqual(d["z"], q.z) + + def test_from_dict_identity(self): + """Test from_dict for identity quaternion.""" + d = {"w": 1.0, "x": 0.0, "y": 0.0, "z": 0.0} + q = Quaternion.from_dict(d) + + self.assertIsInstance(q, Quaternion) + self.assertEqual(q.w, 1.0) + self.assertEqual(q.x, 0.0) + self.assertEqual(q.y, 0.0) + self.assertEqual(q.z, 0.0) + + def test_from_dict_general(self): + """Test from_dict for general quaternion.""" + d = {"w": 1.0, "x": 2.0, "y": 3.0, "z": 4.0} + q = Quaternion.from_dict(d) + + self.assertEqual(q.w, 1.0) + self.assertEqual(q.x, 2.0) + self.assertEqual(q.y, 3.0) + self.assertEqual(q.z, 4.0) + + def test_from_dict_with_integers(self): + """Test from_dict accepts integer values.""" + d = {"w": 1, "x": 2, "y": 3, "z": 4} + q = Quaternion.from_dict(d) + + self.assertEqual(q.w, 1.0) + self.assertEqual(q.x, 2.0) + self.assertEqual(q.y, 3.0) + self.assertEqual(q.z, 4.0) + + def test_roundtrip_serialization_identity(self): + """Test roundtrip: quaternion -> dict -> quaternion for identity.""" + q_original = Quaternion.identity() + d = q_original.to_dict() + q_recovered = Quaternion.from_dict(d) + + self.assertTrue(q_original.isclose(q_recovered)) + + def test_roundtrip_serialization_general(self): + """Test roundtrip: quaternion -> dict -> quaternion for general quaternion.""" + q_original = Quaternion.from_components(1.0, 2.0, 3.0, 4.0) + d = q_original.to_dict() + q_recovered = Quaternion.from_dict(d) + + self.assertTrue(q_original.isclose(q_recovered)) + + def test_roundtrip_serialization_rotation(self): + """Test roundtrip for rotation quaternion.""" + q_original = Quaternion.from_axis_angle(np.array([1, 1, 1]) / np.sqrt(3), np.pi / 3) + d = q_original.to_dict() + q_recovered = Quaternion.from_dict(d) + + self.assertTrue(q_original.isclose(q_recovered)) + + def test_roundtrip_with_negatives(self): + """Test roundtrip with negative components.""" + q_original = Quaternion.from_components(-0.5, -0.5, 0.5, 0.5) + d = q_original.to_dict() + q_recovered = Quaternion.from_dict(d) + + self.assertEqual(q_recovered.w, -0.5) + self.assertEqual(q_recovered.x, -0.5) + self.assertEqual(q_recovered.y, 0.5) + self.assertEqual(q_recovered.z, 0.5) + + def test_to_dict_returns_new_dict(self): + """Test that to_dict returns a new dictionary each time.""" + q = Quaternion.from_components(1.0, 2.0, 3.0, 4.0) + d1 = q.to_dict() + d2 = q.to_dict() + + self.assertIsNot(d1, d2) + self.assertEqual(d1, d2) + + def test_dict_modification_doesnt_affect_quaternion(self): + """Test that modifying the dict doesn't affect the quaternion.""" + q = Quaternion.from_components(1.0, 2.0, 3.0, 4.0) + d = q.to_dict() + + # Modify the dictionary + d["w"] = 999.0 + d["x"] = 999.0 + + # Original quaternion should be unchanged + self.assertEqual(q.w, 1.0) + self.assertEqual(q.x, 2.0) + + def test_inheritance_from_quaternion_type(self): + """Test that Quaternion properly inherits from QuaternionType.""" + from foundationTypes.mathTypes.QuaternionType import QuaternionType + + q = Quaternion.from_components(1.0, 2.0, 3.0, 4.0) + self.assertIsInstance(q, QuaternionType) + + def test_inheritance_from_data_model_helper(self): + """Test that Quaternion inherits from DataModelHelper.""" + from foundationTypes.dataModelHelper import DataModelHelper + + q = Quaternion.from_components(1.0, 2.0, 3.0, 4.0) + self.assertIsInstance(q, DataModelHelper) + + def test_has_serialization_methods(self): + """Test that Quaternion has the required serialization methods.""" + q = Quaternion.identity() + + # Check methods exist + self.assertTrue(hasattr(q, "to_dict")) + self.assertTrue(hasattr(Quaternion, "from_dict")) + + # Check they're callable + self.assertTrue(callable(q.to_dict)) + self.assertTrue(callable(Quaternion.from_dict)) + + +if __name__ == "__main__": + unittest.main() diff --git a/tests/testpyMathTools.py b/tests/testpyMathTools.py deleted file mode 100644 index 5bb39ca..0000000 --- a/tests/testpyMathTools.py +++ /dev/null @@ -1,11 +0,0 @@ -import pytest - -# from pyMathTools import pyMathTools - - -def test_content(): - """Sample pytest test function.""" - - -if __name__ == "__main__": - test_content() From 3a5bffc4f9504000e74c598df654b0982d6a3110 Mon Sep 17 00:00:00 2001 From: kopecn Date: Tue, 2 Dec 2025 14:34:32 -0800 Subject: [PATCH 15/22] cleanup and adding quant visualizers. --- examples/plotQuatUnitCircles.py | 2 +- .../generators/sphericalGenerators.py | 44 ++- ...alSmallCircles.py => plotUnitSpherical.py} | 254 +++++++++--------- 3 files changed, 174 insertions(+), 126 deletions(-) rename src/pyMathToolsPlotHelpers/{plotUnitSphericalSmallCircles.py => plotUnitSpherical.py} (81%) diff --git a/examples/plotQuatUnitCircles.py b/examples/plotQuatUnitCircles.py index 5ac5ea2..d59374a 100644 --- a/examples/plotQuatUnitCircles.py +++ b/examples/plotQuatUnitCircles.py @@ -32,7 +32,7 @@ quat_to_unitSphericalSmallCircle, ) -from pyMathToolsPlotHelpers.plotUnitSphericalSmallCircles import ( +from pyMathToolsPlotHelpers.plotUnitSpherical import ( plot_unit_spherical_multiplot, ) diff --git a/src/pyMathTools/generators/sphericalGenerators.py b/src/pyMathTools/generators/sphericalGenerators.py index b01be06..07d53bb 100644 --- a/src/pyMathTools/generators/sphericalGenerators.py +++ b/src/pyMathTools/generators/sphericalGenerators.py @@ -17,7 +17,7 @@ from foundationTypes.mathTypes.UnitSphericalSmallCircle import UnitSphericalSmallCircle from foundationTypes.mathTypes.UnitSphericalArc import UnitSphericalArc -from pyMathTools.hints import FloatNDArray +from pyMathTools.hints import FloatNDArray, FloatArray3 def generate_spherical_small_circle_points( @@ -224,3 +224,45 @@ def generate_spherical_arc_points( azimuth_points[i] = arctan2(y, x) # azimuth return azimuth_points, polar_points + + +def quaternion_to_spherical_vector(quaternion) -> Tuple[FloatArray3, float, float]: + """ + Convert a quaternion to a unit vector in spherical coordinates. + + The quaternion is applied to a reference direction (positive Z-axis) + to get the pointing direction, which is then converted to spherical coordinates. + + Parameters: + ----------- + quaternion : Quaternion + A quaternion object representing a rotation + + Returns: + -------- + cartesian_point : FloatArray3 + The Cartesian coordinates [x, y, z] of the rotated unit vector + azimuth : float + Azimuth angle in radians (-π to π) + polar : float + Polar angle in radians (0 to π, physics convention: 0 = north pole) + """ + # Reference direction: positive Z-axis (north pole) + reference_direction = array([0.0, 0.0, 1.0]) + + # Apply quaternion rotation to the reference direction + rotated_vector = quaternion.rotate_vector(reference_direction) + + # Normalize to ensure it's a unit vector + x, y, z = rotated_vector + r = sqrt(x**2 + y**2 + z**2) + if r > 1e-10: + x, y, z = x / r, y / r, z / r + + # Convert to spherical coordinates (physics convention) + polar = arccos(clip(z, -1, 1)) # colatitude (0 to π) + azimuth = arctan2(y, x) # azimuth (-π to π) + + cartesian_point = array([x, y, z]) + + return cartesian_point, azimuth, polar diff --git a/src/pyMathToolsPlotHelpers/plotUnitSphericalSmallCircles.py b/src/pyMathToolsPlotHelpers/plotUnitSpherical.py similarity index 81% rename from src/pyMathToolsPlotHelpers/plotUnitSphericalSmallCircles.py rename to src/pyMathToolsPlotHelpers/plotUnitSpherical.py index 8dedb97..cbbaaa1 100644 --- a/src/pyMathToolsPlotHelpers/plotUnitSphericalSmallCircles.py +++ b/src/pyMathToolsPlotHelpers/plotUnitSpherical.py @@ -7,10 +7,12 @@ from foundationTypes.mathTypes.UnitSphericalSmallCircle import UnitSphericalSmallCircle from foundationTypes.mathTypes.UnitSphericalArc import UnitSphericalArc +from pyMathTools.spatial.Quaternion import Quaternion from pyMathTools.transforms.sphericalTransforms import plate_carree_transform from pyMathTools.generators.sphericalGenerators import ( generate_spherical_small_circle_points, generate_spherical_arc_points, + quaternion_to_spherical_vector, ) deg45 = np.deg2rad(45) @@ -19,129 +21,10 @@ deg22_5 = np.pi / 8 -def plot_unit_spherical( - circles: List[UnitSphericalSmallCircle] | None = None, - arcs: List[UnitSphericalArc] | None = None, - num_points: int = 120, - figsize: Tuple[int, int] = (10, 6), - title: str = "Small Circles on Sphere (Plate Carrée Projection)", - show_legend: bool = True, - show_plot: bool = False, -) -> Tuple[Figure, Axes]: - """ - Plot small circles and arcs on a sphere using plate carrée (equirectangular) projection. - - Parameters: - ----------- - circles : List[UnitSphericalSmallCircle], optional - List of UnitSphericalSmallCircle objects. - Each circle represents a small circle on the unit sphere - arcs : List[UnitSphericalArc], optional - List of UnitSphericalArc objects. - Each arc starts at (azimuth, polar) and extends along a great circle - for the specified arc_length (positive or negative) - num_points : int, optional - Number of points to generate for each circle/arc (default: 120) - figsize : tuple, optional - Figure size (default: (10, 6)) - title : str, optional - Plot title - show_legend : bool, optional - Whether to show the legend (default: True) - show_plot : bool, optional - Whether to display the plot (default: False) - - Returns: - -------- - fig, ax : matplotlib figure and axes objects - """ - fig, ax = plt.subplots(figsize=figsize) - - # Generate points for each circle - if circles is not None: - for idx, circle in enumerate(circles): - - # Generate true spherical small circle points - azimuth_points, polar_points = generate_spherical_small_circle_points( - circle, num_points - ) - - # Apply Plate Carrée projection - x_points, y_points = plate_carree_transform(azimuth_points, polar_points) - - # Convert radius to degrees for legend - radius_deg = np.degrees(circle.radius_angle) - azimuth_deg = np.degrees(circle.azimuth) - polar_deg = np.degrees(circle.polar) - label = f"Circle {idx + 1}: A:{azimuth_deg:.1f}°, P:{polar_deg:.1f}°, R:{radius_deg:.1f}°" - - # Plot the circle - ax.scatter( - np.degrees(x_points), - np.degrees(y_points), - s=10, - alpha=0.6, - label=label, - ) - - # Generate points for each arc - if arcs is not None: - for idx, arc in enumerate(arcs): - - # Generate spherical arc points - azimuth_points, polar_points = generate_spherical_arc_points( - arc, num_points - ) - - # Apply Plate Carrée projection - x_points, y_points = plate_carree_transform(azimuth_points, polar_points) - - # Convert to degrees for legend - arc_length_deg = np.degrees(arc.arc_length) - azimuth_deg = np.degrees(arc.azimuth) - polar_deg = np.degrees(arc.polar) - label = f"Arc {idx + 1}: A:{azimuth_deg:.1f}°, P:{polar_deg:.1f}°, L:{arc_length_deg:.1f}°" - - # Plot the arc as a line - ax.plot( - np.degrees(x_points), - np.degrees(y_points), - linewidth=2, - alpha=0.8, - label=label, - marker="o", - markersize=3, - ) - - # Set up axes - ax.set_xlabel("Azimuth (deg)", fontsize=12) - ax.set_ylabel("Latitude (deg)", fontsize=12) - ax.set_xlim(-180, 180) - - # Add grid - ax.grid(True, alpha=0.3) - - # Set x-axis ticks - ax.set_xticks([-180, -90, 0, 90, 180]) - ax.set_xticklabels(["-180°", "-90°", "0°", "90°", "180°"]) - - ax.set_title(title, fontsize=14) - - # Add legend if requested - if show_legend: - ax.legend() - - plt.tight_layout() - - if show_plot: - plt.show() - - return fig, ax - - def plot_unit_spherical_advanced( circles: List[UnitSphericalSmallCircle] | None = None, arcs: List[UnitSphericalArc] | None = None, + quaternions: List[Quaternion] | None = None, num_points: int = 120, figsize: Tuple[int, int] = (10, 6), title: str = "Small Circles on Sphere (Flattened)", @@ -162,6 +45,9 @@ def plot_unit_spherical_advanced( List of UnitSphericalArc objects. Each arc starts at (azimuth, polar) and extends along a great circle for the specified arc_length (positive or negative) + quaternions : List[Quaternion], optional + List of Quaternion objects. + Each quaternion is visualized as a unit vector from origin to its pointing direction num_points : int, optional Number of points to generate around each circle figsize : tuple, optional @@ -258,6 +144,51 @@ def plot_unit_spherical_advanced( markersize=3, ) + # Plot quaternions as vectors + if quaternions is not None: + for idx, quat in enumerate(quaternions): + # Convert quaternion to spherical coordinates + _, azimuth, polar = quaternion_to_spherical_vector(quat) + + # Apply Plate Carrée projection to the endpoint + x_end, y_end = plate_carree_transform(azimuth, polar) + + # Origin is at (0, 0) in this projection + x_origin, y_origin = plate_carree_transform(0, np.pi / 2) + + # Determine color and label + color = colors[idx] if colors is not None and idx < len(colors) else None + label = labels[idx] if labels is not None and idx < len(labels) else None + + if label is None: + azimuth_deg = np.degrees(azimuth) + polar_deg = np.degrees(polar) + label = f"Quat {idx + 1}: w:{quat.w:.1f}, x:{quat.x:.1f}, y:{quat.x:.1f}, z:{quat.x:.1f}" + + # Plot vector as an arrow from origin to the point + ax.annotate( + "", + xy=(np.degrees(x_end), np.degrees(y_end)), + xytext=(np.degrees(x_origin), np.degrees(y_origin)), + arrowprops=dict( + arrowstyle="->", + lw=2, + color=color if color else f"C{idx}", + alpha=0.8, + ), + ) + + # Plot endpoint + ax.scatter( + np.degrees(x_end), + np.degrees(y_end), + s=100, + marker="*", + c=color, + label=label, + zorder=10, + ) + # Set up axes ax.set_xlabel("Azimuth (deg)", fontsize=12) ax.set_ylabel("Latitude (deg)", fontsize=12) @@ -287,6 +218,7 @@ def plot_unit_spherical_advanced( def plot_unit_spherical_polar( circles: List[UnitSphericalSmallCircle] | None = None, arcs: List[UnitSphericalArc] | None = None, + quaternions: List[Quaternion] | None = None, num_points: int = 120, figsize: Tuple[int, int] = (8, 8), title: str = "Small Circles on Sphere (Polar View)", @@ -407,7 +339,8 @@ def plot_unit_spherical_polar( arc_length_deg = np.degrees(arc.arc_length) azimuth_deg = np.degrees(arc.azimuth) polar_deg = np.degrees(arc.polar) - label = f"Arc {idx + 1}: A:{azimuth_deg:.1f}°, P:{polar_deg:.1f}°, L:{arc_length_deg:.1f}°" + orient_deg = np.degrees(arc.orient) + label = f"Arc {idx + 1}: A:{azimuth_deg:.1f}°, P:{polar_deg:.1f}°, L:{arc_length_deg:.1f}°, O:{orient_deg:.1f}°" # Plot the arc ax.plot( @@ -420,6 +353,30 @@ def plot_unit_spherical_polar( markersize=3, ) + # Plot quaternions as vectors + if quaternions is not None: + for idx, quat in enumerate(quaternions): + # Convert quaternion to spherical coordinates + _, azimuth, polar = quaternion_to_spherical_vector(quat) + + # Map polar angle to radius (0 to 1) for polar plot + radius = (polar + np.pi / 2) / np.pi + + # Plot vector as an arrow from origin + ax.annotate( + "", + xy=(azimuth, radius), + xytext=(0, 0), + arrowprops=dict(arrowstyle="->", lw=2, color=f"C{idx}", alpha=0.8), + ) + + # Plot endpoint + azimuth_deg = np.degrees(azimuth) + polar_deg = np.degrees(polar) + label = f"Quat {idx + 1}: w:{quat.w:.1f}, x:{quat.x:.1f}, y:{quat.x:.1f}, z:{quat.x:.1f}" + + ax.scatter(azimuth, radius, s=150, marker="*", label=label, zorder=10) + ax.set_title(title, fontsize=14, pad=20) if labels is not None: @@ -439,6 +396,7 @@ def plot_unit_spherical_polar( def plot_unit_spherical_3d( circles: List[UnitSphericalSmallCircle] | None = None, arcs: List[UnitSphericalArc] | None = None, + quaternions: List[Quaternion] | None = None, num_points: int = 120, figsize: Tuple[int, int] = (10, 10), title: str = "Small Circles on Unit Sphere (3D)", @@ -538,7 +496,8 @@ def plot_unit_spherical_3d( arc_length_deg = np.degrees(arc.arc_length) azimuth_deg = np.degrees(arc.azimuth) polar_deg = np.degrees(arc.polar) - label = f"Arc {idx + 1}: A:{azimuth_deg:.1f}°, P:{polar_deg:.1f}°, L:{arc_length_deg:.1f}°" + orient_deg = np.degrees(arc.orient) + label = f"Arc {idx + 1}: A:{azimuth_deg:.1f}°, P:{polar_deg:.1f}°, L:{arc_length_deg:.1f}°, O:{orient_deg:.1f}°" # Plot the arc as a line ax.plot( @@ -552,6 +511,35 @@ def plot_unit_spherical_3d( markersize=3, ) + # Plot quaternions as vectors + if quaternions is not None: + for idx, quat in enumerate(quaternions): + # Convert quaternion to Cartesian coordinates + cartesian, azimuth, polar = quaternion_to_spherical_vector(quat) + x, y, z = cartesian + + # Plot vector from origin to the point + ax.quiver( + 0, + 0, + 0, + x, + y, + z, + length=1.0, + arrow_length_ratio=0.15, + color=f"C{idx}", + linewidth=2.5, + alpha=0.9, + ) + + # Plot endpoint + azimuth_deg = np.degrees(azimuth) + polar_deg = np.degrees(polar) + label = f"Quat {idx + 1}: w:{quat.w:.1f}, x:{quat.x:.1f}, y:{quat.x:.1f}, z:{quat.x:.1f}" + + ax.scatter(x, y, z, s=150, marker="*", label=label, zorder=10) + ax.set_xlabel("X") ax.set_ylabel("Y") ax.set_zlabel("Z") @@ -580,6 +568,7 @@ def plot_unit_spherical_3d( def plot_unit_spherical_multiplot( circles: List[UnitSphericalSmallCircle] | None = None, arcs: List[UnitSphericalArc] | None = None, + quaternions: List[Quaternion] | None = None, num_points: int = 120, figsize: Tuple[int, int] = (24, 6), title: str = "Small Circles on Unit Sphere - Multiple Views", @@ -618,6 +607,7 @@ def plot_unit_spherical_multiplot( plot_unit_spherical_advanced( circles=circles, arcs=arcs, + quaternions=quaternions, num_points=num_points, title="Plate Carrée Projection", show_centers=False, @@ -632,6 +622,7 @@ def plot_unit_spherical_multiplot( plot_unit_spherical_polar( circles=circles, arcs=arcs, + quaternions=quaternions, num_points=num_points, title="Top View (Down Z-Axis)", ax=ax2, @@ -645,6 +636,7 @@ def plot_unit_spherical_multiplot( plot_unit_spherical_3d( circles=circles, arcs=arcs, + quaternions=quaternions, num_points=num_points, title="3D View", show_sphere=True, @@ -732,9 +724,23 @@ def demo() -> None: UnitSphericalArc(orient=0, azimuth=deg90, polar=deg45, arc_length=-deg22_5), ] + quats: List[Quaternion] = [ + Quaternion.from_components(w=1, x=0, y=0, z=0), + Quaternion.from_components(w=0, x=1, y=0, z=0), + Quaternion.from_components(w=0, x=0, y=1, z=0), + Quaternion.from_components(w=0, x=0, y=0, z=1), + Quaternion.from_components(w=0, x=0.707, y=0, z=0.707), + Quaternion.from_components(w=0, x=0.707, y=0.707, z=0), + Quaternion.from_components(w=0, x=0, y=0.707, z=0.707), + Quaternion.from_components(w=0.707, x=0.707, y=0, z=0), + Quaternion.from_components(w=0.707, x=0, y=0.707, z=0), + Quaternion.from_components(w=0.707, x=0, y=0, z=0.707), + ] + _ = plot_unit_spherical_multiplot( - circles=circles, - arcs=arcs, + # circles=circles, + # arcs=arcs, + quaternions=quats, show_plot=True, ) From 2c13c753240ef7599dbcc15cbc4176d6dd803d57 Mon Sep 17 00:00:00 2001 From: kopecn Date: Tue, 2 Dec 2025 23:07:49 -0800 Subject: [PATCH 16/22] cont... --- .../constructUnitSphericalSmallCircle.py | 46 ++++- .../generators/sphericalGenerators.py | 176 +++++++++++++++--- .../transforms/sphericalTransforms.py | 45 +++-- .../plotUnitSpherical.py | 7 +- 4 files changed, 231 insertions(+), 43 deletions(-) diff --git a/src/pyMathTools/constructors/constructUnitSphericalSmallCircle.py b/src/pyMathTools/constructors/constructUnitSphericalSmallCircle.py index 5cbffab..e657853 100644 --- a/src/pyMathTools/constructors/constructUnitSphericalSmallCircle.py +++ b/src/pyMathTools/constructors/constructUnitSphericalSmallCircle.py @@ -1,3 +1,10 @@ +""" +Constructor functions for UnitSphericalSmallCircle objects. + +This module provides utilities for creating UnitSphericalSmallCircle instances +from various representations, using ISO Physics Convention for spherical coordinates. +""" + import numpy as np from numpy import atan2, acos, pi @@ -9,17 +16,46 @@ def quat_to_unitSphericalSmallCircle( q: Quaternion, radius_angle: float = pi / 4 ) -> UnitSphericalSmallCircle: - # Extract vector part + """ + Convert a quaternion's vector part to a UnitSphericalSmallCircle. + + Uses ISO Physics Convention for spherical coordinates: + - φ (azimuth): angle in xy-plane from +x axis, range [0, 2π) + - θ (polar): angle from +z axis (colatitude), range [0, π] + Parameters: + ----------- + q : Quaternion + Quaternion whose vector part (x, y, z) defines the center direction + radius_angle : float + Angular radius of the small circle (default: π/4) + + Returns: + -------- + UnitSphericalSmallCircle + Small circle centered at the direction defined by the quaternion's vector part + + Notes: + ------ + - If vector part is zero, defaults to north pole (azimuth=0, polar=π/2) + - Uses standard spherical-to-Cartesian formulas: + x = r*sin(θ)*cos(φ) + y = r*sin(θ)*sin(φ) + z = r*cos(θ) + """ + # Extract vector part (this represents a direction in Cartesian space) r = np.sqrt(q.x**2 + q.y**2 + q.z**2) if r == 0: - # Default to north pole if vector part is zero + # Default to equator at +x axis if vector part is zero azimuth = 0.0 polar = pi / 2 else: - polar = atan2(q.y, q.x) - azimuth = acos(q.z / r) + # Convert Cartesian (x,y,z) to spherical (azimuth, polar) + # ISO Physics Convention: + azimuth = atan2(q.y, q.x) # φ = atan2(y, x) - angle in xy-plane from +x + polar = acos(q.z / r) # θ = acos(z/r) - angle from +z axis + return UnitSphericalSmallCircle( - azimuth=polar, polar=azimuth, radius_angle=radius_angle + azimuth=azimuth, polar=polar, radius_angle=radius_angle ) diff --git a/src/pyMathTools/generators/sphericalGenerators.py b/src/pyMathTools/generators/sphericalGenerators.py index 07d53bb..a547f61 100644 --- a/src/pyMathTools/generators/sphericalGenerators.py +++ b/src/pyMathTools/generators/sphericalGenerators.py @@ -1,3 +1,41 @@ +""" +Generators for spherical geometry primitives. + +This module provides functions for generating points on spherical circles, arcs, +and converting quaternions to spherical vectors. + +ISO Physics Convention for Spherical Coordinates +------------------------------------------------- +All functions in this module use the ISO physics convention: + +Notation: (r, θ, φ) where r=1 for unit sphere + +θ (theta) - Polar Angle (Colatitude): + - Measured from +z axis downward + - Range: [0, π] + - θ = 0: North pole (+z axis) + - θ = π/2: Equator (xy-plane) + - θ = π: South pole (-z axis) + +φ (phi) - Azimuthal Angle: + - Measured from +x axis in xy-plane, counterclockwise when viewed from above + - Range: [0, 2π) + - φ = 0: +x axis + - φ = π/2: +y axis + - φ = π: -x axis + - φ = 3π/2: -y axis + +Cartesian Conversion: + x = r*sin(θ)*cos(φ) + y = r*sin(θ)*sin(φ) + z = r*cos(θ) + +Inverse Conversion: + θ = arccos(z/r) + φ = arctan2(y, x) + r = sqrt(x² + y² + z²) +""" + from typing import Tuple from numpy import ( clip, @@ -30,24 +68,24 @@ def generate_spherical_small_circle_points( A small circle is the set of all points at a fixed angular distance (radius_angle) from a center point on the sphere. + Uses ISO Physics Convention (see module docstring for details). + Parameters: ----------- - circle.azimuth : float or NDArray - Azimuth angle of circle center (radians, -2π to 2π) - circle.polar : float or NDArray - Polar angle of circle center (radians, physics convention: 0 to π) - where 0 = north pole, π/2 = equator, π = south pole - circle.radius_angle : float or NDArray - Angular radius of the circle (radians) + circle : UnitSphericalSmallCircle + Small circle specification with: + - azimuth (φ): angle in xy-plane from +x axis [0, 2π) + - polar (θ): angle from +z axis (colatitude) [0, π] + - radius_angle: angular radius of the circle num_points : int Number of points to generate around the circle Returns: -------- azimuth_points : NDArray - Azimuth angles of points on the circle (radians) + Azimuth angles (φ) of points on the circle in radians polar_points : NDArray - Polar angles of points on the circle (physics convention: 0 to π, colatitude) + Polar angles (θ) of points on the circle in radians """ # Use physics convention directly: polar ∈ [0, π], where 0 = north pole, π/2 = equator # Convert center to Cartesian coordinates (physics convention) @@ -131,24 +169,26 @@ def generate_spherical_arc_points( vector (from origin through the start point), and your fingers curl in the positive direction of rotation (positive orient). + Uses ISO Physics Convention (see module docstring for details). + Parameters: ----------- arc : UnitSphericalArc Arc specification with: - - azimuth: starting point longitude (radians) - - polar: starting point colatitude in physics convention (0 to π, 0 = north pole) - - orient: rotation angle about the radial vector (radians, -π to π) - Following right-hand rule: thumb along radial, fingers curl in +orient direction - - arc_length: extension along great circle (radians, can be negative) + - azimuth (φ): starting point angle in xy-plane from +x axis [0, 2π) + - polar (θ): starting point angle from +z axis (colatitude) [0, π] + - orient: rotation angle about the radial vector (radians) + Following right-hand rule: thumb along radial, fingers curl in +orient direction + - arc_length: extension along great circle (radians, can be negative) num_points : int Number of points to generate along the arc Returns: -------- azimuth_points : NDArray - Azimuth angles of points along the arc (radians) + Azimuth angles (φ) of points along the arc in radians polar_points : NDArray - Polar angles of points along the arc (physics convention: 0 to π, colatitude) + Polar angles (θ) of points along the arc in radians """ # Use physics convention directly: polar ∈ [0, π], where 0 = north pole, π/2 = equator # Convert starting point to Cartesian coordinates (physics convention) @@ -230,8 +270,10 @@ def quaternion_to_spherical_vector(quaternion) -> Tuple[FloatArray3, float, floa """ Convert a quaternion to a unit vector in spherical coordinates. - The quaternion is applied to a reference direction (positive Z-axis) - to get the pointing direction, which is then converted to spherical coordinates. + The quaternion is applied to a reference direction (+X axis) to get the + pointing direction, which is then converted to spherical coordinates. + + Uses ISO Physics Convention (see module docstring for details). Parameters: ----------- @@ -243,12 +285,21 @@ def quaternion_to_spherical_vector(quaternion) -> Tuple[FloatArray3, float, floa cartesian_point : FloatArray3 The Cartesian coordinates [x, y, z] of the rotated unit vector azimuth : float - Azimuth angle in radians (-π to π) + Azimuth angle (φ) in radians, range (-π, π] + Angle in xy-plane from +x axis polar : float - Polar angle in radians (0 to π, physics convention: 0 = north pole) + Polar angle (θ) in radians, range [0, π] + Angle from +z axis (colatitude) + + Notes: + ------ + Reference direction is +X axis: [1, 0, 0] + - Identity quaternion (1,0,0,0) points to +X: azimuth=0, polar=π/2 + - Quaternion (0,0,1,0) is 180° rotation about Y, points to -X: azimuth=π, polar=π/2 """ - # Reference direction: positive Z-axis (north pole) - reference_direction = array([0.0, 0.0, 1.0]) + # Reference direction: positive X-axis + # Identity quaternion (1,0,0,0) will point to +X + reference_direction = array([1.0, 0.0, 0.0]) # Apply quaternion rotation to the reference direction rotated_vector = quaternion.rotate_vector(reference_direction) @@ -266,3 +317,84 @@ def quaternion_to_spherical_vector(quaternion) -> Tuple[FloatArray3, float, floa cartesian_point = array([x, y, z]) return cartesian_point, azimuth, polar + + +def compute_spherical_arc_endpoint( + arc: UnitSphericalArc, +) -> Tuple[float, float]: + """ + Compute the endpoint of a UnitSphericalArc on a unit sphere. + + Parameters + ---------- + arc : UnitSphericalArc + The spherical arc containing: + - azimuth: Starting azimuth angle in radians (0 to 2π) + - polar: Starting polar angle in radians (colatitude/zenith angle) + measured from vertical +Z axis (0 at north pole, + π/2 at equator, π at south pole) + - arc_length: Length of the arc in radians (angular distance to travel) + - orient: Orientation/bearing of the arc in radians (-π to π) + This is the direction to travel from the start point + + Returns + ------- + Tuple[float, float] + (azimuth_end, polar_end) - endpoint in spherical coordinates + polar_end uses the same colatitude convention + + Notes + ----- + This uses spherical trigonometry to compute great circle navigation. + The polar angle uses colatitude convention (angle from +Z axis): + - polar = 0 at north pole (+z axis) + - polar = π/2 at equator + - polar = π at south pole (-z axis) + + Examples + -------- + >>> arc = UnitSphericalArc( + ... arc_length=deg2rad(45), + ... azimuth=0.0, + ... orient=0.0, + ... polar=deg2rad(90) # Starting at equator + ... ) + >>> azimuth_end, polar_end = compute_spherical_arc_endpoint(arc) + """ + # polar is already in colatitude convention (angle from +Z axis) + theta_start = arc.polar + + # Using spherical trigonometry formulas for great circle navigation: + # cos(theta_end) = cos(theta_start)*cos(arc_length) + + # sin(theta_start)*sin(arc_length)*cos(orient) + cos_theta_end = cos(theta_start) * cos(arc.arc_length) + sin(theta_start) * sin( + arc.arc_length + ) * cos(arc.orient) + theta_end = arccos(clip(cos_theta_end, -1.0, 1.0)) + + # Compute the change in azimuth using: + # sin(Δazimuth) = sin(arc_length)*sin(orient) / sin(theta_end) + # cos(Δazimuth) = (cos(arc_length) - cos(theta_start)*cos(theta_end)) / + # (sin(theta_start)*sin(theta_end)) + + if abs(sin(theta_start)) < 1e-10: + # Special case: starting point at pole + # When starting from a pole, the ending azimuth is simply the bearing direction + azimuth_end = arc.orient + elif abs(sin(theta_end)) < 1e-10: + # Special case: endpoint at pole + azimuth_end = arc.azimuth + else: + sin_delta_az = sin(arc.arc_length) * sin(arc.orient) / sin(theta_end) + cos_delta_az = (cos(arc.arc_length) - cos(theta_start) * cos(theta_end)) / ( + sin(theta_start) * sin(theta_end) + ) + delta_azimuth = arctan2(sin_delta_az, cos_delta_az) + azimuth_end = arc.azimuth + delta_azimuth + + # Normalize azimuth to [0, 2π) + azimuth_end = azimuth_end % (2 * pi) + + polar_end = theta_end + + return azimuth_end, polar_end diff --git a/src/pyMathTools/transforms/sphericalTransforms.py b/src/pyMathTools/transforms/sphericalTransforms.py index a4fe0fa..e2b0396 100644 --- a/src/pyMathTools/transforms/sphericalTransforms.py +++ b/src/pyMathTools/transforms/sphericalTransforms.py @@ -4,6 +4,14 @@ This module provides functions for converting to and from spherical coordinates, including projections and other mathematical transformations commonly used in geometry, mapping, and scientific computing. + +ISO Physics Convention +----------------------- +All functions use the ISO physics convention for spherical coordinates: +- θ (theta/polar): angle from +z axis (colatitude), range [0, π] +- φ (phi/azimuth): angle in xy-plane from +x axis, range [0, 2π) + +See pyMathTools.generators.sphericalGenerators module docstring for complete details. """ from numpy import pi @@ -14,30 +22,41 @@ def plate_carree_transform( azimuth: FloatOrNDArray, polar: FloatOrNDArray ) -> tuple[FloatOrNDArray, FloatOrNDArray]: """ - Apply plate carrée (equirectangular) projection to spherical coordinates (physics convention). + Apply plate carrée (equirectangular) projection to spherical coordinates. This is a simple cylindrical projection that preserves angular measurements, - making it suitable for robot positioning where actual angles in radians are needed. + making it suitable for applications where actual angles in radians are needed. + + Uses ISO Physics Convention (see module docstring for details). Parameters: ----------- azimuth : float or ndarray - Azimuth angle(s) in radians (0 to 2π), rotation about Z axis + Azimuth angle(s) φ in radians + Angle in xy-plane from +x axis, range [0, 2π) polar : float or ndarray - Polar angle(s) in radians (0 to π), angle from north pole + Polar angle(s) θ in radians + Angle from +z axis (colatitude), range [0, π] + - θ = 0: north pole (+z) + - θ = π/2: equator + - θ = π: south pole (-z) Returns: -------- - - x : float or ndarray - Azimuth in radians (0 to 2π) - - y : float or ndarray - Latitude in radians (-π/2 to π/2), where: - - π/2 = north pole - - 0 = equator - - -π/2 = south pole + x : float or ndarray + Projected x-coordinate = azimuth (φ) in radians [0, 2π) + y : float or ndarray + Projected y-coordinate = latitude in radians [-π/2, π/2] + - y = π/2: north pole + - y = 0: equator + - y = -π/2: south pole + + Notes: + ------ + Conversion formula: latitude = π/2 - polar """ x = azimuth - # Convert physics polar angle (theta, 0 to π) to latitude (-π/2 to π/2) - # lat = π/2 - theta + # Convert physics polar angle (θ, 0 to π) to latitude (-π/2 to π/2) + # latitude = π/2 - θ y = pi / 2 - polar return x, y diff --git a/src/pyMathToolsPlotHelpers/plotUnitSpherical.py b/src/pyMathToolsPlotHelpers/plotUnitSpherical.py index cbbaaa1..57ceb84 100644 --- a/src/pyMathToolsPlotHelpers/plotUnitSpherical.py +++ b/src/pyMathToolsPlotHelpers/plotUnitSpherical.py @@ -163,7 +163,7 @@ def plot_unit_spherical_advanced( if label is None: azimuth_deg = np.degrees(azimuth) polar_deg = np.degrees(polar) - label = f"Quat {idx + 1}: w:{quat.w:.1f}, x:{quat.x:.1f}, y:{quat.x:.1f}, z:{quat.x:.1f}" + label = f"Quat {idx + 1}: w:{quat.w:.1f}, x:{quat.x:.1f}, y:{quat.y:.1f}, z:{quat.z:.1f}" # Plot vector as an arrow from origin to the point ax.annotate( @@ -373,7 +373,7 @@ def plot_unit_spherical_polar( # Plot endpoint azimuth_deg = np.degrees(azimuth) polar_deg = np.degrees(polar) - label = f"Quat {idx + 1}: w:{quat.w:.1f}, x:{quat.x:.1f}, y:{quat.x:.1f}, z:{quat.x:.1f}" + label = f"Quat {idx + 1}: w:{quat.w:.1f}, x:{quat.x:.1f}, y:{quat.y:.1f}, z:{quat.z:.1f}" ax.scatter(azimuth, radius, s=150, marker="*", label=label, zorder=10) @@ -536,7 +536,7 @@ def plot_unit_spherical_3d( # Plot endpoint azimuth_deg = np.degrees(azimuth) polar_deg = np.degrees(polar) - label = f"Quat {idx + 1}: w:{quat.w:.1f}, x:{quat.x:.1f}, y:{quat.x:.1f}, z:{quat.x:.1f}" + label = f"Quat {idx + 1}: w:{quat.w:.1f}, x:{quat.x:.1f}, y:{quat.y:.1f}, z:{quat.z:.1f}" ax.scatter(x, y, z, s=150, marker="*", label=label, zorder=10) @@ -735,6 +735,7 @@ def demo() -> None: Quaternion.from_components(w=0.707, x=0.707, y=0, z=0), Quaternion.from_components(w=0.707, x=0, y=0.707, z=0), Quaternion.from_components(w=0.707, x=0, y=0, z=0.707), + Quaternion.from_components(w=0, x=0.577, y=-0.577, z=0.577), ] _ = plot_unit_spherical_multiplot( From 995ae5490522bfd2e800b9852645b33f051f6b96 Mon Sep 17 00:00:00 2001 From: kopecn Date: Wed, 3 Dec 2025 17:26:11 -0800 Subject: [PATCH 17/22] adding color arguements --- .pylintrc | 4 + .../plotUnitSpherical.py | 150 ++++++++++++++++-- 2 files changed, 137 insertions(+), 17 deletions(-) diff --git a/.pylintrc b/.pylintrc index 2dc5dc7..f32cab5 100644 --- a/.pylintrc +++ b/.pylintrc @@ -28,3 +28,7 @@ inlinevar-rgx=(?:(?P[a-z_][a-z0-9_]*)|(?P[a-z]+(?:[A-Z][a-z0-9]*)* [design] max-locals = 24 + +[MESSAGES CONTROL] +disable=W0621 +disable=W0102 diff --git a/src/pyMathToolsPlotHelpers/plotUnitSpherical.py b/src/pyMathToolsPlotHelpers/plotUnitSpherical.py index 57ceb84..0042003 100644 --- a/src/pyMathToolsPlotHelpers/plotUnitSpherical.py +++ b/src/pyMathToolsPlotHelpers/plotUnitSpherical.py @@ -33,6 +33,9 @@ def plot_unit_spherical_advanced( show_centers: bool = True, show_plot: bool = False, ax: Optional[Axes] = None, + arc_colors: Optional[List[str]] = ["red", "green", "blue"], + circle_colors: Optional[List[str]] = ["red", "green", "blue"], + quaternion_colors: Optional[List[str]] = ["red", "green", "blue"], ) -> Tuple[Figure, Axes]: """ Advanced version with more customization options. @@ -55,7 +58,7 @@ def plot_unit_spherical_advanced( title : str, optional Plot title colors : list, optional - List of colors for each circle + List of colors for each circle (deprecated, use circle_colors instead) labels : list, optional List of labels for each circle show_centers : bool, optional @@ -64,11 +67,29 @@ def plot_unit_spherical_advanced( Whether to show the plot (default: False, ignored if ax is provided) ax : matplotlib.axes.Axes, optional Axes to plot on. If None, creates new figure and axes + arc_colors : List[str], optional + List of colors for arcs. Defaults to ["red", "green", "blue"]. + If insufficient colors, uses matplotlib's color cycle. + circle_colors : List[str], optional + List of colors for circles. Defaults to ["red", "green", "blue"]. + If insufficient colors, uses matplotlib's color cycle. + quaternion_colors : List[str], optional + List of colors for quaternions. Defaults to ["red", "green", "blue"]. + If insufficient colors, uses matplotlib's color cycle. Returns: -------- fig, ax : matplotlib figure and axes objects """ + # Helper function to get color with fallback to matplotlib's color cycle + def get_color(color_list: Optional[List[str]], idx: int) -> str: + """Get color from list or fallback to matplotlib's color cycle.""" + if color_list is None or len(color_list) == 0: + return f"C{idx}" + if idx < len(color_list): + return color_list[idx] + return f"C{idx}" + # Track whether we created the axes or it was provided created_ax = ax is None @@ -88,8 +109,11 @@ def plot_unit_spherical_advanced( # Apply Plate Carrée projection x_points, y_points = plate_carree_transform(azimuth_points, polar_points) - # Determine color and label - color = colors[idx] if colors is not None and idx < len(colors) else None + # Determine color and label (prioritize new color parameters over deprecated colors) + if colors is not None and idx < len(colors): + color = colors[idx] + else: + color = get_color(circle_colors, idx) label = labels[idx] if labels is not None and idx < len(labels) else None # Plot the circle @@ -112,7 +136,7 @@ def plot_unit_spherical_advanced( np.degrees(y_center), s=100, marker="x", - c=color if color else "red", + c=color, linewidths=2, ) @@ -129,7 +153,7 @@ def plot_unit_spherical_advanced( x_points, y_points = plate_carree_transform(azimuth_points, polar_points) # Determine color and label - color = colors[idx] if colors is not None and idx < len(colors) else None + color = get_color(arc_colors, idx) label = labels[idx] if labels is not None and idx < len(labels) else None # Plot the arc as a line @@ -157,12 +181,10 @@ def plot_unit_spherical_advanced( x_origin, y_origin = plate_carree_transform(0, np.pi / 2) # Determine color and label - color = colors[idx] if colors is not None and idx < len(colors) else None + color = get_color(quaternion_colors, idx) label = labels[idx] if labels is not None and idx < len(labels) else None if label is None: - azimuth_deg = np.degrees(azimuth) - polar_deg = np.degrees(polar) label = f"Quat {idx + 1}: w:{quat.w:.1f}, x:{quat.x:.1f}, y:{quat.y:.1f}, z:{quat.z:.1f}" # Plot vector as an arrow from origin to the point @@ -173,7 +195,7 @@ def plot_unit_spherical_advanced( arrowprops=dict( arrowstyle="->", lw=2, - color=color if color else f"C{idx}", + color=color, alpha=0.8, ), ) @@ -225,6 +247,9 @@ def plot_unit_spherical_polar( labels: Optional[List[str]] = None, ax: Optional[Axes] = None, show_plot: bool = False, + arc_colors: Optional[List[str]] = ["red", "green", "blue"], + circle_colors: Optional[List[str]] = ["red", "green", "blue"], + quaternion_colors: Optional[List[str]] = ["red", "green", "blue"], ) -> Tuple[Figure, Axes]: """ Plot small circles and arcs on a sphere using polar projection. @@ -244,11 +269,35 @@ def plot_unit_spherical_polar( Figure size (default: (8, 8)) title : str, optional Plot title + labels : list, optional + List of labels for each element + ax : matplotlib.axes.Axes, optional + Axes to plot on. If None, creates new figure and axes + show_plot : bool, optional + Whether to show the plot (default: False, ignored if ax is provided) + arc_colors : List[str], optional + List of colors for arcs. Defaults to ["red", "green", "blue"]. + If insufficient colors, uses matplotlib's color cycle. + circle_colors : List[str], optional + List of colors for circles. Defaults to ["red", "green", "blue"]. + If insufficient colors, uses matplotlib's color cycle. + quaternion_colors : List[str], optional + List of colors for quaternions. Defaults to ["red", "green", "blue"]. + If insufficient colors, uses matplotlib's color cycle. Returns: -------- fig, ax : matplotlib figure and axes objects """ + # Helper function to get color with fallback to matplotlib's color cycle + def get_color(color_list: Optional[List[str]], idx: int) -> str: + """Get color from list or fallback to matplotlib's color cycle.""" + if color_list is None or len(color_list) == 0: + return f"C{idx}" + if idx < len(color_list): + return color_list[idx] + return f"C{idx}" + # Track whether we created the axes or it was provided created_ax = ax is None @@ -303,8 +352,11 @@ def plot_unit_spherical_polar( polar_deg = np.degrees(circle.polar) label = f"Circle {idx + 1}: A:{azimuth_deg:.1f}°, P:{polar_deg:.1f}°, R:{radius_deg:.1f}°" + # Get color for this circle + color = get_color(circle_colors, idx) + # Plot in polar coordinates using plot instead of scatter for continuous lines - ax.plot(azimuth_plot, radius_plot, linewidth=2, alpha=0.7, label=label) + ax.plot(azimuth_plot, radius_plot, linewidth=2, alpha=0.7, label=label, color=color) # Generate points for each arc if arcs is not None: @@ -342,6 +394,9 @@ def plot_unit_spherical_polar( orient_deg = np.degrees(arc.orient) label = f"Arc {idx + 1}: A:{azimuth_deg:.1f}°, P:{polar_deg:.1f}°, L:{arc_length_deg:.1f}°, O:{orient_deg:.1f}°" + # Get color for this arc + color = get_color(arc_colors, idx) + # Plot the arc ax.plot( azimuth_plot, @@ -351,6 +406,7 @@ def plot_unit_spherical_polar( label=label, marker="o", markersize=3, + color=color, ) # Plot quaternions as vectors @@ -362,12 +418,15 @@ def plot_unit_spherical_polar( # Map polar angle to radius (0 to 1) for polar plot radius = (polar + np.pi / 2) / np.pi + # Get color for this quaternion + color = get_color(quaternion_colors, idx) + # Plot vector as an arrow from origin ax.annotate( "", xy=(azimuth, radius), xytext=(0, 0), - arrowprops=dict(arrowstyle="->", lw=2, color=f"C{idx}", alpha=0.8), + arrowprops=dict(arrowstyle="->", lw=2, color=color, alpha=0.8), ) # Plot endpoint @@ -375,7 +434,7 @@ def plot_unit_spherical_polar( polar_deg = np.degrees(polar) label = f"Quat {idx + 1}: w:{quat.w:.1f}, x:{quat.x:.1f}, y:{quat.y:.1f}, z:{quat.z:.1f}" - ax.scatter(azimuth, radius, s=150, marker="*", label=label, zorder=10) + ax.scatter(azimuth, radius, s=150, marker="*", label=label, c=color, zorder=10) ax.set_title(title, fontsize=14, pad=20) @@ -404,6 +463,9 @@ def plot_unit_spherical_3d( show_legend: bool = True, show_plot: bool = False, ax: Optional[Axes] = None, + arc_colors: Optional[List[str]] = ["red", "green", "blue"], + circle_colors: Optional[List[str]] = ["red", "green", "blue"], + quaternion_colors: Optional[List[str]] = ["red", "green", "blue"], ) -> Tuple[Figure, Axes]: """ Plot small circles and arcs on a 3D unit sphere. @@ -430,11 +492,30 @@ def plot_unit_spherical_3d( Whether to show the plot (default: False, ignored if ax is provided) ax : matplotlib.axes.Axes, optional 3D axes to plot on. If None, creates new figure and axes + arc_colors : List[str], optional + List of colors for arcs. Defaults to ["red", "green", "blue"]. + If insufficient colors, uses matplotlib's color cycle. + circle_colors : List[str], optional + List of colors for circles. Defaults to ["red", "green", "blue"]. + If insufficient colors, uses matplotlib's color cycle. + quaternion_colors : List[str], optional + List of colors for quaternions. Defaults to ["red", "green", "blue"]. + If insufficient colors, uses matplotlib's color cycle. Returns: -------- fig, ax : matplotlib figure and axes objects """ + + # Helper function to get color with fallback to matplotlib's color cycle + def get_color(color_list: Optional[List[str]], idx: int) -> str: + """Get color from list or fallback to matplotlib's color cycle.""" + if color_list is None or len(color_list) == 0: + return f"C{idx}" + if idx < len(color_list): + return color_list[idx] + return f"C{idx}" + # Track whether we created the axes or it was provided created_ax = ax is None @@ -475,8 +556,11 @@ def plot_unit_spherical_3d( polar_deg = np.degrees(circle.polar) label = f"Circle {idx + 1}: A:{azimuth_deg:.1f}°, P:{polar_deg:.1f}°, R:{radius_deg:.1f}°" + # Get color for this circle + color = get_color(circle_colors, idx) + # Plot the circle - ax.scatter(x, y, z, s=10, alpha=0.8, label=label) + ax.scatter(x, y, z, s=10, alpha=0.8, label=label, c=color) # Plot each arc if arcs is not None: @@ -499,6 +583,9 @@ def plot_unit_spherical_3d( orient_deg = np.degrees(arc.orient) label = f"Arc {idx + 1}: A:{azimuth_deg:.1f}°, P:{polar_deg:.1f}°, L:{arc_length_deg:.1f}°, O:{orient_deg:.1f}°" + # Get color for this arc + color = get_color(arc_colors, idx) + # Plot the arc as a line ax.plot( x, @@ -509,6 +596,7 @@ def plot_unit_spherical_3d( label=label, marker="o", markersize=3, + color=color, ) # Plot quaternions as vectors @@ -518,6 +606,9 @@ def plot_unit_spherical_3d( cartesian, azimuth, polar = quaternion_to_spherical_vector(quat) x, y, z = cartesian + # Get color for this quaternion + color = get_color(quaternion_colors, idx) + # Plot vector from origin to the point ax.quiver( 0, @@ -528,7 +619,7 @@ def plot_unit_spherical_3d( z, length=1.0, arrow_length_ratio=0.15, - color=f"C{idx}", + color=color, linewidth=2.5, alpha=0.9, ) @@ -538,7 +629,7 @@ def plot_unit_spherical_3d( polar_deg = np.degrees(polar) label = f"Quat {idx + 1}: w:{quat.w:.1f}, x:{quat.x:.1f}, y:{quat.y:.1f}, z:{quat.z:.1f}" - ax.scatter(x, y, z, s=150, marker="*", label=label, zorder=10) + ax.scatter(x, y, z, s=150, marker="*", label=label, c=color, zorder=10) ax.set_xlabel("X") ax.set_ylabel("Y") @@ -573,6 +664,9 @@ def plot_unit_spherical_multiplot( figsize: Tuple[int, int] = (24, 6), title: str = "Small Circles on Unit Sphere - Multiple Views", show_plot: bool = False, + arc_colors: Optional[List[str]] = ["red", "green", "blue"], + circle_colors: Optional[List[str]] = ["red", "green", "blue"], + quaternion_colors: Optional[List[str]] = ["red", "green", "blue"], ) -> Figure: """ Create a multiplot showing Plate Carrée, Top View, and 3D views of spherical small circles. @@ -581,6 +675,10 @@ def plot_unit_spherical_multiplot( ----------- circles : List[UnitSphericalSmallCircle] or array-like List of UnitSphericalSmallCircle objects or array of [azimuth, polar, radius] lists + arcs : List[UnitSphericalArc], optional + List of UnitSphericalArc objects + quaternions : List[Quaternion], optional + List of Quaternion objects num_points : int, optional Number of points to generate around each circle (default: 120) figsize : tuple, optional @@ -589,6 +687,15 @@ def plot_unit_spherical_multiplot( Overall title show_plot : bool, optional Whether to show the plot (default: False) + arc_colors : List[str], optional + List of colors for arcs. Defaults to ["red", "green", "blue"]. + If insufficient colors, uses matplotlib's color cycle. + circle_colors : List[str], optional + List of colors for circles. Defaults to ["red", "green", "blue"]. + If insufficient colors, uses matplotlib's color cycle. + quaternion_colors : List[str], optional + List of colors for quaternions. Defaults to ["red", "green", "blue"]. + If insufficient colors, uses matplotlib's color cycle. Returns: -------- @@ -612,6 +719,9 @@ def plot_unit_spherical_multiplot( title="Plate Carrée Projection", show_centers=False, ax=ax1, + arc_colors=arc_colors, + circle_colors=circle_colors, + quaternion_colors=quaternion_colors, ) # Customize for multiplot ax1.set_xlabel("Azimuth (deg)", fontsize=10) @@ -626,6 +736,9 @@ def plot_unit_spherical_multiplot( num_points=num_points, title="Top View (Down Z-Axis)", ax=ax2, + arc_colors=arc_colors, + circle_colors=circle_colors, + quaternion_colors=quaternion_colors, ) # Customize for multiplot ax2.set_xlabel("Azimuth (deg)", fontsize=10) @@ -642,6 +755,9 @@ def plot_unit_spherical_multiplot( show_sphere=True, show_legend=True, ax=ax3, + arc_colors=arc_colors, + circle_colors=circle_colors, + quaternion_colors=quaternion_colors, ) # Customize for multiplot ax3.set_xlabel("X", fontsize=10) @@ -739,8 +855,8 @@ def demo() -> None: ] _ = plot_unit_spherical_multiplot( - # circles=circles, - # arcs=arcs, + circles=circles, + arcs=arcs, quaternions=quats, show_plot=True, ) From b988c0b1abcf08d3e7681189f1b8fb9975bb0e75 Mon Sep 17 00:00:00 2001 From: kopecn Date: Wed, 3 Dec 2025 17:52:44 -0800 Subject: [PATCH 18/22] humph... --- src/pyMathTools/generators/sphericalGenerators.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/pyMathTools/generators/sphericalGenerators.py b/src/pyMathTools/generators/sphericalGenerators.py index a547f61..8808281 100644 --- a/src/pyMathTools/generators/sphericalGenerators.py +++ b/src/pyMathTools/generators/sphericalGenerators.py @@ -56,6 +56,7 @@ from foundationTypes.mathTypes.UnitSphericalArc import UnitSphericalArc from pyMathTools.hints import FloatNDArray, FloatArray3 +from pyMathTools.spatial import Quaternion def generate_spherical_small_circle_points( @@ -266,7 +267,9 @@ def generate_spherical_arc_points( return azimuth_points, polar_points -def quaternion_to_spherical_vector(quaternion) -> Tuple[FloatArray3, float, float]: +def quaternion_to_spherical_vector( + quaternion: Quaternion, +) -> Tuple[FloatArray3, float, float]: """ Convert a quaternion to a unit vector in spherical coordinates. From fa4d0775f018bb6da323cfe76f660d4182e69a28 Mon Sep 17 00:00:00 2001 From: kopecn Date: Wed, 3 Dec 2025 18:45:03 -0800 Subject: [PATCH 19/22] updates... --- .../generators/sphericalGenerators.py | 139 +----------------- src/pyMathTools/pyMathTools.py | 1 - src/pyMathTools/spatial/Quaternion.py | 62 ++++++++ .../transforms/sphericalTransforms.py | 122 ++++++++++++++- .../plotUnitSpherical.py | 29 ++-- 5 files changed, 202 insertions(+), 151 deletions(-) delete mode 100644 src/pyMathTools/pyMathTools.py diff --git a/src/pyMathTools/generators/sphericalGenerators.py b/src/pyMathTools/generators/sphericalGenerators.py index 8808281..6310ba3 100644 --- a/src/pyMathTools/generators/sphericalGenerators.py +++ b/src/pyMathTools/generators/sphericalGenerators.py @@ -55,8 +55,7 @@ from foundationTypes.mathTypes.UnitSphericalSmallCircle import UnitSphericalSmallCircle from foundationTypes.mathTypes.UnitSphericalArc import UnitSphericalArc -from pyMathTools.hints import FloatNDArray, FloatArray3 -from pyMathTools.spatial import Quaternion +from pyMathTools.hints import FloatNDArray def generate_spherical_small_circle_points( @@ -265,139 +264,3 @@ def generate_spherical_arc_points( azimuth_points[i] = arctan2(y, x) # azimuth return azimuth_points, polar_points - - -def quaternion_to_spherical_vector( - quaternion: Quaternion, -) -> Tuple[FloatArray3, float, float]: - """ - Convert a quaternion to a unit vector in spherical coordinates. - - The quaternion is applied to a reference direction (+X axis) to get the - pointing direction, which is then converted to spherical coordinates. - - Uses ISO Physics Convention (see module docstring for details). - - Parameters: - ----------- - quaternion : Quaternion - A quaternion object representing a rotation - - Returns: - -------- - cartesian_point : FloatArray3 - The Cartesian coordinates [x, y, z] of the rotated unit vector - azimuth : float - Azimuth angle (φ) in radians, range (-π, π] - Angle in xy-plane from +x axis - polar : float - Polar angle (θ) in radians, range [0, π] - Angle from +z axis (colatitude) - - Notes: - ------ - Reference direction is +X axis: [1, 0, 0] - - Identity quaternion (1,0,0,0) points to +X: azimuth=0, polar=π/2 - - Quaternion (0,0,1,0) is 180° rotation about Y, points to -X: azimuth=π, polar=π/2 - """ - # Reference direction: positive X-axis - # Identity quaternion (1,0,0,0) will point to +X - reference_direction = array([1.0, 0.0, 0.0]) - - # Apply quaternion rotation to the reference direction - rotated_vector = quaternion.rotate_vector(reference_direction) - - # Normalize to ensure it's a unit vector - x, y, z = rotated_vector - r = sqrt(x**2 + y**2 + z**2) - if r > 1e-10: - x, y, z = x / r, y / r, z / r - - # Convert to spherical coordinates (physics convention) - polar = arccos(clip(z, -1, 1)) # colatitude (0 to π) - azimuth = arctan2(y, x) # azimuth (-π to π) - - cartesian_point = array([x, y, z]) - - return cartesian_point, azimuth, polar - - -def compute_spherical_arc_endpoint( - arc: UnitSphericalArc, -) -> Tuple[float, float]: - """ - Compute the endpoint of a UnitSphericalArc on a unit sphere. - - Parameters - ---------- - arc : UnitSphericalArc - The spherical arc containing: - - azimuth: Starting azimuth angle in radians (0 to 2π) - - polar: Starting polar angle in radians (colatitude/zenith angle) - measured from vertical +Z axis (0 at north pole, - π/2 at equator, π at south pole) - - arc_length: Length of the arc in radians (angular distance to travel) - - orient: Orientation/bearing of the arc in radians (-π to π) - This is the direction to travel from the start point - - Returns - ------- - Tuple[float, float] - (azimuth_end, polar_end) - endpoint in spherical coordinates - polar_end uses the same colatitude convention - - Notes - ----- - This uses spherical trigonometry to compute great circle navigation. - The polar angle uses colatitude convention (angle from +Z axis): - - polar = 0 at north pole (+z axis) - - polar = π/2 at equator - - polar = π at south pole (-z axis) - - Examples - -------- - >>> arc = UnitSphericalArc( - ... arc_length=deg2rad(45), - ... azimuth=0.0, - ... orient=0.0, - ... polar=deg2rad(90) # Starting at equator - ... ) - >>> azimuth_end, polar_end = compute_spherical_arc_endpoint(arc) - """ - # polar is already in colatitude convention (angle from +Z axis) - theta_start = arc.polar - - # Using spherical trigonometry formulas for great circle navigation: - # cos(theta_end) = cos(theta_start)*cos(arc_length) + - # sin(theta_start)*sin(arc_length)*cos(orient) - cos_theta_end = cos(theta_start) * cos(arc.arc_length) + sin(theta_start) * sin( - arc.arc_length - ) * cos(arc.orient) - theta_end = arccos(clip(cos_theta_end, -1.0, 1.0)) - - # Compute the change in azimuth using: - # sin(Δazimuth) = sin(arc_length)*sin(orient) / sin(theta_end) - # cos(Δazimuth) = (cos(arc_length) - cos(theta_start)*cos(theta_end)) / - # (sin(theta_start)*sin(theta_end)) - - if abs(sin(theta_start)) < 1e-10: - # Special case: starting point at pole - # When starting from a pole, the ending azimuth is simply the bearing direction - azimuth_end = arc.orient - elif abs(sin(theta_end)) < 1e-10: - # Special case: endpoint at pole - azimuth_end = arc.azimuth - else: - sin_delta_az = sin(arc.arc_length) * sin(arc.orient) / sin(theta_end) - cos_delta_az = (cos(arc.arc_length) - cos(theta_start) * cos(theta_end)) / ( - sin(theta_start) * sin(theta_end) - ) - delta_azimuth = arctan2(sin_delta_az, cos_delta_az) - azimuth_end = arc.azimuth + delta_azimuth - - # Normalize azimuth to [0, 2π) - azimuth_end = azimuth_end % (2 * pi) - - polar_end = theta_end - - return azimuth_end, polar_end diff --git a/src/pyMathTools/pyMathTools.py b/src/pyMathTools/pyMathTools.py deleted file mode 100644 index dd0b80e..0000000 --- a/src/pyMathTools/pyMathTools.py +++ /dev/null @@ -1 +0,0 @@ -"""Main module.""" diff --git a/src/pyMathTools/spatial/Quaternion.py b/src/pyMathTools/spatial/Quaternion.py index 0294127..a946ae5 100644 --- a/src/pyMathTools/spatial/Quaternion.py +++ b/src/pyMathTools/spatial/Quaternion.py @@ -295,6 +295,68 @@ def axis(self) -> FloatArray3: return np.array([0.0, 0.0, 1.0]) # Default axis for identity return vect / norm + @property + def vector_spherical(self) -> tuple[float, float]: + """Return the quaternion's pointing direction in spherical coordinates (ISO physics convention). + + The quaternion is applied to a reference direction (+X axis) to get the + pointing direction, which is then converted to spherical coordinates. + + Returns: + tuple[azimuth, polar] where: + - azimuth (φ): angle in xy-plane from +x axis, range (-π, π] + - polar (θ): angle from +z axis (colatitude), range [0, π] + + Notes: + Reference direction is +X axis: [1, 0, 0] + - Identity quaternion (1,0,0,0) points to +X: azimuth=0, polar=π/2 + """ + # Reference direction: positive X-axis + reference_direction = np.array([1.0, 0.0, 0.0]) + + # Apply quaternion rotation to the reference direction + rotated_vector = self.rotate_vector(reference_direction) + + # Normalize to ensure it's a unit vector + x, y, z = rotated_vector + r = np.sqrt(x**2 + y**2 + z**2) + if r > 1e-10: + x, y, z = x / r, y / r, z / r + + # Convert to spherical coordinates (physics convention) + polar = float(np.arccos(np.clip(z, -1, 1))) # colatitude (0 to π) + azimuth = float(np.arctan2(y, x)) # azimuth (-π to π) + + return (azimuth, polar) + + @property + def vector_cartesian(self) -> tuple[float, float, float]: + """Return the quaternion's pointing direction as a cartesian unit vector. + + The quaternion is applied to a reference direction (+X axis) to get the + pointing direction as a unit vector in cartesian coordinates. + + Returns: + tuple[x, y, z] - unit vector components + + Notes: + Reference direction is +X axis: [1, 0, 0] + - Identity quaternion (1,0,0,0) points to +X: (1, 0, 0) + """ + # Reference direction: positive X-axis + reference_direction = np.array([1.0, 0.0, 0.0]) + + # Apply quaternion rotation to the reference direction + rotated_vector = self.rotate_vector(reference_direction) + + # Normalize to ensure it's a unit vector + x, y, z = rotated_vector + r = np.sqrt(x**2 + y**2 + z**2) + if r > 1e-10: + x, y, z = x / r, y / r, z / r + + return (float(x), float(y), float(z)) + # MARK: - Array Conversion Methods def as_float_array(self) -> FloatArray4: diff --git a/src/pyMathTools/transforms/sphericalTransforms.py b/src/pyMathTools/transforms/sphericalTransforms.py index e2b0396..7607e54 100644 --- a/src/pyMathTools/transforms/sphericalTransforms.py +++ b/src/pyMathTools/transforms/sphericalTransforms.py @@ -14,8 +14,21 @@ See pyMathTools.generators.sphericalGenerators module docstring for complete details. """ -from numpy import pi +from typing import Tuple from pyMathTools.hints import FloatOrNDArray +from numpy import float64 +from numpy.typing import NDArray +from numpy import ( + clip, + pi, + arctan2, + cos, + sin, + arccos, +) +from numpy.linalg import norm + +from foundationTypes.mathTypes.UnitSphericalArc import UnitSphericalArc def plate_carree_transform( @@ -60,3 +73,110 @@ def plate_carree_transform( # latitude = π/2 - θ y = pi / 2 - polar return x, y + + +def cartesian_to_spherical(point: NDArray[float64]) -> tuple: + """ + Convert Cartesian coordinates to spherical coordinates. + + Uses ISO physics convention: + - polar (theta): angle from +z axis (colatitude), range [0, π] + - azimuth (phi): angle in xy-plane from +x axis, range (-π, π] + + Args: + point: Cartesian coordinates [x, y, z] + + Returns: + Tuple of (azimuth, polar) in radians + """ + point = point / norm(point) + x, y, z = point + + # Azimuth: angle in xy-plane from +x axis + azimuth = arctan2(y, x) + + # Polar: angle from +z axis (colatitude) + polar = arccos(clip(z, -1.0, 1.0)) + + return azimuth, polar + + +def compute_spherical_arc_endpoint( + arc: UnitSphericalArc, +) -> Tuple[float, float]: + """ + Compute the endpoint of a UnitSphericalArc on a unit sphere. + + Parameters + ---------- + arc : UnitSphericalArc + The spherical arc containing: + - azimuth: Starting azimuth angle in radians (0 to 2π) + - polar: Starting polar angle in radians (colatitude/zenith angle) + measured from vertical +Z axis (0 at north pole, + π/2 at equator, π at south pole) + - arc_length: Length of the arc in radians (angular distance to travel) + - orient: Orientation/bearing of the arc in radians (-π to π) + This is the direction to travel from the start point + + Returns + ------- + Tuple[float, float] + (azimuth_end, polar_end) - endpoint in spherical coordinates + polar_end uses the same colatitude convention + + Notes + ----- + This uses spherical trigonometry to compute great circle navigation. + The polar angle uses colatitude convention (angle from +Z axis): + - polar = 0 at north pole (+z axis) + - polar = π/2 at equator + - polar = π at south pole (-z axis) + + Examples + -------- + >>> arc = UnitSphericalArc( + ... arc_length=deg2rad(45), + ... azimuth=0.0, + ... orient=0.0, + ... polar=deg2rad(90) # Starting at equator + ... ) + >>> azimuth_end, polar_end = compute_spherical_arc_endpoint(arc) + """ + # polar is already in colatitude convention (angle from +Z axis) + theta_start = arc.polar + + # Using spherical trigonometry formulas for great circle navigation: + # cos(theta_end) = cos(theta_start)*cos(arc_length) + + # sin(theta_start)*sin(arc_length)*cos(orient) + cos_theta_end = cos(theta_start) * cos(arc.arc_length) + sin(theta_start) * sin( + arc.arc_length + ) * cos(arc.orient) + theta_end = arccos(clip(cos_theta_end, -1.0, 1.0)) + + # Compute the change in azimuth using: + # sin(Δazimuth) = sin(arc_length)*sin(orient) / sin(theta_end) + # cos(Δazimuth) = (cos(arc_length) - cos(theta_start)*cos(theta_end)) / + # (sin(theta_start)*sin(theta_end)) + + if abs(sin(theta_start)) < 1e-10: + # Special case: starting point at pole + # When starting from a pole, the ending azimuth is simply the bearing direction + azimuth_end = arc.orient + elif abs(sin(theta_end)) < 1e-10: + # Special case: endpoint at pole + azimuth_end = arc.azimuth + else: + sin_delta_az = sin(arc.arc_length) * sin(arc.orient) / sin(theta_end) + cos_delta_az = (cos(arc.arc_length) - cos(theta_start) * cos(theta_end)) / ( + sin(theta_start) * sin(theta_end) + ) + delta_azimuth = arctan2(sin_delta_az, cos_delta_az) + azimuth_end = arc.azimuth + delta_azimuth + + # Normalize azimuth to [0, 2π) + azimuth_end = azimuth_end % (2 * pi) + + polar_end = theta_end + + return azimuth_end, polar_end diff --git a/src/pyMathToolsPlotHelpers/plotUnitSpherical.py b/src/pyMathToolsPlotHelpers/plotUnitSpherical.py index 0042003..0671455 100644 --- a/src/pyMathToolsPlotHelpers/plotUnitSpherical.py +++ b/src/pyMathToolsPlotHelpers/plotUnitSpherical.py @@ -12,7 +12,6 @@ from pyMathTools.generators.sphericalGenerators import ( generate_spherical_small_circle_points, generate_spherical_arc_points, - quaternion_to_spherical_vector, ) deg45 = np.deg2rad(45) @@ -81,6 +80,7 @@ def plot_unit_spherical_advanced( -------- fig, ax : matplotlib figure and axes objects """ + # Helper function to get color with fallback to matplotlib's color cycle def get_color(color_list: Optional[List[str]], idx: int) -> str: """Get color from list or fallback to matplotlib's color cycle.""" @@ -172,7 +172,7 @@ def get_color(color_list: Optional[List[str]], idx: int) -> str: if quaternions is not None: for idx, quat in enumerate(quaternions): # Convert quaternion to spherical coordinates - _, azimuth, polar = quaternion_to_spherical_vector(quat) + azimuth, polar = quat.vector_spherical # Apply Plate Carrée projection to the endpoint x_end, y_end = plate_carree_transform(azimuth, polar) @@ -289,6 +289,7 @@ def plot_unit_spherical_polar( -------- fig, ax : matplotlib figure and axes objects """ + # Helper function to get color with fallback to matplotlib's color cycle def get_color(color_list: Optional[List[str]], idx: int) -> str: """Get color from list or fallback to matplotlib's color cycle.""" @@ -356,7 +357,14 @@ def get_color(color_list: Optional[List[str]], idx: int) -> str: color = get_color(circle_colors, idx) # Plot in polar coordinates using plot instead of scatter for continuous lines - ax.plot(azimuth_plot, radius_plot, linewidth=2, alpha=0.7, label=label, color=color) + ax.plot( + azimuth_plot, + radius_plot, + linewidth=2, + alpha=0.7, + label=label, + color=color, + ) # Generate points for each arc if arcs is not None: @@ -413,7 +421,7 @@ def get_color(color_list: Optional[List[str]], idx: int) -> str: if quaternions is not None: for idx, quat in enumerate(quaternions): # Convert quaternion to spherical coordinates - _, azimuth, polar = quaternion_to_spherical_vector(quat) + azimuth, polar = quat.vector_spherical # Map polar angle to radius (0 to 1) for polar plot radius = (polar + np.pi / 2) / np.pi @@ -434,7 +442,9 @@ def get_color(color_list: Optional[List[str]], idx: int) -> str: polar_deg = np.degrees(polar) label = f"Quat {idx + 1}: w:{quat.w:.1f}, x:{quat.x:.1f}, y:{quat.y:.1f}, z:{quat.z:.1f}" - ax.scatter(azimuth, radius, s=150, marker="*", label=label, c=color, zorder=10) + ax.scatter( + azimuth, radius, s=150, marker="*", label=label, c=color, zorder=10 + ) ax.set_title(title, fontsize=14, pad=20) @@ -603,8 +613,7 @@ def get_color(color_list: Optional[List[str]], idx: int) -> str: if quaternions is not None: for idx, quat in enumerate(quaternions): # Convert quaternion to Cartesian coordinates - cartesian, azimuth, polar = quaternion_to_spherical_vector(quat) - x, y, z = cartesian + x, y, z = quat.vector_cartesian # Get color for this quaternion color = get_color(quaternion_colors, idx) @@ -625,8 +634,6 @@ def get_color(color_list: Optional[List[str]], idx: int) -> str: ) # Plot endpoint - azimuth_deg = np.degrees(azimuth) - polar_deg = np.degrees(polar) label = f"Quat {idx + 1}: w:{quat.w:.1f}, x:{quat.x:.1f}, y:{quat.y:.1f}, z:{quat.z:.1f}" ax.scatter(x, y, z, s=150, marker="*", label=label, c=color, zorder=10) @@ -855,8 +862,8 @@ def demo() -> None: ] _ = plot_unit_spherical_multiplot( - circles=circles, - arcs=arcs, + # circles=circles, + # arcs=arcs, quaternions=quats, show_plot=True, ) From e0f6234cdb626f35d011421fd5669e4b2fa3e4ef Mon Sep 17 00:00:00 2001 From: kopecn Date: Wed, 3 Dec 2025 19:17:24 -0800 Subject: [PATCH 20/22] moving things around and cleaning up --- examples/sphericalPlotting/plotArcs.py | 45 ++++++ .../plotQuatUnitCircles.py | 2 +- .../constructUnitSphericalSmallCircle.py | 61 -------- src/pyMathTools/generators/__init__.py | 0 .../{constructors => spherical}/__init__.py | 0 src/pyMathTools/spherical/constructors.py | 148 ++++++++++++++++++ .../sphericalGenerators.py | 0 .../sphericalTransforms.py | 0 src/pyMathTools/transforms/__init__.py | 0 .../plotUnitSpherical.py | 4 +- 10 files changed, 196 insertions(+), 64 deletions(-) create mode 100644 examples/sphericalPlotting/plotArcs.py rename examples/{ => sphericalPlotting}/plotQuatUnitCircles.py (96%) delete mode 100644 src/pyMathTools/constructors/constructUnitSphericalSmallCircle.py delete mode 100644 src/pyMathTools/generators/__init__.py rename src/pyMathTools/{constructors => spherical}/__init__.py (100%) create mode 100644 src/pyMathTools/spherical/constructors.py rename src/pyMathTools/{generators => spherical}/sphericalGenerators.py (100%) rename src/pyMathTools/{transforms => spherical}/sphericalTransforms.py (100%) delete mode 100644 src/pyMathTools/transforms/__init__.py diff --git a/examples/sphericalPlotting/plotArcs.py b/examples/sphericalPlotting/plotArcs.py new file mode 100644 index 0000000..991acb2 --- /dev/null +++ b/examples/sphericalPlotting/plotArcs.py @@ -0,0 +1,45 @@ +from typing import List + +from foundationTypes.mathTypes.UnitSphericalArc import UnitSphericalArc + +from pyMathToolsPlotHelpers.plotUnitSpherical import ( + deg45, + deg90, + deg180, + deg22_5, + plot_unit_spherical_multiplot, +) + +from pyMathTools.spherical.constructors import ( + arc_from_two_points, +) +from pyMathTools.spherical.sphericalTransforms import compute_spherical_arc_endpoint + + +def demo() -> None: + """Demonstration of multiplot view showing spherical small circles in three different projections.""" + + arc1 = UnitSphericalArc(orient=0, azimuth=0, polar=0, arc_length=deg90) + arc2 = UnitSphericalArc(orient=deg90, azimuth=0, polar=deg90, arc_length=deg90) + az, po = compute_spherical_arc_endpoint(arc2) + arc3 = arc_from_two_points( + azimuth1=az, + polar1=po, + azimuth2=arc1.azimuth, + polar2=arc1.polar, + ) + + arcs: List[UnitSphericalArc] = [ + arc1, + arc2, + arc3, + ] + + _ = plot_unit_spherical_multiplot( + arcs=arcs, + show_plot=True, + ) + + +if __name__ == "__main__": + demo() diff --git a/examples/plotQuatUnitCircles.py b/examples/sphericalPlotting/plotQuatUnitCircles.py similarity index 96% rename from examples/plotQuatUnitCircles.py rename to examples/sphericalPlotting/plotQuatUnitCircles.py index d59374a..6f2b593 100644 --- a/examples/plotQuatUnitCircles.py +++ b/examples/sphericalPlotting/plotQuatUnitCircles.py @@ -28,7 +28,7 @@ from numpy import pi import quaternion -from pyMathTools.constructors.constructUnitSphericalSmallCircle import ( +from pyMathTools.spherical.constructors import ( quat_to_unitSphericalSmallCircle, ) diff --git a/src/pyMathTools/constructors/constructUnitSphericalSmallCircle.py b/src/pyMathTools/constructors/constructUnitSphericalSmallCircle.py deleted file mode 100644 index e657853..0000000 --- a/src/pyMathTools/constructors/constructUnitSphericalSmallCircle.py +++ /dev/null @@ -1,61 +0,0 @@ -""" -Constructor functions for UnitSphericalSmallCircle objects. - -This module provides utilities for creating UnitSphericalSmallCircle instances -from various representations, using ISO Physics Convention for spherical coordinates. -""" - -import numpy as np -from numpy import atan2, acos, pi - -from foundationTypes.mathTypes.UnitSphericalSmallCircle import UnitSphericalSmallCircle - -from pyMathTools.spatial.Quaternion import Quaternion - - -def quat_to_unitSphericalSmallCircle( - q: Quaternion, radius_angle: float = pi / 4 -) -> UnitSphericalSmallCircle: - """ - Convert a quaternion's vector part to a UnitSphericalSmallCircle. - - Uses ISO Physics Convention for spherical coordinates: - - φ (azimuth): angle in xy-plane from +x axis, range [0, 2π) - - θ (polar): angle from +z axis (colatitude), range [0, π] - - Parameters: - ----------- - q : Quaternion - Quaternion whose vector part (x, y, z) defines the center direction - radius_angle : float - Angular radius of the small circle (default: π/4) - - Returns: - -------- - UnitSphericalSmallCircle - Small circle centered at the direction defined by the quaternion's vector part - - Notes: - ------ - - If vector part is zero, defaults to north pole (azimuth=0, polar=π/2) - - Uses standard spherical-to-Cartesian formulas: - x = r*sin(θ)*cos(φ) - y = r*sin(θ)*sin(φ) - z = r*cos(θ) - """ - # Extract vector part (this represents a direction in Cartesian space) - r = np.sqrt(q.x**2 + q.y**2 + q.z**2) - - if r == 0: - # Default to equator at +x axis if vector part is zero - azimuth = 0.0 - polar = pi / 2 - else: - # Convert Cartesian (x,y,z) to spherical (azimuth, polar) - # ISO Physics Convention: - azimuth = atan2(q.y, q.x) # φ = atan2(y, x) - angle in xy-plane from +x - polar = acos(q.z / r) # θ = acos(z/r) - angle from +z axis - - return UnitSphericalSmallCircle( - azimuth=azimuth, polar=polar, radius_angle=radius_angle - ) diff --git a/src/pyMathTools/generators/__init__.py b/src/pyMathTools/generators/__init__.py deleted file mode 100644 index e69de29..0000000 diff --git a/src/pyMathTools/constructors/__init__.py b/src/pyMathTools/spherical/__init__.py similarity index 100% rename from src/pyMathTools/constructors/__init__.py rename to src/pyMathTools/spherical/__init__.py diff --git a/src/pyMathTools/spherical/constructors.py b/src/pyMathTools/spherical/constructors.py new file mode 100644 index 0000000..84a59c2 --- /dev/null +++ b/src/pyMathTools/spherical/constructors.py @@ -0,0 +1,148 @@ +""" +Constructor functions for UnitSphericalSmallCircle objects. + +This module provides utilities for creating UnitSphericalSmallCircle instances +from various representations, using ISO Physics Convention for spherical coordinates. +""" + +import numpy as np +from numpy import atan2, acos, pi + +from foundationTypes.mathTypes.UnitSphericalSmallCircle import UnitSphericalSmallCircle +from foundationTypes.mathTypes.UnitSphericalArc import UnitSphericalArc + +from pyMathTools.spatial.Quaternion import Quaternion + + +def quat_to_unitSphericalSmallCircle( + q: Quaternion, radius_angle: float = pi / 4 +) -> UnitSphericalSmallCircle: + """ + Convert a quaternion's vector part to a UnitSphericalSmallCircle. + + Uses ISO Physics Convention for spherical coordinates: + - φ (azimuth): angle in xy-plane from +x axis, range [0, 2π) + - θ (polar): angle from +z axis (colatitude), range [0, π] + + Parameters: + ----------- + q : Quaternion + Quaternion whose vector part (x, y, z) defines the center direction + radius_angle : float + Angular radius of the small circle (default: π/4) + + Returns: + -------- + UnitSphericalSmallCircle + Small circle centered at the direction defined by the quaternion's vector part + + Notes: + ------ + - If vector part is zero, defaults to north pole (azimuth=0, polar=π/2) + - Uses standard spherical-to-Cartesian formulas: + x = r*sin(θ)*cos(φ) + y = r*sin(θ)*sin(φ) + z = r*cos(θ) + """ + # Extract vector part (this represents a direction in Cartesian space) + r = np.sqrt(q.x**2 + q.y**2 + q.z**2) + + if r == 0: + # Default to equator at +x axis if vector part is zero + azimuth = 0.0 + polar = pi / 2 + else: + # Convert Cartesian (x,y,z) to spherical (azimuth, polar) + # ISO Physics Convention: + azimuth = atan2(q.y, q.x) # φ = atan2(y, x) - angle in xy-plane from +x + polar = acos(q.z / r) # θ = acos(z/r) - angle from +z axis + + return UnitSphericalSmallCircle( + azimuth=azimuth, polar=polar, radius_angle=radius_angle + ) + + +def arc_from_two_points( + azimuth1: float, + polar1: float, + azimuth2: float, + polar2: float, + isPositive: bool = True, +) -> UnitSphericalArc: + """ + Create a UnitSphericalArc connecting two points on a unit sphere. + + Uses ISO Physics Convention for spherical coordinates: + - φ (azimuth): angle in xy-plane from +x axis + - θ (polar): angle from +z axis (colatitude), range [0, π] + + The arc follows a great circle path between the two points. The direction of + travel is determined by the isPositive parameter. + + Parameters: + ----------- + azimuth1 : float + Starting point azimuth angle (φ₁) in radians + polar1 : float + Starting point polar angle (θ₁) in radians (colatitude from +z) + azimuth2 : float + Ending point azimuth angle (φ₂) in radians + polar2 : float + Ending point polar angle (θ₂) in radians (colatitude from +z) + isPositive : bool, optional + If True (default), creates arc with positive arc_length + If False, creates arc with negative arc_length (opposite direction) + + Returns: + -------- + UnitSphericalArc + Arc starting at (azimuth1, polar1) that reaches toward (azimuth2, polar2) + + Notes: + ------ + Uses spherical trigonometry formulas: + - arc_length: Great circle distance using spherical law of cosines + - orient: Initial bearing computed relative to north direction + + Examples: + --------- + >>> # Arc from north pole to equator at 0° azimuth + >>> arc = arc_from_two_points(0, 0, 0, np.pi/2, isPositive=True) + """ + # Great circle distance using spherical law of cosines + # In ISO physics convention: polar = θ (colatitude from +z axis) + # γ = arccos(cos(θ1)*cos(θ2) + sin(θ1)*sin(θ2)*cos(φ2 - φ1)) + + cos_theta1 = np.cos(polar1) + sin_theta1 = np.sin(polar1) + cos_theta2 = np.cos(polar2) + sin_theta2 = np.sin(polar2) + + delta_phi = azimuth2 - azimuth1 + cos_delta_phi = np.cos(delta_phi) + + # Arc length (great circle distance) + cos_arc_length = cos_theta1 * cos_theta2 + sin_theta1 * sin_theta2 * cos_delta_phi + arc_length = np.arccos(np.clip(cos_arc_length, -1.0, 1.0)) + + # Initial bearing from point 1 to point 2 (relative to north) + # This gives the orient parameter + # bearing = atan2(sin(Δφ)*sin(θ2), cos(θ1)*sin(θ2) - sin(θ1)*cos(θ2)*cos(Δφ)) + + sin_delta_phi = np.sin(delta_phi) + + y_component = sin_delta_phi * sin_theta2 + x_component = cos_theta1 * sin_theta2 - sin_theta1 * cos_theta2 * cos_delta_phi + + orient = atan2(y_component, x_component) + + # Apply isPositive parameter + if not isPositive: + arc_length = -arc_length + + return UnitSphericalArc( + azimuth=azimuth1, + polar=polar1, + orient=orient, + arc_length=arc_length, + ) diff --git a/src/pyMathTools/generators/sphericalGenerators.py b/src/pyMathTools/spherical/sphericalGenerators.py similarity index 100% rename from src/pyMathTools/generators/sphericalGenerators.py rename to src/pyMathTools/spherical/sphericalGenerators.py diff --git a/src/pyMathTools/transforms/sphericalTransforms.py b/src/pyMathTools/spherical/sphericalTransforms.py similarity index 100% rename from src/pyMathTools/transforms/sphericalTransforms.py rename to src/pyMathTools/spherical/sphericalTransforms.py diff --git a/src/pyMathTools/transforms/__init__.py b/src/pyMathTools/transforms/__init__.py deleted file mode 100644 index e69de29..0000000 diff --git a/src/pyMathToolsPlotHelpers/plotUnitSpherical.py b/src/pyMathToolsPlotHelpers/plotUnitSpherical.py index 0671455..2782940 100644 --- a/src/pyMathToolsPlotHelpers/plotUnitSpherical.py +++ b/src/pyMathToolsPlotHelpers/plotUnitSpherical.py @@ -8,8 +8,8 @@ from foundationTypes.mathTypes.UnitSphericalArc import UnitSphericalArc from pyMathTools.spatial.Quaternion import Quaternion -from pyMathTools.transforms.sphericalTransforms import plate_carree_transform -from pyMathTools.generators.sphericalGenerators import ( +from pyMathTools.spherical.sphericalTransforms import plate_carree_transform +from pyMathTools.spherical.sphericalGenerators import ( generate_spherical_small_circle_points, generate_spherical_arc_points, ) From 79fdfa69fa158e559aeb4235601f9c5f5773340a Mon Sep 17 00:00:00 2001 From: kopecn Date: Wed, 3 Dec 2025 19:26:42 -0800 Subject: [PATCH 21/22] adding inverse cartesian_to_spherical --- .../spherical/sphericalTransforms.py | 28 +++++++++++++++++++ 1 file changed, 28 insertions(+) diff --git a/src/pyMathTools/spherical/sphericalTransforms.py b/src/pyMathTools/spherical/sphericalTransforms.py index 7607e54..b202fbd 100644 --- a/src/pyMathTools/spherical/sphericalTransforms.py +++ b/src/pyMathTools/spherical/sphericalTransforms.py @@ -75,6 +75,34 @@ def plate_carree_transform( return x, y +def spherical_to_cartesian(azimuth: float, polar: float) -> tuple[float, float, float]: + """ + Convert spherical coordinates to Cartesian coordinates on the unit sphere. + + Uses ISO physics convention: + - polar (theta): angle from +z axis (colatitude), range [0, π] + - azimuth (phi): angle in xy-plane from +x axis, range (-π, π] or [0, 2π) + + Args: + azimuth: Azimuth angle (φ) in radians - angle in xy-plane from +x axis + polar: Polar angle (θ) in radians - angle from +z axis (colatitude) + + Returns: + Tuple with Cartesian coordinates [x, y, z] on unit sphere + + Notes: + Conversion formulas (ISO physics convention): + - x = sin(θ) * cos(φ) + - y = sin(θ) * sin(φ) + - z = cos(θ) + """ + x = sin(polar) * cos(azimuth) + y = sin(polar) * sin(azimuth) + z = cos(polar) + + return (x, y, z) + + def cartesian_to_spherical(point: NDArray[float64]) -> tuple: """ Convert Cartesian coordinates to spherical coordinates. From b524277e6753177a429c142cee9b6021e1f9de22 Mon Sep 17 00:00:00 2001 From: kopecn Date: Wed, 3 Dec 2025 20:16:31 -0800 Subject: [PATCH 22/22] that looks good. --- .../sphericalPlotting/plotQuatUnitCircles.py | 37 +++-- src/pyMathTools/spatial/Quaternion.py | 143 ++++++++++++++++++ src/pyMathTools/spherical/constructors.py | 51 ------- .../plotUnitSpherical.py | 12 +- 4 files changed, 176 insertions(+), 67 deletions(-) diff --git a/examples/sphericalPlotting/plotQuatUnitCircles.py b/examples/sphericalPlotting/plotQuatUnitCircles.py index 6f2b593..9f12860 100644 --- a/examples/sphericalPlotting/plotQuatUnitCircles.py +++ b/examples/sphericalPlotting/plotQuatUnitCircles.py @@ -25,12 +25,12 @@ """ +from typing import List from numpy import pi -import quaternion -from pyMathTools.spherical.constructors import ( - quat_to_unitSphericalSmallCircle, -) +from foundationTypes.mathTypes.UnitSphericalSmallCircle import UnitSphericalSmallCircle + +from pyMathTools.spatial.Quaternion import Quaternion from pyMathToolsPlotHelpers.plotUnitSpherical import ( plot_unit_spherical_multiplot, @@ -48,13 +48,30 @@ def main(): All circles use a radius angle of pi/8 except y and z, which use the default. """ + + quat4 = Quaternion.from_unit_x_to_vector([1, 0, 0]) + quat5 = Quaternion.from_unit_x_to_vector([0, 1, 0]) + quat6 = Quaternion.from_unit_x_to_vector([0, 0, 1]) + + quat1 = Quaternion.from_components(0, 1, 0, 0) + quat2 = Quaternion.from_components(0, 0, 1, 0) + quat3 = Quaternion.from_components(0, 0, 0, 1) + + circles: List[UnitSphericalSmallCircle] = [ + quat1.to_unitSphericalSmallCircle(radius_angle=pi / 8), + quat2.to_unitSphericalSmallCircle(radius_angle=pi / 8), + quat3.to_unitSphericalSmallCircle(), + quat4.to_unitSphericalSmallCircle(radius_angle=pi / 16), + quat5.to_unitSphericalSmallCircle(radius_angle=pi / 16), + quat6.to_unitSphericalSmallCircle(radius_angle=pi / 16), + ] + + quats: List[Quaternion] = [quat1, quat2, quat3, quat4, quat5, quat6] + + # from_vector_part _ = plot_unit_spherical_multiplot( - circles=[ - quat_to_unitSphericalSmallCircle(q=quaternion.one, radius_angle=pi / 8), - quat_to_unitSphericalSmallCircle(q=quaternion.x, radius_angle=pi / 8), - quat_to_unitSphericalSmallCircle(q=quaternion.y), - quat_to_unitSphericalSmallCircle(q=quaternion.z), - ], + circles=circles, + quaternions=quats, show_plot=True, ) diff --git a/src/pyMathTools/spatial/Quaternion.py b/src/pyMathTools/spatial/Quaternion.py index a946ae5..966ee3c 100644 --- a/src/pyMathTools/spatial/Quaternion.py +++ b/src/pyMathTools/spatial/Quaternion.py @@ -34,6 +34,8 @@ RotationMatrix, FloatOrQuaternion, ) +from foundationTypes.mathTypes.UnitSphericalSmallCircle import UnitSphericalSmallCircle + T = TypeVar("T", bound="Quaternion") @@ -620,6 +622,140 @@ def from_axis_angle(cls: Type[T], axis: FloatArray3, angle: float) -> T: z=axis_normalized[2] * sin_half, ) + @classmethod + def _from_unit_direction_to_vector( + cls: Type[T], + start_direction: FloatArray3, + target_vector: FloatArray3, + perpendicular_axis: FloatArray3, + ) -> T: + """Helper method to create a quaternion rotating from a unit direction to target. + + Args: + start_direction: Starting unit direction vector + target_vector: Target direction vector (will be normalized) + perpendicular_axis: Axis to use for 180° rotation case + + Returns: + Unit quaternion representing the rotation + """ + # Normalize target vector + target = np.array(target_vector, dtype=float) + target_norm = np.linalg.norm(target) + if target_norm < 1e-10: + # Zero vector, return identity + return cls.identity() + target = target / target_norm + + # Compute dot product + dot = np.dot(start_direction, target) + + # Check for special cases + if dot > 0.9999: + # Vectors are parallel (same direction) + return cls.identity() + elif dot < -0.9999: + # Vectors are opposite (180° rotation needed) + # Use the provided perpendicular axis + perp = np.array(perpendicular_axis, dtype=float) + perp = perp / np.linalg.norm(perp) + return cls.from_components(w=0.0, x=perp[0], y=perp[1], z=perp[2]) + + # General case: compute rotation axis and angle + # Rotation axis is perpendicular to both vectors + axis = np.cross(start_direction, target) + axis = axis / np.linalg.norm(axis) + + # Compute rotation angle + angle = np.arccos(np.clip(dot, -1.0, 1.0)) + + # Create quaternion from axis-angle + return cls.from_axis_angle(axis, angle) + + @classmethod + def from_unit_x_to_vector(cls: Type[T], target_vector: FloatArray3) -> T: + """Create a quaternion that rotates the +X axis to point toward target_vector. + + This constructor computes the rotation that transforms the unit X direction + [1, 0, 0] to align with the given target vector. + + Args: + target_vector: Target direction vector (will be normalized) + + Returns: + Unit quaternion representing the rotation from +X to target_vector + + Notes: + - If target_vector points in +X direction, returns identity quaternion + - If target_vector points in -X direction, returns 180° rotation about +Z axis + - For all other directions, uses the shortest rotation path + + Examples: + >>> # Rotate +X to point toward +Y + >>> q = Quaternion.from_unit_x_to_vector([0, 1, 0]) + >>> # Rotate +X to point toward +Z + >>> q = Quaternion.from_unit_x_to_vector([0, 0, 1]) + """ + start = np.array([1.0, 0.0, 0.0]) + perpendicular = np.array([0.0, 0.0, 1.0]) # +Z for 180° case + return cls._from_unit_direction_to_vector(start, target_vector, perpendicular) + + @classmethod + def from_unit_y_to_vector(cls: Type[T], target_vector: FloatArray3) -> T: + """Create a quaternion that rotates the +Y axis to point toward target_vector. + + This constructor computes the rotation that transforms the unit Y direction + [0, 1, 0] to align with the given target vector. + + Args: + target_vector: Target direction vector (will be normalized) + + Returns: + Unit quaternion representing the rotation from +Y to target_vector + + Notes: + - If target_vector points in +Y direction, returns identity quaternion + - If target_vector points in -Y direction, returns 180° rotation about +Z axis + - For all other directions, uses the shortest rotation path + + Examples: + >>> # Rotate +Y to point toward +X + >>> q = Quaternion.from_unit_y_to_vector([1, 0, 0]) + >>> # Rotate +Y to point toward +Z + >>> q = Quaternion.from_unit_y_to_vector([0, 0, 1]) + """ + start = np.array([0.0, 1.0, 0.0]) + perpendicular = np.array([0.0, 0.0, 1.0]) # +Z for 180° case + return cls._from_unit_direction_to_vector(start, target_vector, perpendicular) + + @classmethod + def from_unit_z_to_vector(cls: Type[T], target_vector: FloatArray3) -> T: + """Create a quaternion that rotates the +Z axis to point toward target_vector. + + This constructor computes the rotation that transforms the unit Z direction + [0, 0, 1] to align with the given target vector. + + Args: + target_vector: Target direction vector (will be normalized) + + Returns: + Unit quaternion representing the rotation from +Z to target_vector + + Notes: + - If target_vector points in +Z direction, returns identity quaternion + - If target_vector points in -Z direction, returns 180° rotation about +X axis + - For all other directions, uses the shortest rotation path + + Examples: + >>> # Rotate +Z to point toward +X + >>> q = Quaternion.from_unit_z_to_vector([1, 0, 0]) + >>> # Rotate +Z to point toward +Y + >>> q = Quaternion.from_unit_z_to_vector([0, 1, 0]) + """ + start = np.array([0.0, 0.0, 1.0]) + perpendicular = np.array([1.0, 0.0, 0.0]) # +X for 180° case + return cls._from_unit_direction_to_vector(start, target_vector, perpendicular) + @classmethod def from_dict(cls, obj: Any) -> "Quaternion": """Create a quaternion instance from a dictionary representation. @@ -637,3 +773,10 @@ def from_dict(cls, obj: Any) -> "Quaternion": # Note: This calls the concrete implementation's from_components # cls ensures the correct subclass type is returned return cls.from_components(w, x, y, z) + + def to_unitSphericalSmallCircle( + q: Quaternion, radius_angle: float = np.pi / 4 + ) -> UnitSphericalSmallCircle: + """ """ + a, p = q.vector_spherical + return UnitSphericalSmallCircle(azimuth=a, polar=p, radius_angle=radius_angle) diff --git a/src/pyMathTools/spherical/constructors.py b/src/pyMathTools/spherical/constructors.py index 84a59c2..2a28fe3 100644 --- a/src/pyMathTools/spherical/constructors.py +++ b/src/pyMathTools/spherical/constructors.py @@ -8,59 +8,8 @@ import numpy as np from numpy import atan2, acos, pi -from foundationTypes.mathTypes.UnitSphericalSmallCircle import UnitSphericalSmallCircle from foundationTypes.mathTypes.UnitSphericalArc import UnitSphericalArc -from pyMathTools.spatial.Quaternion import Quaternion - - -def quat_to_unitSphericalSmallCircle( - q: Quaternion, radius_angle: float = pi / 4 -) -> UnitSphericalSmallCircle: - """ - Convert a quaternion's vector part to a UnitSphericalSmallCircle. - - Uses ISO Physics Convention for spherical coordinates: - - φ (azimuth): angle in xy-plane from +x axis, range [0, 2π) - - θ (polar): angle from +z axis (colatitude), range [0, π] - - Parameters: - ----------- - q : Quaternion - Quaternion whose vector part (x, y, z) defines the center direction - radius_angle : float - Angular radius of the small circle (default: π/4) - - Returns: - -------- - UnitSphericalSmallCircle - Small circle centered at the direction defined by the quaternion's vector part - - Notes: - ------ - - If vector part is zero, defaults to north pole (azimuth=0, polar=π/2) - - Uses standard spherical-to-Cartesian formulas: - x = r*sin(θ)*cos(φ) - y = r*sin(θ)*sin(φ) - z = r*cos(θ) - """ - # Extract vector part (this represents a direction in Cartesian space) - r = np.sqrt(q.x**2 + q.y**2 + q.z**2) - - if r == 0: - # Default to equator at +x axis if vector part is zero - azimuth = 0.0 - polar = pi / 2 - else: - # Convert Cartesian (x,y,z) to spherical (azimuth, polar) - # ISO Physics Convention: - azimuth = atan2(q.y, q.x) # φ = atan2(y, x) - angle in xy-plane from +x - polar = acos(q.z / r) # θ = acos(z/r) - angle from +z axis - - return UnitSphericalSmallCircle( - azimuth=azimuth, polar=polar, radius_angle=radius_angle - ) - def arc_from_two_points( azimuth1: float, diff --git a/src/pyMathToolsPlotHelpers/plotUnitSpherical.py b/src/pyMathToolsPlotHelpers/plotUnitSpherical.py index 2782940..8ac8fad 100644 --- a/src/pyMathToolsPlotHelpers/plotUnitSpherical.py +++ b/src/pyMathToolsPlotHelpers/plotUnitSpherical.py @@ -570,7 +570,7 @@ def get_color(color_list: Optional[List[str]], idx: int) -> str: color = get_color(circle_colors, idx) # Plot the circle - ax.scatter(x, y, z, s=10, alpha=0.8, label=label, c=color) + ax.scatter(x, y, z, s=2, alpha=0.8, label=label, c=color) # Plot each arc if arcs is not None: @@ -601,11 +601,11 @@ def get_color(color_list: Optional[List[str]], idx: int) -> str: x, y, z, - linewidth=2, + linewidth=1.5, alpha=0.8, label=label, marker="o", - markersize=3, + markersize=1, color=color, ) @@ -627,16 +627,16 @@ def get_color(color_list: Optional[List[str]], idx: int) -> str: y, z, length=1.0, - arrow_length_ratio=0.15, + arrow_length_ratio=0.05, color=color, - linewidth=2.5, + linewidth=1.5, alpha=0.9, ) # Plot endpoint label = f"Quat {idx + 1}: w:{quat.w:.1f}, x:{quat.x:.1f}, y:{quat.y:.1f}, z:{quat.z:.1f}" - ax.scatter(x, y, z, s=150, marker="*", label=label, c=color, zorder=10) + ax.scatter(x, y, z, s=50, marker="*", label=label, c=color, zorder=10) ax.set_xlabel("X") ax.set_ylabel("Y")