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

rotation of plane 'α'  about axis 'w' with angle 'θ', version 1.0
peter.vlasschaert@gmail.com, 06/10/2018.

r0: [x0,y0,z0] ,position vector from O to L, 'see Figure 1 '

def :  'α' plane ( = HIJK), 3D
======================

p1: a*x+b*y+c*z+d = 0;

coef: 'p1' , x,y,z,constant 
--------------------------------------

p1a:coeff(lhs(p1),x);
p1b:coeff(lhs(p1),y);
p1c:coeff(lhs(p1),z);
p1d:subst([x=0,y=0,z=0],lhs(p1));

point '  L', 'x[0],y[0],z[0] 'with vector 'r0'
=================================

p2:[x=x0,y=y0,z=z0];
p2a:subst(p2,p1);

How to find Normal (= LM) to the Plane ('α')
====================================

p3:p1-p2a;

coef: 'p3' , a,b,c
--------------------------------------

p3a:coeff(lhs(p3),a);
p3b:coeff(lhs(p3),b);
p3c:coeff(lhs(p3),c);

p4:p3a*a+p3b*b+p3c*c = rhs(p3);

p4a:[a,b,c];
p4b:transpose([p3a,p3b,p3c]);

scalar product : 0 , see 'result'
----------------------------------------------

p4c:p4a.p4b =0 ;

result : vector 'p4a' perpendicular to vector plane 'p4b'.
-----------------------------------------------------------------------------------------
p4a =>  vector (LM) , see Figure 1 ' transpose('p4a')'
p4b =>  vector (p4b), see Figure 1 , vector inside the plan 'α'
------------------------------------------------------------------------------------------
rem : vector(LM) call the normal vector the plane. ('p4d')

p4d:transpose(p4a);

How to find rotation axis  'w' , 'see Figure1'
====================================

solution : through the cross product : 'p4d' X unit vector (=ez)

1e) step : unit vectors form unit matrix

p5a:[1,0,0];
p5b:[0,1,0];
p5c:[0,0,1];
p5:unit[matrix] =matrix(p5a,p5b,p5c);

call : ex,ey,ez , 'see Figure 1'

p5a1:ex=transpose(p5a);

p5a1:ex=transpose(p5a);
p5b1:ey=transpose(p5b);
p5c1:ez=transpose(p5c);

result :columns of unit matrix

2e) cross product , maxima , do computation 

a1)  Maxima
-----------------

p6:axis= p4a  ~ p5c;

p6a:express(p6);

with use , package 'load(vect)'
 b  component  'ex'
-a  component  'ey'
0   component  'ez'

a2) Math
-----------

p7:axis=matrix(
 [lhs(p5a1),lhs(p5b1),lhs(p5b1)], 
 [p1a,p1b,p1c], 
 [0,0,1]
);

p7a:lhs(p7)=determinant(rhs(p7));

a1 and a2 ,same 'result'

How to find angle 'θ' , angle between 'z axis' and 'normal to the plane' (see Figure1)
======================================================================

Formula :  scalar product
=====================
cos(θ) =  ?

1e) unit normal vector to the plane 'α'
-------------------------------------------------------

p8:q = sqrt((p4a[1])^2+(p4a[2])^2+(p4a[3])^2);

p8a:subst(p8,p4a/q);

2e) ez vector (unit vector)
--------------------------------------

p9:p5c;

3e) result : cos(θ)
---------------------------

p10:cos(θ) =p8a.p5c;

4e) solve θ : ?
---------------------

p10a:solve(p10,θ);

p10b:p10a[1];

How use : given  a,b,c,d  from equation of the plane 'α'
----------------------------------------------------------------------------------
=>  a*x+b*y+c*z+d = 0

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