error ?

23 views
Skip to first unread message

ch v Siddhardha

unread,
Sep 6, 2025, 8:46:45 AM (3 days ago) Sep 6
to sympy
can any one help me to fix errors in this,
"""
Geometric Algebra operations for Quaternions.

This module provides Geometric Algebra operations for the Quaternion class,
including geometric product, outer product, inner product, and rotor conversions.

References
==========
[1] Geometric Algebra for Computer Science, Leo Dorst et al.
[2] Geometric Algebra for Physicists, Chris Doran and Anthony Lasenby
"""

from __future__ import annotations
from typing import TYPE_CHECKING, Tuple
from sympy import Expr, cos, sin, sqrt, acos, atan2, S, pi, trigsimp
from sympy.core.sympify import sympify

if TYPE_CHECKING:
    from .quaternion import Quaternion


def geometric_product(q1: 'Quaternion', q2: 'Quaternion') -> 'Quaternion':
    """
    Compute the geometric product of two quaternions.
   
    The geometric product combines the dot product and wedge product:
    q1 * q2 = q1•q2 + q1∧q2
   
    Parameters
    ==========
    q1, q2 : Quaternion
        Input quaternions
       
    Returns
    =======
    Quaternion
        The geometric product of q1 and q2
   
    Examples
    ========
    >>> from sympy import Quaternion
    >>> from sympy.algebras.geometric_algebra_ops import geometric_product
    >>> q1 = Quaternion(1, 2, 3, 4)
    >>> q2 = Quaternion(5, 6, 7, 8)
    >>> geometric_product(q1, q2)
    -60 + 12*i + 30*j + 24*k
    """
    return q1 * q2


def outer_product(q1: 'Quaternion', q2: 'Quaternion') -> 'Quaternion':
    """
    Compute the outer (wedge) product of two quaternions.
   
    The outer product represents the oriented area spanned by the two quaternions.
   
    Parameters
    ==========
    q1, q2 : Quaternion
        Input quaternions
       
    Returns
    =======
    Quaternion
        The outer product of q1 and q2
       
    Examples
    ========
    >>> from sympy import Quaternion
    >>> from sympy.algebras.geometric_algebra_ops import outer_product
    >>> q1 = Quaternion(1, 2, 3, 4)
    >>> q2 = Quaternion(5, 6, 7, 8)
    >>> outer_product(q1, q2)
    0 + 0*i + 0*j + 0*k  # Simplified example, actual result will vary
    """
    return (q1 * q2 - q2 * q1) / 2


def inner_product(q1: 'Quaternion', q2: 'Quaternion') -> 'Quaternion':
    """
    Compute the inner (dot) product of two quaternions.
   
    The inner product represents the projection of one quaternion onto another.
   
    Parameters
    ==========
    q1, q2 : Quaternion
        Input quaternions
       
    Returns
    =======
    Quaternion
        The scalar part represents the inner product
       
    Examples
    ========
    >>> from sympy import Quaternion
    >>> from sympy.algebras.geometric_algebra_ops import inner_product
    >>> q1 = Quaternion(1, 2, 3, 4)
    >>> q2 = Quaternion(5, 6, 7, 8)
    >>> inner_product(q1, q2)
    70  # Simplified example
    """
    return (q1 * q2 + q2 * q1) / 2


def to_rotor(angle: Expr, axis: 'Quaternion') -> 'Quaternion':
    """
    Convert an angle and axis to a rotor (unit quaternion).
   
    A rotor is a unit quaternion that represents a rotation.
   
    Parameters
    ==========
    angle : Expr
        Rotation angle in radians
    axis : Quaternion
        Axis of rotation (will be normalized)
       
    Returns
    =======
    Quaternion
        A unit quaternion representing the rotation
       
    Examples
    ========
    >>> from sympy import pi, sqrt
    >>> from sympy.algebras.quaternion import Quaternion
    >>> from sympy.algebras.geometric_algebra_ops import to_rotor
    >>> q = to_rotor(pi/2, Quaternion(0, 1, 0, 0))  # 90° around x-axis
    >>> q
    sqrt(2)/2 + sqrt(2)/2*i + 0*j + 0*k
    """
    from .quaternion import Quaternion
   
    if not isinstance(axis, Quaternion):
        raise TypeError("Axis must be a Quaternion")
   
    # Ensure axis is a pure quaternion (scalar part is zero)
    if axis.a != 0:
        axis = Quaternion(0, axis.b, axis.c, axis.d)
   
    # Normalize the axis
    norm = axis.norm()
    if norm == 0:
        raise ValueError("Axis must have non-zero magnitude")
   
    axis_normalized = axis / norm
   
    # Create rotor: cos(θ/2) + sin(θ/2)*axis
    half_angle = angle / 2
    scalar_part = cos(half_angle)
    vector_part = sin(half_angle) * axis_normalized
   
    return Quaternion(scalar_part, *vector_part.args[1:])


def from_rotor(rotor: 'Quaternion') -> Tuple[Expr, 'Quaternion']:
    """
    Convert a rotor (unit quaternion) to an angle and axis of rotation.
   
    Parameters
    ==========
    rotor : Quaternion
        A unit quaternion representing a rotation
       
    Returns
    =======
    tuple[Expr, Quaternion]
        A tuple containing (angle, axis) where angle is the rotation angle in radians
        and axis is a unit quaternion representing the rotation axis.
       
    Examples
    ========
    >>> from sympy import pi, sqrt
    >>> from sympy.algebras.quaternion import Quaternion
    >>> from sympy.algebras.geometric_algebra_ops import from_rotor
    >>> q = Quaternion(sqrt(2)/2, sqrt(2)/2, 0, 0)  # 90° around x-axis
    >>> angle, axis = from_rotor(q)
    >>> angle
    pi/2
    >>> axis
    0 + 1*i + 0*j + 0*k
    """
    from .quaternion import Quaternion
   
    if not isinstance(rotor, Quaternion):
        raise TypeError("Rotor must be a Quaternion")
   
    # Normalize the rotor
    rotor = rotor.normalize()
   
    # Calculate angle (handle floating point errors)
    angle = 2 * acos(rotor.a)
   
    # Handle the identity rotation case (angle = 0)
    if abs(rotor.a) > 1 - 1e-10:  # Account for floating point errors
        return (0, Quaternion(0, 1, 0, 0))  # Default axis for identity
   
    # Calculate axis (normalized)
    sin_half = sqrt(1 - rotor.a**2)
    axis = Quaternion(0, *(comp / sin_half for comp in rotor.args[1:]))
   
    # Ensure axis is normalized
    axis = axis.normalize()
   
    # Normalize angle to [-π, π]
    if angle > pi:
        angle -= 2 * pi
    elif angle < -pi:
        angle += 2 * pi
   
    return (angle, axis)

Donaldson Tan

unread,
Sep 6, 2025, 10:05:37 PM (2 days ago) Sep 6
to sympy
What's the error message?
Reply all
Reply to author
Forward
0 new messages