E=[1 0 ;0 1], I=-i*[0 1;1 0], J=[0 -1;1 0], K=[-i 0;0 i]

% Plus or minus these four matrices, together, form a group of order 8, which is a representation of the quaternion group.

% We are interested in real linear combinations of these 4 matrices.

% This collective of all of these is a 4 dimensional real vector subspace of the 8 dimensional real vector space of 2x2 complex matrices. 

% This 4 dimensional subspace is an embodiment of a structure called "the quaternions." Every mathematician and physicist used to have to know about them. The structure was invented by William Rowan Hamilton in 1843. The first edition of Maxwell's "A Treatise on Electricity And Magnetism" was written using quaternions. It is only the final (posthumous, prepared by J. J. Thomson) edition that switched to vector notation. Gibbs and Yale and vectors had, finally, defeated Hamilton and Trinity College and quaternions.

% We will call this 4 dimensional real vector space of matrices  "Quaternion Space."

% The real numbers are identified with real multiples of the identity matrix E, and linear combinations of I, J and K are called "pure" quaternions and identified with three dimensional space.

% Here is the info you need to create a multiplication table for quaternions.

I^2
J^2
K^2

I*J-K
J*K-I
K*I-J

I*J+J*I
J*K+K*J
K*I+I*K

% Any 4 dimensional space can be identified with R^4 and v2q and q2v below implement that identification. They are inverse functions to each other: isomorphisms between quaternion space and R^4. 

% Just as complex numbers can be "thought of" as points in R^2, so too quaternions can be "thought of" as points in R^4. But in that representation the quaternion multiplication is similar to but more complicated than complex multiplication of points in the complex plane. 

% In quaternion space quaternion multiplication is just matrix multiplication, quite efficient and simple to carry out.

v2q=@(v) v(1)*E+v(2)*I+v(3)*J+v(4)*K

v=[3;5;7;11]

W=v2q(v)

q2v=@(W) [ real(W(1,1)); -imag(W(2,1)); real(W(2,1)); -imag(W(1,1))]

q2v(W)

% When you represent a quaternion in R^4, the pure quaternion (or imaginary) part "is" a vector in R^3:

qim=@(v) v(2:4)

qim(v)

% When you represent a quaternion in R^4, the real part "is" a real number:

qre=@(v) v(1)

qre(v)

% We emphasize that real 3-d vectors (i.e. vectors in "space") are identified with the 3-d subspace of pure quaternions. 

% Cross product and dot product can be applied to the space part of quaternions when they are represented in R^4.

qcross=@(v,w) cross(v(2:4),w(2:4))

qdot=@(v,w) v(2:4)'*w(2:4)

% Consider the quaternions X, Y

X=3*E-7*I+3*J-5*K

Y=1*E+2*I+3*J+4*K

% and associated 4-vectors A,B:

A=q2v(X)
B=q2v(Y)

% The quaternion product of these is 

X*Y

% with vector representation 

C=q2v(X*Y)

% It is a little hard to see the relationship between the vector representations of X and Y and their quaternion product C. However, there IS a formula involving cross product and dot product in combinations. Here it is.

[A(1)*B(1) - qdot(A,B) ; A(1)*qim(B) + B(1)*qim(A) + qcross(A,B)]

% If you examine it carefully, you can see that it is only the cross product term in the pure quaternion part of the product formula that causes quaternion product to be non-commutative. It is also easy to see why it took Hamilton so long to invent the quaternion product. He was NOT working with a matrix representation of the quaternions (and linear algebra did not exist in modern form then) just with points in R^4. Of course, the real inspiration was realizing it might be useful to use R^4 in the first place.

norm(A)
norm(X)

inv(X)

X*(X'/norm(X)^2)

(X'/norm(X)^2)*X

% The quaternion inverse is matrix inverse in quaternion space. Every quaternion except the 0 quaternion has a multiplicative inverse, just like in the real numbers and the complex numbers. 

% Quaternion multiplication is matrix multiplication of 2x2 matrices.

X*inv(X)


% If v is a real 3-d vector p=Versor(theta,v), defined below, is a quaternion that will rotate a real 3-d vector w by angle theta around axis v when calculated as the quaternion product q=p*v2q([0;w])*p'. 

% This quaternion will have real part zero and the rotated vector in R^3 is then qim(q2v(q)).

% For a proof of that, see my extended notes. But if you believe it (trust me) it is quite straightforward to actually use.

% The following is the vector form of the versor: not that useful.

versor= @(theta,v) [cos(theta/2);sin(theta/2)*v/norm(v)]

% Here is the quaternion form of the versor: very useful.

Versor= @(theta,v) v2q([cos(theta/2);sin(theta/2)*v/norm(v)])


% We are creating the rotated arrows (around axis [1;0.2;0.2]) below:

for k=1:21
theta=(k-1)*pi/20;
p=Versor(theta,[1;0.2;0.2]);
z(:,k)=q2v(p*v2q([0;1;0;0])*(p'));
end


figure(40)
quiver3(0, 0, 0, 1, 0.2, 0.2)
% The following line will guarantee the rotation is not "distorted" 
% by different units on different axes.
set(gca,'DataAspectRatio',[1 1 1])
hold on
for k=1:21
quiver3(0,0,0, z(2,k), z(3,k), z(4,k))
end
hold off


% Now we'll create a 3-d object (the convex hull of some points, a polyhedron) and rotate that by 2*pi/3 around axis [10;2; 2].

W=[1 3 5 ; 2 5 8 ; 9 10 11; 12 8 10 ]

p=Versor(2*pi/3,[10;2; 2]);

clear z
for k=1:4
z(:,k)=q2v(p*v2q([0;W(k,:)'])*(p'));
end
Q=(z(2:4,:))'

% Look down the axis of rotation (so the tip of the arrow is pointing at you) to see the rotation of the polyhedron most clearly.


DT = delaunayTriangulation(W);
[C,v] = convexHull(DT);
trisurf(C,DT.Points(:,1),DT.Points(:,2),DT.Points(:,3),'FaceColor','cyan')
hold on

DT2 = delaunayTriangulation(Q);
[C2,v2] = convexHull(DT2);
trisurf(C2,DT2.Points(:,1),DT2.Points(:,2),DT2.Points(:,3),'FaceColor','red')

quiver3(0,0,0, 10, 2, 2)

set(gca,'DataAspectRatio',[1 1 1])

hold off


% As an interesting aside, there is an important basis of the 8 dimensional real vector space of the 2x2 complex matrices. It is generated by

s0=E
s1=i*I
s2=i*J
s3=i*K

% These matrices are hermitian, which means they are their own conjugate transpose, and so correspond to ``observables'' in quantum mechanics. 

% They represent spin with respect to the coordinate axes in quantum mechanical descriptions of certain particles. 

% Any 2x2 hermitian matrix can be found (in one way) as a real combination of these Pauli spin matrices. The intersection of this 4 dimensional real vector apace with Quaternion Space is just the one dimensional multiples of the identity E=s0.

% These spin matrices have a different multiplication table from the quaternions: note the first and second groups below.

s1*s1
s2*s2
s3*s3

s1*s2-i*s3
s2*s3-i*s1
s3*s1-i*s2

s1*s2+s2*s1
s2*s3+s3*s2
s3*s1+s1*s3

% The real vector space of hermitian matrices do not form a real algebra: they are not closed under multiplication, as can be seen from the second grouping of relations above.

% (+/-)1 and (+/-)i times these four matrices, together, form a group of order 16, which is a representation of the Pauli spin group.

% And the 8 with positive choice form a basis of the (8 dimensional) real vector space of 2x2 matrices with complex entries.