The Algebra of Space (G3)

This notebook is an adaptation from the clifford python documentation.

Import Grassmann and instantiate a three dimensional geometric algebra

julia> using Grassmann
julia> basis"3"(⟨111⟩, v, v₁, v₂, v₃, v₁₂, v₁₃, v₂₃, v₁₂₃)

Given a three dimensional GA with the orthonormal basis $v_i\cdot v_j = \delta_{ij}$, the basis consists of scalars, three vectors, three bivectors, and a trivector.

\[\{\underbrace{v}_{\text{scalar}},\qquad\underbrace{v_1,v_2,v_3}_{\text{vectors}},\qquad\underbrace{v_{12},v_{23},v_{13}}_{\text{bivectors}},\qquad\underbrace{v_{123}}_{\text{trivector}}\}\]

The @basis macro declares the algebra and assigns the Submanifold elements to local variables. The Basis can also be assigned to G3 as

julia> G3 = Λ(3)DirectSum.Basis{⟨111⟩,8}(v, v₁, v₂, v₃, v₁₂, v₁₃, v₂₃, v₁₂₃)

You may wish to explicitly assign the blades to variables like so,

e1 = G3.v1
e2 = G3.v2
# etc ...

Or, if you're lazy you can use the macro with different local names

julia> @basis ℝ^3 E e(⟨+++⟩, v, v₁, v₂, v₃, v₁₂, v₁₃, v₂₃, v₁₂₃)
julia> e3, e123(v₃, v₁₂₃)

Basics

The basic products are available

julia> v1 * v2 # geometric productv₁₂
julia> v1 | v2 # inner product𝟎
julia> v1 ∧ v2 # exterior productv₁₂
julia> v1 ∧ v2 ∧ v3 # even more exterior productsv₁₂₃

Multivectors can be defined in terms of the basis blades. For example, you can construct a rotor as a sum of a scalar and a bivector, like so

julia> θ = π/40.7853981633974483
julia> R = cos(θ) + sin(θ)*v230.7071067811865476 + 0.7071067811865475v₂₃
julia> R = exp(θ*v23)0.7071067811865476 + 0.7071067811865475v₂₃

You can also mix grades without any reason

julia> A = 1 + 2v1 + 3v12 + 4v1231 + 2v₁ + 3v₁₂ + 4v₁₂₃

The reversion operator is accomplished with the tilde ~ in front of the Multivector on which it acts

julia> ~A1 + 2v₁ - 3v₁₂ - 4v₁₂₃

Taking a projection into a specific grade of a Multivector is usually written $\langle A\rangle_n$ and can be done using the soft brackets, like so

julia> A(0)1v
julia> A(1)2v₁ + 0v₂ + 0v₃
julia> A(2)3v₁₂ + 0v₁₃ + 0v₂₃

Using the reversion and grade projection operators, we can define the magnitude of A as $|A|^2 = \langle\tilde A A\rangle$

julia> ~A*A30 + 4v₁ + 12v₂ + 24v₃
julia> scalar(ans)30v

This is done in the abs and abs2 operators

julia> abs2(A)30 + 4v₁ + 12v₂ + 24v₃
julia> scalar(ans)30v

The dual of a multivector A can be defined as $\tilde AI$, where I is the pseudoscalar for the geometric algebra. In G3, the dual of a vector is a bivector:

julia> a = 1v1 + 2v2 + 3v31v₁ + 2v₂ + 3v₃
julia> ⋆a3v₁₂ - 2v₁₃ + 1v₂₃

Reflections

Reflecting a vector $c$ about a normalized vector $n$ is pretty simple, $c\mapsto -ncn$

julia> c = v1+v2+v3 # a vector1v₁ + 1v₂ + 1v₃
julia> n = v1 # the reflectorv₁
julia> -n*c*n # reflect a in hyperplane normal to n-1v₁ + 1v₂ + 1v₃ + 0v₁₂₃

Because we have the inv available, we can equally well reflect in un-normalized vectors using $a\mapsto n^{-1}an$

julia> a = v1+v2+v3 # the vector1v₁ + 1v₂ + 1v₃
julia> n = 3v1 # the reflector3v₁
julia> inv(n)*a*n1.0v₁ - 1.0v₂ - 1.0v₃ + 0.0v₁₂₃
julia> n\a*n1.0v₁ - 1.0v₂ - 1.0v₃ + 0.0v₁₂₃

Reflections can also be made with respect to the hyperplane normal to the vector, in which case the formula is negated.

A vector can be rotated using the formula $a\mapsto \tilde R aR$, where R is a rotor. A rotor can be defined by multiple reflections, $R = mn$ or by a plane and an angle $R = e^{\theta B/2}$. For example,

julia> R = exp(π/4*v12)0.7071067811865476 + 0.7071067811865475v₁₂
julia> ~R*v1*R2.22045e-16v₁ + 1.0v₂ + 0.0v₃ + 0.0v₁₂₃

Maybe we want to define a function which can return rotor of some angle $\theta$ in the $v_{12}$-plane, $R_{12} = e^{\theta v_{12}/2}$

R12(θ) = exp(θ/2*v12)

And use it like this

julia> R = R12(π/2)0.7071067811865476 + 0.7071067811865475v₁₂
julia> a = v1+v2+v31v₁ + 1v₂ + 1v₃
julia> ~R*a*R-1.0v₁ + 1.0v₂ + 1.0v₃ + 0.0v₁₂₃

You might as well make the angle argument a bivector, so that you can control the plane of rotation as well as the angle

R_B(B) = exp(B/2)

Then you could do

julia> Rxy = R_B(π/4*v12)0.9238795325112867 + 0.3826834323650898v₁₂
julia> Ryz = R_B(π/5*v23)0.9510565162951535 + 0.3090169943749474v₂₃

julia> R_B(π/6*(v23+v12))0.93224 + 0.255859v₁₂ + 0.0v₁₃ + 0.255859v₂₃

Maybe you want to define a function which returns a function that enacts a specified rotation, $f(B) = a\mapsto e^{B/2}\\ae^{B/2}$. This just saves you having to write out the sandwich product, which is nice if you are cascading a bunch of rotors, like so

R_factory(B) = (R = exp(B/2); a -> ~R*a*R)
Rxy = R_factory(π/3*v12)
Ryz = R_factory(π/3*v23)
Rxz = R_factory(π/3*v13)

Then you can do things like

julia> R = R_factory(π/6*(v23+v12)) # this returns a function#1 (generic function with 1 method)
julia> R(a) # which acts on a vector0.522956v₁ + 0.738144v₂ + 1.47704v₃ + 0.0v₁₂₃
julia> Rxy(Ryz(Rxz(a)))0.408494v₁ - 0.658494v₂ + 1.54904v₃ + 0.0v₁₂₃

To make cascading a sequence of rotations as concise as possible, we could define a function which takes a list of bivectors $A,B,C,...$, and enacts the sequence of rotations which they represent on some vector $x$.

julia> R_seq(args...) = (R = prod(exp.(args./2)); a -> ~R*a*R)R_seq (generic function with 1 method)
julia> R = R_seq(π/2*v23, π/2*v12, v1)#3 (generic function with 1 method)
julia> R(v1)2.22045e-16 + 3.33067e-16v₁ + 1.0v₂ + 5.55112e-17v₂₃

Barycentric Coordinates

We can find the barycentric coordinates of a point in a triangle using area ratios.

julia> function barycoords(p, a, b, c)
         ab = b-a
         ca = a-c
         bc = c-b
         A = -ab∧ca
         (bc∧(p-b)/A, ca∧(p-c)/A, ab∧(p-a)/A)
       endbarycoords (generic function with 1 method)
julia> barycoords(0.25v1+0.25v2, 0v1, 1v1, 1v2)(0.5 + 0.0v₁₂ + 0.0v₁₃ - 0.0v₂₃, 0.25 + 0.0v₁₂ + 0.0v₁₃ - 0.0v₂₃, 0.25 + 0.0v₁₂ + 0.0v₁₃ - 0.0v₂₃)

The Algebra of Space (G3)

Basics

Reflections

Rotations

Barycentric Coordinates