# Representing rotations by quaternions

The most usual ways to represent a rotation in three-dimensional space are rotations matrices or Euler angles. In this package, nevertheless, they are represented using unit quaternions.

Quaternions are a four-dimensional non-commutative algebra with the property that any rotation may be represented by a quaternion of norm 1. Some of their advantages consist in that we only need 4 parameters to represent a quaternion (as opposed to the 9 needed by rotation matrices) and that they do not suffer from the phenomenon of gimbal lock (as opposed to Euler angles).

If a quaternion $q$ represents a rotation matrix $R$, their action on a vector is defined as

This package assumes that your coordinate system follows the right-hand rule. If, for whatever reason, you want to use a left-handed system, $q$ and $q^{-1}$ must be transposed on the above formula.

## The `rotate`

function

All rotations on `DeformableBodies.jl`

are done by the function `rotate`

. Its signature consists of

`rotate(v::Vector, q::Quaternion, center::Vector)`

This function rotates the vector `v`

by the rotation represented by `q`

in the frame of reference centered on `center`

.

If the third argument is left blank, it defaults to the origin.

```
julia> v = [3,0,0]
3-element Array{Int64,1}:
3
0
0
julia> q = Quaternion(1,2,3,4)
1 + 2i + 3j + 4k
julia> rotate(v, q, [0,0,0])
3-element Array{Float64,1}:
-2.0
2.0
1.0
julia> rotate(v, q)
3-element Array{Float64,1}:
-2.0
2.0
1.0
```

## Representing the identity rotation

The identity rotation is the matrix $I$ satisfying $I v = v$ for all $v$. This is represented by the quaternion $1$, which can be gotten using Julia's multiple dispatch via the expression `one(Quaternion)`

.

```
julia> q = one(Quaternion)
1 + 0i + 0j + 0k
julia> rotate([1,2,3], q)
3-element Array{Float64,1}:
1.0
2.0
3.0
```

## Composition of rotations

The composition of rotations translates to the quaternion world as the product of quaternions. Thus, it is the same to rotate a vector by $q_1$ and then by $q_2$ and to rotate by $q_2\,q_1$.

```
julia> v = [9, 0, 0]
3-element Array{Int64,1}:
9
0
0
julia> q1 = Quaternion(1,2,3,4)
1 + 2i + 3j + 4k
julia> q2 = Quaternion(4,3,2,1)
4 + 3i + 2j + 1k
julia> rotate(v, q2*q1)
3-element Array{Float64,1}:
-1.0
-4.0
8.0
julia> rotate(rotate(v, q1), q2)
3-element Array{Float64,1}:
-1.0
-4.0
8.0
```

Rotations are non-commutative. Applying $q_2$ after $q_1$ is the same as multiplying $q_2$ to the **left** of $q_1$.

## Axis-angle representation

An important property of three-dimensional space it that every rotation fixes a line. Therefore, they accept an axis-angle representation. That is, a rotation $R$ is defined by a unit vector $\hat{n}$ and an angle $\theta$ such that $R$ rotates a vector counterclockwise by an amount of $\theta$ around the line defined by $\hat{n}$.

To get the quaternion representing an rotation of `θ`

around `n`

, use the method `axistoquaternion`

. You may pass an unormalized vector and the method will normalize it.

```
julia> q = axistoquaternion(axis = [0,0,2], angle = π/2)
0.7071067811865476 + 0.0i + 0.0j + 0.7071067811865475k
julia> rotate([1,0,0], q)
3-element Array{Float64,1}:
2.220446049250313e-16
1.0
0.0
```

In fact, this combination is so useful that the function `rotate`

is overloaded to directly convert from an axis-angle pair.

```
julia> rotate([1,0,0], axis=[0,0,1], angle=π/2)
3-element Array{Float64,1}:
2.220446049250313e-16
1.0
0.0
```

This version of `rotate`

also accepts the optional argument `center`

, which defaults to the origin.

`rotate([1,0,0], axis=[0,0,1], angle=π/2, center=[0,1,0])`

## Converting between matrices and quaternions

Sometimes a rotation may already come represented as a matrix or the transition rotations from `solve!`

may be needed in matrix form some application. To deal with these cases, the package provides two auxiliary functions `matrixtoquaternion`

and `quaterniontomatrix`

.

To convert a rotation matrix to a unit quaternion, do

```
julia> R = [cos(π/4) -sin(π/4) 0;
sin(π/4) cos(π/4) 0;
0 0 1]
3×3 Array{Float64,2}:
0.707107 -0.707107 0.0
0.707107 0.707107 0.0
0.0 0.0 1.0
julia> R * [sqrt(2), 0, 0]
3-element Array{Float64,1}:
1.0000000000000002
1.0
0.0
julia> q = matrixtoquaternion(R)
0.9238795325112867 + 0.0i + 0.0j + 0.3826834323650898k
julia> rotate([sqrt(2), 0, 0], q)
3-element Array{Float64,1}:
0.9999999999999999
1.0
0.0
```

Some minor differences may occur due to floating-point rounding errors.

The function `matrixtoquaternion`

assumes that the input is a *rotation matrix* but, for efficency reasons, no check is done in this regard. If you do not make sure beforehand that the matrix is orthogonal, bad things may happen.

To convert a quaternion to a matrix simply do

```
julia> q = Quaternion(1,2,3,4)
1 + 2i + 3j + 4k
julia> R = quaterniontomatrix(q)
3×3 Array{Float64,2}:
-0.666667 0.133333 0.733333
0.666667 -0.333333 0.666667
0.333333 0.933333 0.133333
julia> rotate([3,0,0], q)
3-element Array{Float64,1}:
-2.0
2.0
1.0
julia> R * [3,0,0]
3-element Array{Float64,1}:
-2.0
2.0
1.0
```

The function `quaterniontomatrix`

works for every quaternion, and does not require the input to be a unit quaternion.

There is a **unique** rotation matrix representing a given quaternion but there are **two** unit quaternions representing the same matrix.

This means that `quaterniontomatrix ∘ matrixtoquaternion`

equals the identity (disconsidering floating-point rouding errors) but the opposite is not in general true. For a simple example.

```
q = Quaternion(-1.0)
(matrixtoquaternion ∘ quaterniontomatrix)(q)
```

Nevertheless, both these quaternions produce the same rotations.

## Rotating a PointMass

All the previous functionalities only require the submodule `Quaternions`

and work directly with vectors. Nevertheless, the models on `DeformableBodies.jl`

are constructed with respect to the type `PointMass`

. To help with that, the function rotate is overloaded to directly rotate the position of a `PointMass`

without interfering with its mass.

```
rotate(p::PointMass, q::Quaternion, center::Vector)
rotate(p::PointMass; axis::Vector, angle::Real, center::Vector)
```

The usage is identical to the version for vectors including the fact that the argument `center`

must be a Vector and not another PointMass.

An usual application consists in rotating a body around its center of mass.

```
julia> body = [ PointMass(rand(), rand(3)) for i in 1:5 ]
5-element Array{PointMass{Float64},1}:
PointMass{Float64}(0.5887392979926158, [0.9851461509508204, 0.5871472754686167, 0.18988173568255662])
PointMass{Float64}(0.16447626657374115, [0.15917411711754115, 0.21781301182062407, 0.6154574686908942])
PointMass{Float64}(0.6294684337802217, [0.6007476400909275, 0.8297934043579711, 0.011699945759655606])
PointMass{Float64}(0.3674998635545208, [0.24634732086407585, 0.2279264662827194, 0.7397194326716405])
PointMass{Float64}(0.14298635160866868, [0.6206875249356907, 0.26670397043925376, 0.6435613516386134])
julia> center_of_mass(body)
3-element Array{Float64,1}:
0.6146350376396904
0.54180467556897
0.30860986988838707
julia> q = Quaternion(1,2,3,4)
1 + 2i + 3j + 4k
julia> rotated = [ rotate(p, q, center_of_mass(body)) for p in body ]
5-element Array{PointMass{Float64},1}:
PointMass{Float64}(0.5887392979926158, [0.2866060103346143, 0.6945457950059541, 0.4586029163376565])
PointMass{Float64}(0.16447626657374115, [1.1000983352765157, 0.550726349005324, -0.104689643277118])
PointMass{Float64}(0.6294684337802217, [0.4445578554829961, 0.23861021818764017, 0.5331822276913696])
PointMass{Float64}(0.3674998635545208, [1.1344567669596526, 0.6883119760028129, -0.049624422666217916])
PointMass{Float64}(0.14298635160866868, [0.8195510387085607, 0.8608408899763599, 0.0985269050993488])
julia> center_of_mass(rotated)
3-element Array{Float64,1}:
0.6146350376396904
0.5418046755689702
0.308609869888387
```