introduction of rotation with quaternions (part 2),version  1.0
peter.vlasschaert@gmail.com, 02/11/2018

load packages.
=============

kill(all)$
load(vect)$
load(lrats)$

declare([a1,a2,a3,a4,b1,b2,b3,b4,c1,c2,c3,c4],scalar);

  definition : multiplication, coef  ' quaternions '
========================================

multiplication : =>   h = qa.qb  , {qa,qb,h } : quaternion
===============

qmult2(a1,a2,a3,a4,b1,b2,b3,b4,h):=(qa:a1+a2*i+a3*j+a4*k,qb:b1+b2*i+b3*j+b4*k,h1:expand(qa.qb),
                                                                h2:[i.j=k,j.i=-k,j.k=i,k.j=-i,k.i=j,i.k=-j,i.i=-1,j.j=-1,k.k=-1],h:subst(h2,h1))$

rem : i.j=k,j.i=-k,j.k=i,k.j=-i,k.i=j,i.k=-j  , cyclic permutation :  ' i => j => k => i ',same rules cross product.

coef                 : =>    given :  q = a1+a2*i+a3*j+a4*k  => [a1,a2,a3,a4]
===============

coefq(qq,ss):=(s1:subst([i=0,j=0,k=0],qq),s2:coeff(qq,i),s3:coeff(qq,j),s4:coeff(qq,k),ss:[s1,s2,s3,s4])$

new : qmult2,coefq  
rem : . = non commutative multiplication.

How to convert quaternion : with matrix respresentation of the quaternion.
==============================================================

see : complex number (q, j=k=0) , introduction of rotation with quaternions (part 1),version  1.0
---------------------------------------------------------------------------------------------------------------------------------------------
def : z = complex number = q = a1+a2*i  ,   i^2 =-1
----------------------------------------------------------------------------------------------------------------------------------------------

qa = a1+a2*i+a3*j+a4*k ;  ' first quaternion        '
qb = b1+b2*i+b3*j+b4*k;  ' second quaternion '

p1:qmult2(a1,a2,a3,a4,b1,b2,b3,b4,h);

result :  h  ?
==========

p2:coefq(p1,ss);

coef b1:
=======

p3a:coeff(p2,b1);

p3a1:transpose(p3a);

coef b2:
=======

p3b:coeff(p2,b2);

p3b1:transpose(p3b);

coef b3:
=======

p3c:coeff(p2,b3);

p3c1:transpose(p3c);

coef b4:
=======

p3d:coeff(p2,b4);

p3d1:transpose(p3d);

matrix representation ,row form.
===========================

p4:matrix(p3a,p3b,p3c,p3d);

matrix representation  , ' qa'
=======================
qa = a1+a2*i+a3*j+a4*k ;  ' first quaternion        '

p4a:transpose(p4);

p4a[1,2];

vector representation , 'qb'
=======================

p4b:matrix([b1],[b2],[b3],[b4]);

check :  see relation between 'p5' and 'p2'
                ============================

'-------------------------------------------------------------------------------------------------------------------------------------'

p5:p4a.p4b;

p6a:p5[1]=[p2[1]];

is(p6a);

p6b:is(p5[2]=[p2[2]]);
p6c:is(p5[3]=[p2[3]]);
p6d:is(p5[4]=[p2[4]]);

'-------------------------------------------------------------------------------------------------------------------------------------'

rem : is(... = ....), use for check 

What are matrix representation  1, i , j ,k
=================================
see :  'p4a' , matrix representation  'qa'
qa = a1+ a2*i+a3*j+a4*k
------------------------------------------------------------

a)  1   =>  a1=1,a2=a3=a4=0
======================

p7a:ev(p4a,a1=1,a2=0,a3=0,a4=0);

rem : identity matrix (4);

b)  i   =>  a2 =1  ,a1=a3=a4=0
======================

p7b:ev(p4a,a2=1,a1=0,a3=0,a4=0);

c)  j   =>  a3 =1  ,a1=a2=a4=0
=======================

p7c:ev(p4a,a3=1,a2=0,a1=0,a4=0);

d)  k   => a4 =1  ,a2=a3=a1=0
=======================

p7d:ev(p4a,a4=1,a2=0,a3=0,a1=0);

representation of vector  : qa(a1,[a2,a3,a4]) ,a1=0
=========================================
vector = [a2,a3,a4]

p7e:ev(p4a,a1=0);

rem :main diagonal ,see ones , 'p7a',reflection around main diagonal , ai => -ai  'or' -ai => ai
i = 1,2,3,4
rem :  %i , 'complex number Maxima ' =>  i ,for general quaternion representation.
( introduction of rotation with quaternions (part 1),version  1.0)

cqa :   conjugate : qa ,  ' Matrix  representation '
========================================================
complex number :  z = a + b*i  , complex conjugate : z=a-b*i
Ned : " complex getal : z =a+b*i , complex toegevoegde : z=a-b*i  "
========================================================
cqa = conjugate quaternion  'qa'
-----------------------------------------------------------------------------------------------------
use : 'p4a'

cqa :  qa =>  cqa = a1-a2*i-a3*j-a4*k   ' Matrix representation ';
===================================================

p8:ev(p4a,a2=-a2,a3=-a3,a4=-a4);

rem : definition pure quaternion , qa =a1+a2*i+a3*j+a4*k ,when a1=0
         =====================================================
result : vector =[a2,a3,a4].
rem : vector will be rotated ,  q = x*i+y*j+z*k  , {x,y,z are cartesian coordinates}
result : vector =[x,y,z]
---------------------------------------------------------------------------------------------------------------------
total result : vector = pure quaternion. ( no real part ).
---------------------------------------------------------------------------------------------------------------------

rotation : quaternion  from one coordinate to other coordinate 
====================================================
----------------------------------------------------------------------------------------------
result : qar, ' (q)uaternion (a)fter (r)otation '.
----------------------------------------------------------------------------------------------

definition rotation , ' qa' 
===================

qar = q1a.(x*i+y*j+z*k).cq1a

--------------------------------------------------------------------------
qa     = a1+a2*i+a3*j+a4*k   =>  (|q1a|    = 1)
cqa   = a1 -a2*i-a3*j-a4*k     =>  (|cq1a|  = 1)
vector =[x,y,z] , pure quaternion
--------------------------------------------------------------------------
result  : 'qar'  will be pure quaternion
rem     :   | |      = norm. ( two norm = distance norm)
--------------------------------------------------------------------------

a) norm of quaternion
===================

declare([q1,q2,q3,q4],scalar);

definition : norm of quaternion  'normq'
=================================

normq(a1,a2,a3,a4,tt):=( tt:sqrt((a1)^2+(a2)^2+(a3)^2+(a4)^2))$

p9:normq(a1,a2,a3,a4,tt);

definition : convert quaternion, 'qa' to unit quaternion 'q1a'
=================================================

convertunitq(a1,a2,a3,a4,ll):=(tt:sqrt((a1)^2+(a2)^2+(a3)^2+(a4)^2),ll:[a1/tt,a2/tt,a3/tt,a4/tt])$

p10:convertunitq(a1,a2,a3,a4,ll);

p10a:p10[1]=q1;
p10b:p10[2]=q2;
p10c:p10[3] =q3;
p10d:p10[4] =q4;

b) make  unit quaternions
======================

quaternion :
-----------------

p11:q1a=[rhs(p10a),rhs(p10b),rhs(p10c),rhs(p10d)];

conjugate quaternion :
----------------------------------

p11a:cq1a=[rhs(p10a),-rhs(p10b),-rhs(p10c),-rhs(p10d)];

p11a1:rhs(p11a)[1];
p11a2:rhs(p11a)[2];
p11a3:rhs(p11a)[3];
p11a4:rhs(p11a)[4];

vector : ' pure quaternion'
-----------

declare([x,y,z],scalar);

p11b:vec=[0,x,y,z];

c) Formula   ' qar', see above.
==========

first part : qh1 = q1a.vec
------------------------------------

p12a:qmult2(p11a1,-p11a2,-p11a3,-p11a4,0,x,y,z,qh1);

p12b:coefq(p12a,ss1);

p12b1:p12b[1];
p12b2:p12b[2];
p12b3:p12b[3];
p12b4:p12b[4];

second part : qh2 = qh1.cq1a
--------------------------------------------

p13a:qmult2(p12b1,p12b2,p12b3,p12b4,p11a1,p11a2,p11a3,p11a4,qh2);

d) coef from 'qar'
==============

coef: i,j,k
========

p13a1:coeff(p13a,i);
p13a2:coeff(p13a,j);
p13a3:coeff(p13a,k);

coef : x,y,z  from p13a1 and p13a2 and p13a3
======================================

p13a1 :?  x,y,z
---------

p13a1x:coeff(p13a1,x);
p13a1y:coeff(p13a1,y);
p13a1z:coeff(p13a1,z);

p13a2 :? x,y,z
---------

p13a2x:coeff(p13a2,x);
p13a2y:coeff(p13a2,y);
p13a2z:coeff(p13a2,z);

p13a3 :? x,y,z
---------

p13a3x:coeff(p13a3,x);
p13a3y:coeff(p13a3,y);
p13a3z:coeff(p13a3,z);

e)  matrix ' representation' , qar
==========================

rem : |q| =q1^2+q2^2+q3^2+q4^2 = 1

check : true 

p14:p10a^2+p10b^2+p10c^2+p10d^2;

p14a:ratsimp(p14);

p14a1:rhs(p14a)=lhs(p14a);

fill in matrix  'M' ,menu above.
How to :  Menu: Algebra : Enter Matrix  : Rows => 3
                                                                             Colummns => 3
                                                                             general
                                                                             M

M: matrix(
 [p13a1x,p13a1y,p13a1z], 
 [p13a2x,p13a2y,p13a2z], 
 [p13a3x,p13a3y,p13a3z]
);

subst : use  ' p14a1'

p15a:[-q4^2-q3^2=q2^2+q1^2-1];
p15b:[-q4^2-q2^2=q1^2+q3^2-1];
p15c:[-q3^2-q2^2=q4^2+q1^2-1];

subst , main diagonal of  M

M[1,1];

p16a:lratsubst(p15a,M[1,1]);
p16b:lratsubst(p15b,M[2,2]);
p16c:lratsubst(p15c,M[3,3]);

M[1,1]:p16a;
M[2,2]:p16b;
M[3,3]:p16c;

final matrix : MM = M
=================

M;

p17:fullratsimp(M/2);

MM:2*p17;

to find trace matrix : sum elements of main diagonal  or load packages.

a) load packages 
==============

load(functs)$

p18:tracematrix(M);

b)  by sum of the elements on main diagonal.
======================================

p19:sum(M[i,i],i,1,length(M));

check
---------

is(p18=p19);

TM = trace matrix (M) , by definition.
==============================

TM:p19;

'--------------------------------------------------------------------------------------------------------------------------------------'