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

flight over rotating spherical earth (part1:rotation),version 1.0
peter.vlasschaert@gmail.com, 12/09/2018.

part a: 2D rotation , [,]
==================
==================

We need 2 unit vectors : 1e) a long x-axis : [1,0]  = e[x]
                                               2e) a long y-axis : [0,1]  = e[y]

unit : vector :
------------------

p1a:[1,0];
p1b:[0,1];

rotation : e[x]  = > e'[x]  ,  angle between = α = β
-----------------------------------------------------------------------

projection : e'[x] on the x-axis

p1a1: cos(α);

projection : e'[x] on the y-axis

p1a2:cos(π/2-α);

result :
======

p1a3:[p1a1,p1a2];

transpose :
-----------------

p1a4:transpose(p1a3);

rotation : e[x]  = > e'[x]  ,  angle between = α = β
-----------------------------------------------------------------------

projection : e'[y] on the x-axis , a long  negative : x-axis

p1b1: -sin(α);

projection : e'[y] on the y-axis , a long  positive : y-axis

p1b2: cos(α);

result :
======

p1b3:[p1b1,p1b2];

transpose :
-----------------

p1b4:transpose(p1b3);

unit        matrix : 2D :
-----------------------------

p1:matrix(p1a,p1b);

rotation matrix : 2D :
-----------------------------

p2:[[p1a4[1,1],p1b4[1,1]],[p1a4[2,1],p1b4[2,1]]];

p3:matrix(p2[1],p2[2]);

p3a1:p3[1,1];
p3a2:p3[1,2];
p3a3:p3[2,1];
p3a4:p3[2,2];

coordinates : initial vector (2D),'p4'. ( choose :vector)
-----------------------------------------------------
begin point : vector =>  (0,0) 
end    point  : vector =>  (x0,y0) , arrow 

p4:matrix([x0],[y0]);

p4a:transpose(p4);

coordinates : final vector (2D),'p4'. 
-----------------------------------------------------
begin point : vector =>  (0,0) 
end    point  : vector =>  (x1,y1) , arrow 

p5:matrix([x1],[y1]);

p5a:transpose(p5);

from : initial vector => final vector
============================

p5a1:p5[1][1];
p5a2:p5[2][1];

p6:p5a1= p3a1*x0+p3a2*y0;
p6a:p5a2= p3a3*x0+p3a4*y0;

Maxima : 
-------------

negdistrib;

p7:coefmatrix([p6,p6a],[x0,y0]);

Math : get same : 'p3'  = (-1)*Maxima.
-------------------------------------------------------

p7a:(-1)*p7;

find : determinant 'p7a' ?
=====================

p7b:determinant(p7a);
p7c:trigsimp(p7b);

           check : matrix 'p7a' :  [1,0]  => [1,0]
========================================
use : 'p6','p6a'. => result  x1=1,y1=0, when α = 0

p8:lhs(p6) =subst([x0=1,y0=0,α=0],rhs(p6));
p8a:lhs(p6a) =subst([x0=1,y0=0,α=0],rhs(p6a));

rotation matrix : 'p7a'
================================
inverse matrix ( p7a)  = transpose(p7a)

inverse matrix : 'p7a' ?

p9:invert(p7a);

p9a:trigsimp(p9);

transpose matrix : 'p7a' ?

p10:transpose(p7a);

check : is 'maxima' => true,false  ?

p11:p9a=p10;

p11a:is(p11);

rem : determinant (inversematrix= 'p10')

p12:determinant(p10);
p12a:trigsimp(p12);

part b: 3D rotation , [,,]
===================
===================

unit vectors 3D :
-------------------------

p13a:[1,0,0];
p13b:[0,1,0];
p13c:[0,0,1];

unit matrix 3D :
-----------------------

p14a:matrix(p13a,p13b,p13c);

rotation around x-axis

rotation around x-axis : e[y]  = > e'[y]  ,  angle between = α = β
--------------------------------------------------------------------------------------------
= > rotation  [0,1,0] around the  x-axis

projection : e'[y] on the y-axis

p15a1:cos(α);

projection : e'[y] on the z-axis

p15a2:cos(π/2-α);

result :
======

p15a3:[0,p15a1,p15a2];

transpose :

p15a4:transpose(p15a3);

rotation around x-axis : e[z]  = > e'[z]  ,  angle between = α = β
--------------------------------------------------------------------------------------------
= > rotation  [0,0,1] around the  x-axis

projection : e'[z] on the  y-axis

p16a1:-cos(π/2-α);

projection : e'[z] on the z-axis

p16a2:cos(α);

result :
======

p16a3:[0,p16a1,p16a2];

transpose :

p16a4:transpose(p16a3);

total : result

p17a1:p13a;
p17a2:p15a3;
p17a3:p16a3;

 total : rotation matrix: "rotation about  'x-axis' with angle = α"
<====================================================>

matrix(
 [0,0,0], 
 [0,0,0], 
 [0,0,0]
);

build rows : matrix
-----------------------------

p18a1:p17a1;
p18a2:[0,p17a2[2],p17a3[2]];
p18a3:[0,p17a2[3],p17a3[3]];

p19:matrix(p18a1,p18a2,p18a3);

rem : orientation, by cross product
e[x] cross e[ey ] into the plane, e[z]  is above the plane

<=================================================>

rotation about :
 y -axis,z-axis by angle = α , counter clockwise
---------------------------------------------------------------------------------------------

rotation :  about  y-axis with angle α : e'[x] =ex,e'[z]=ez
---------------------------------------------------------------------------------

p20a1:ex =[cos(α),0,-sin(α)];
p20a2:ey=[0,1,0];
p20a3:ez=[sin(α),0,cos(α)];

matrix : notation ,find rows :

p20a4:matrix(rhs(p20a1),rhs(p20a2),rhs(p20a3))$
p20a5:transpose(p20a4);

rotation :  about  z-axis with angle α : e'[x] =ex,e'[y]=ez
---------------------------------------------------------------------------------

p21a1:ex =[cos(α),sin(α),0];
p21a2:ey=[-sin(α),cos(α),0];
p21a3:ez=[0,0,1];

matrix : notation ,find rows :

p21a4:matrix(rhs(p21a1),rhs(p21a2),rhs(p21a3))$
p21a5:transpose(p21a4);

result : rot(α) x,rot(α) y,rot(α),z
==========================

how use : ''   ?  with function

p22:rota[x]=p19;
p23:rota[y]=p20a5;
p24:rota[z]=p21a5;

combination :rotations, principle 'stacked' use period 
====================================================================
first         rotation : rota[z] with α = α1
second  rotation : rota[y] with α = α2
third       rotation : rota[z] with α = α3
---------------------------------------------------------------------------------------------------------------------------
=>  first is last
rem : zyz ,called eulerian transformation between angles( = euler angles,α1,α2,α3)
---------------------------------------------------------------------------------------------------------------------------

p25a:subst([α = α1],rhs(p24));
p25b:subst([α = α2],rhs(p23));
p25c:subst([α = α3],rhs(p24));

total rotation :   tot[rotation]
=======================
use period = .  not star * ,matrix multiplication.

p26:tot[rota]= p25c.p25b.p25a;

use definition, '' not double quotes
=============================

p27:rotx(α):=''p19;
p28:roty(α):=''p20a5;
p29:rotz(α):=''p21a5;

result :example . (zyz) use period = .

p30:tot[rota]=rotz(α3).roty(α2).rotz(α1);

check : 'p26' and ' p30' are the same :true or false
=========================================

p31:rhs(p26)=rhs(p30)$
p32:is(p31);

rem :commands ' maxima ': transpose,invert.period,determinant,matrix,coefmatrix

------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------