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

crossproduct and unit vectors , version 1.0
peter.vlasschaert@gmail.com,02/09/2018

declare([A,B,C],nonscalar);

unit vectors :  'see below'
===========

p1:[e[x],e[y],e[z]];
p2:subst([x=θ,y=r],p1);
p3:subst([x=t,y=n],p1);

use : operator  ' ~ ' ?
=================================
math        : cross-product ( = outer product)  =      p1 cross p2
maxima   : cross-product ( = outer product)  =  - (p2 cross p1) or (p1 cross p2)

see : example 
---------------------

p1a:p1~p2;

p1b:p1~p3;

p1c:-1*(p3~p1);

rem : minus sign not found ' express' , gives the components 

p1a1:express(p1a);

p1b1:express(p1b);
p1c1:express(p1c);

why : r => n, t => θ, z=z,   relative position hold , see  'p2,p3'

p2a:express(A);
p2b:subst([x=θ,y=r],p2a);
p2c:subst([x=t,y=n],p2a);

vector A : 
========

system : x,y,z  , (cartesian coordinates); 
------------------------------------------------------------
unit vector =>  [e[x],e[y],e[z]]

p3: A=sum(express(A)[i]*p1[i],i,1,3);

scalar product ( = outer product )

p31:A^2 = rhs(p3)*rhs(p3);

p32:lhs(p31) = expand(rhs(p31));

p33a:let(e[x]*e[x],1);
p33b:let(e[y]*e[y],1);
p33c:let(e[z]*e[z],1);
p33d:let(e[x]*e[y],0);
p33e:let(e[y]*e[z],0);
p33f:let(e[x]*e[z],0);

p34: letsimp(p32);

p341:ratsimp(p34);

p342:expand(p341);

p342a:sqrt(p342);

p351:[p3*e[z],p3*e[y],p3*e[x]];
p352a:p351[1];
p352b:p351[2];
p352c:p351[3];
p353:ratsimp(p351);
p354:expand(p353);
p355:letsimp(p354);
p356:ratsimp(p355);
p357:solve(p356,[e[z],e[y],e[x]])[1];

p358:A*subst(p357,p3);

p359:sqrt(ratsimp(p358));

p360:A=rhs(p359);

p361:p357;

p362a:A*p361[1];
p362b:A*p361[2];
p362c:A*p361[3];

p363a:rhs(p362a)=subst([p360],lhs(p362a));
solve(p363a,e[z]);
p363b:rhs(p362b)=subst([p360],lhs(p362b));
solve(p363b,e[y]);
p363c:rhs(p362c)=subst([p360],lhs(p362c));
solve(p363c,e[x]);

p364:p3;

rem  : distance = A = abs(A).

system: θ,r,z , (polar coordinates)
---------------------------------------------------
unit vector => [e[θ],e[r],e[z]]

p3a:A=sum(p2b[i]*p2[i],i,1,3);

system: n,t,z , (normal,tangential system)
--------------------------------------------------------------
unit vector => [e[t],e[n],e[z]]

p3b:A=sum(p2c[i]*p2[i],i,1,3);

rem : appendix :  2D applications , A[z] > 0 , above the plane
                                                                  A[z] < 0, below the plane
==================================================

let(cross(e[x],e[y]),e[z]);
let(cross(e[y],e[x]),-1*e[z]);
let(cross(e[x],e[x]),0);
let(cross(e[y],e[y]),0);

p4:cross(e[x],e[y])+cross(e[y],e[x]);

p4a:letsimp(p4);

use :cross1 ,  same  as 'let'  with  cross
=================================

cross1(a,b) :=block([c1:e[z]],if a=b  then c:0,if a ≠ b and (a=e[x],b=e[y])  then c:c1 else c:-c1 ,return(c));

p5:cross1(d,d);

p5a:cross1(a,a);

p5b:cross1(e[x],e[y]);

p5c:cross1(e[y],e[x]);

p5d:cross1(e[x],e[y])+cross1(e[y],e[x]);

definition :cross2
===============

cross2(a,b,c):=block([], 
           if a=e[x] and b=e[x]  then c:0,  
           if b=e[x] and a=e[x]  then c:0,
           if a=e[y] and b=e[y]  then c:0 ,
           if b=e[y] and a=e[y]  then c:0,
           if a=e[z] and b=e[z]  then c:0 ,  
           if b=e[z] and a=e[z]  then c:0,
           if a=e[x] and b=e[y]  then c:e[z], 
           if a=e[y] and b=e[z]  then c:e[x],  
           if a=e[z] and b=e[x]  then c:e[y] ,         
           if a=e[y] and b=e[x]  then c:-e[z], 
           if a=e[z] and b=e[y]  then c:-e[x],  
           if a=e[x] and b=e[z]  then c:-e[y] ,  
        return(c));
                   
                  

p6:cross2(e[x],e[y],c);

p6a:cross2(e[y],e[x],c);

p6b:cross2(e[y],e[y],c);

p6c:cross2(e[y],e[z],c);

p6d:cross2(e[x],e[z],c);

definition :cross3
===============

cross3(a,b):=block([], 
if (a=b and b=a) then c:0,  
if a=e[x] and b=e[y]  then c:e[z],  
if a=e[y] and b=e[z]  then c:e[x] , 
if a=e[z] and b=e[x]  then c:e[y] ,           
if a=e[y] and b=e[x]  then c:-e[z], 
if a=e[z] and b=e[y]  then c:-e[x],  
if a=e[x] and b=e[z]  then c:-e[y] ,  
        return(c));                  

p7:cross3(e[y],e[y]);

p7a:cross3(e[z],e[y]);

p7b:cross3(e[x],e[z]);

definition :cross4
===============
x=> y=> z=> x , cyclic permutation
result :
------------------------------
x => y , y => z, z=> x

cross4(a,b):=block([], 
if (a=b and b=a) then 
         c:0,  
if (a ≠ b and b ≠ a)  then   
        if a=e[x] and b=e[y]  then c:e[z]
               elseif       
        a=e[y] and b=e[x]  then c:-e[z],  
     
       if a=e[y] and b=e[z]  then c:e[x]         
               elseif      
       a=e[z] and b=e[y]  then c:-e[x],  
          
       if a=e[z] and b=e[x]  then c:e[y]
               elseif
       a=e[x] and b=e[z]  then c:-e[y],
        return(c));          

p8:cross4(e[x],e[z]);

p8a:cross4(e[y],e[z]);

p8b:cross4(e[z],e[z]);

p8c:cross4(e[y],e[x]);

p8d:cross4(e[z],e[x]);

definition :cross5 
===============
principle : cyclic permutation
========================
d=1,  system :[e[x],e[y],e[z]], x => y , y => z, z=> x
d=2,  system  :[e[θ],e[r],e[z]] ,θ => r,  r => z, z => θ
d=3,  system  :[e[t],e[n],e[z]], t =>  n, n => z,z=> t

cross5(a,b,d):=block([], 
if (a=b and b=a) then 
         c:0,  

    if d=1 then    
        if a=e[x] and b=e[y]  then c:e[z]
               elseif       
        a=e[y] and b=e[x]  then c:-e[z],  
     
       if a=e[y] and b=e[z]  then c:e[x]         
               elseif      
       a=e[z] and b=e[y]  then c:-e[x],  
          
       if a=e[z] and b=e[x]  then c:e[y]
               elseif
       a=e[x] and b=e[z]  then c:-e[y],          
if d=2 then
       if a=e[θ] and b=e[r]  then c:e[z]
                     elseif       
        a=e[r] and b=e[θ]  then c:-e[z],  
     
       if a=e[r] and b=e[z]  then c:e[θ]         
               elseif      
       a=e[z] and b=e[r]  then c:-e[θ],  
          
       if a=e[z] and b=e[θ]  then c:e[r]
               elseif
       a=e[θ] and b=e[z]  then c:-e[r],
if d=3 then   
         if a=e[t] and b=e[n]  then c:e[z]
                     elseif       
        a=e[n] and b=e[t]  then c:-e[z],  
     
       if a=e[n] and b=e[z]  then c:e[t]         
               elseif      
       a=e[z] and b=e[n]  then c:-e[t],  
          
       if a=e[z] and b=e[t]  then c:e[n]
               elseif
       a=e[t] and b=e[z]  then c:-e[n],
        
    
        return(c));          

p9:cross5(e[x],e[z],1);

p9a:cross5(e[z],e[z],1);

p9b:cross5(e[θ],e[r],2);

p9c:cross5(e[θ],e[θ],2);

p9d:cross5(e[r],e[θ],2);

p9e:cross5(e[z],e[t],3);

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