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

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

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

see:introduction of rotation with quaternions (part 3),version  1.0 
use M :matrix. relation rotation matrix for  Mα,Mβ,Mγ

Euler angles  => quaternion
=======================

a)  Mα matrix : rotation ,about x axes with angle = α 
------------------------------------------------------------------------------

p1a:[1,0,0];
p1b:[0,cos(α),-sin(α)];
p1c:[0,sin(α),cos(α)];

Mα:matrix(p1a,p1b,p1c);

TMα:tracematrix(Mα);

b) find q1α,q2α,q3α,q4α
-----------------------------------

introduction of rotation with quaternions (part 3),version  1.0

M:matrix(
        [2*q2^2+2*q1^2-1,   2*q2*q3-2*q1*q4,    2*q2*q4+2*q1*q3],
        [2*q1*q4+2*q2*q3,   2*q3^2+2*q1^2-1,    2*q3*q4-2*q1*q2],
        [2*q2*q4-2*q1*q3,   2*q3*q4+2*q1*q2,    2*q4^2+2*q1^2-1]
    );

[q1,q2,q3,q4]  => [q1α,q2α,q3α,q4α]
TM                     => TMα
M[i,i]                  =>  Mα[i,i]  ,  i=1,2,3    elements of the main diagonal

q1α:sqrt(TMα+1)/2,radcan;
q2α:sqrt(-TMα+2*Mα[1,1]+1)/2;
q3α:sqrt(-TMα+2*Mα[2,2]+1)/2;
q4α:sqrt(-TMα+2*Mα[3,3]+1)/2;

simplify : [q1α,q2α,q3α,q4α]
------------------------------------------

double angles
============

p2:'cos(2*α) = trigexpand(cos(2*α)) ;

matchdeclare([s],true);
let(sin(s)^2,1-cos(s)^2);
p3:letsimp(p2);

p3a:subst([α=α/2],p3);

p3b:solve(p3a,cos(α/2));

p3c:p3b[2];

p3d:subst([cos(α)^2=1-sin(α)^2],p2);

p3e:subst([α=α/2],p3d);

p3f:solve(p3e,sin(α/2));

p3g:p3f[2];

rem:  sqrt(1-cos(α))/sqrt(2) =sqrt(2-2*cos(α))/2.

result : quaternion , qα
===================

qα:[q1α,q2α,q3α,q4α];
qα:subst([q1α =lhs(p3c),q2α=lhs(p3g)],qα);

p4:'qα =qα[1]+qα[2]*i+qα[3]*j+qα[4]*k;

analog : Mα,Mβ,Mγ,see website rotation
=================================

rotation : about x-axis   ,angle = α  =>  'i' component
------------------------------------------------------------------------------

Mα;

TMα;

p4a: p4; 

rotation : about y-axis   ,angle = β  =>  'j' component
------------------------------------------------------------------------------

Mβ:matrix([cos(β),0,sin(β)],[0,1,0],[-sin(β),0,cos(β)]);

TMβ : tracematrix(Mβ);

p4b:qβ =cos(β/2)+j*sin(β/2);

rotation : about z-axis   ,angle = γ  =>  'k' component
------------------------------------------------------------------------------

Mγ:matrix([cos(γ),-sin(γ),0],[sin(γ),cos(γ),0],[0,0,1]);

TMγ : tracematrix(Mγ);

p4c:qγ =cos(γ/2)+k*sin(γ/2);

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