sympy-ga-cl30+.py

   # Description: Tutorial code for Cl(3, 0)+

if True:

### Initialization

    # Import all symbols from the `sympy` and `sympy.galgebra` packages.
    from sympy.galgebra import *
    from sympy import *

### Defining a basis

    # Theorem:
    # Cl(3, 0)+ is algebra-isomorphic to quaternions

    # Set up Cl_{3, 0}.
    e1, e2, e3 = MV.setup('e1 e2 e3', metric='[1, 1, 1]')
    i = e1 * e2
    j = e1 * e3
    k = e2 * e3

    print
    print 'Algebra of quaternions Cl(3, 0)+'
    print '--------------------------------'
    print

### Quaternion multiplication

    # Set up an arbitrary quaternion.
    a0, a1, a2, a3 = symbols('a0 a1 a2 a3')
    a = a0 + a1 * i + a2 * j + a3 * k

    # Set up another arbitrary quaternion.
    b0, b1, b2, b3 = symbols('b0 b1 b2 b3')
    b = b0 + b1 * i + b2 * j + b3 * k

    # What is the product of `a` and `b`?
    c = a * b

    # Print one component per line.
    print 'Quaternion multiplication:'
    print c.Fmt(3)
    print

    # I don't know how to make the package to print 
    # `i`, `j`, and `k`.

### Quaternion rotation

    # Define a rotation axis.
    v = 9 * k

    # Define a rotation angle.
    alpha = pi / 2

    # Construct a quaternion representing
    # rotation around `v` with angle `alpha`.
    q = ((alpha / 2) * (v / v.norm())).exp()

    # This is the same as
    # q = cos(alpha / 2) + sin(alpha / 2) * (v / v.norm())

    # Define a vector to rotate.
    v = 3 * i + 6 * j - 7 * k

    # Rotate the vector.
    rv = q * v * inv(q)

    # Output the rotated vector.
    print 'Rotate', v, 'around k', '=', rv

### Quaternion interpolation

    # Define two rotation axes.
    v1 = 5 * i - 3 * j + 9 * k
    v2 = 3 * i + 4 * j + k

    # Define two rotation angles.
    alpha1 = pi / 8
    alpha2 = pi / 2

    # Compute the rotation quaternions.
    q1 = ((alpha1 / 2) * (v1 / v1.norm())).exp()
    q2 = ((alpha2 / 2) * (v2 / v2.norm())).exp()

    # Interpolate midpoint between the quaternions.
    t = 0.5

    # This is what should be done:
    # q = exp((1 - t) * log(q1) + t * log(q2))
    # Unfortunately, the package does not support a logarithm.

    # Alternatively, one could do
    # q = q1 * (inv(q1) * q2)**t * q2
    # Unfortunately, the package does not support a power either.

    # To create an interesting animation, you would
    # linearly interpolate t from 0 to 1, and then
    # observe how the object rotates from one
    # position to another.

### Quaternion reflection

    # Define a vector along which to reflect.
    q = k

    # Define a vector to reflect.
    v = 3 * i + 6 * j - 7 * k

    # Reflect along q.
    rv = -q * v * inv(q)

    # Note that this is a reflection along `q`.
    # We can also define reflection through `q`;
    # just remove the minus. 

    # Output the reflected vector.
    print 'Reflect', v, 'along k', '=', rv