share/contrib/quaternion.mac

   1 /* written by Gosei Furuya <go.maxima@gmail.com>
   2 #  This program is free software; you can redistribute it and/or modify
   3 #  it under the terms of the GNU General Public License as published by
   4 #  the Free Software Foundation; either version 2 of the License, or
   5 #  (at your option) any later version.
   6 */
   7
   8 dotdistrib:true$
   9 dotexptsimp:false $
  10 dotscrules:true $
  11 dotassoc:false$
  12 declare([i,j,k],nonscalar)$
  13
  14 defrule(f1,j.j,-1)$
  15 defrule(f2,i.i,-1)$
  16 defrule(f3,k.k,-1)$
  17 defrule(f4,i.j,k)$
  18 defrule(f5,j.i,-k)$
  19 defrule(f6,j.k,i)$
  20 defrule(f7,k.j,-i)$
  21 defrule(f8,k.i,j)$
  22 defrule(f9,i.k,-j)$
  23
  24 /*coeff4(a+b*i+c*j+d*k)--->[a,b,c,d]*/
  25 expand4(_ex):=block([_eex],_eex:apply2(expand(_ex),f1,f2,f3,f4,f5,f6,f7,f8,f9),
  26 _eex:apply2(_eex,f1,f2,f3,f4,f5,f6,f7,f8,f9),_eex)$
  27
  28 conj4(_ex1):=block([_eex1],_eex1:expand4(_ex1),
  29 _ex1:coeff(_eex1,i)*i+coeff(_eex1,j)*j+coeff(_eex1,k)*k,ratsimp(_eex1-2*_ex1))$
  30
  31 norm4(_ex1):= sqrt(expand4(_ex1 . conj4(_ex1)))$
  32
  33 inv4(_ex):= block( if (norm4(_ex)#0) then conj4(_ex)/norm4(_ex)^2 else false)$
  34
  35 scalarpart4(_ex1):=block([_eex1],_eex1:expand4(_ex1),
  36 _ex1:coeff(_eex1,i)*i+coeff(_eex1,j)*j+coeff(_eex1,k)*k,ratsimp(_eex1-_ex1))$
  37
  38 vectorpart4(_ex1):=block([_eex1],_eex1:expand4(_ex1),
  39 coeff(_eex1,i)*i+coeff(_eex1,j)*j+coeff(_eex1,k)*k)$
  40
  41 decomp4(_ex1):=block([_nn,_mm],_nn:norm4(_ex1),_mm:vectorpart4(_ex1),[acos(scalarpart4(_ex1)/_nn),_mm/norm4(_mm),_nn])$
  42
  43 invdecomp4(_exlist):=block([_th,_pureq,_r],_th:_exlist[1],_pureq:_exlist[2],_r:_exlist[3],expand4(_r*cos(_th)+_r*sin(_th)*_pureq))$
  44
  45 /*for matrix representation
  46 # matrix_element_mult:lambda([x,y],expand4(x.y))$
  47 # matrix([1,2*j],[3*i,4]).matrix([i,0],[j,k]);
  48 # result MATRIX([i-2,2*i],[4*j-3,4*k])
  49 # in this case a general inverse matrix exists.
  50 # for purpose of representing some matrix with basematrix of quaternion
  51 # a:matrix([1,i,j,k],[i,-1,-k,j],[j,k,-1,-i],[k,-j,i,-1])$
  52 # a.a.a is  MATRIX([-8,0,0,0],[0,-8,0,0],[0,0,-8,0],[0,0,0,-8])
  53 # so inverse a is -1/8*(a.a)
  54 # inva:-1/8*(a.a)$
  55 # we obtain coeff by inva.b,b is a matrix
  56 */
  57
  58 basematrix1:[matrix([0,-%i],[-%i,0]),matrix([0,-1],[1,0]),matrix([-%i,0],[0,%i])];
  59 basematrix2:[matrix([0,0,-1,0],[0,0,0,1],[1,0,0,0],[0,-1,0,0]),matrix([0,-1,0,0],[1,0,0,0],[0,0,0,-1],[0,0,1,0]),matrix([0,0,0,1],[0,0,1,0],[0,-1,0,0],[-1,0,0,0])];
  60 matcoeff4(_xx,_basematrix):=block([xx,bmat:_basematrix],expand(subst([xx=_xx,i=bmat[1],j=bmat[2],k=bmat[3]],[xx/4-(k . xx . k)/4-(j . xx . j)/4-(i . xx . i)/4,-(xx . i)/4-(k . xx . j)/4+(j . xx . k)/4-(i . xx)/4,-(xx . j)/4+(k . xx . i)/4-(j . xx)/4-(i . xx . k)/4,-(xx . k)/4-(k . xx)/4-(j . xx . i)/4+(i . xx . j)/4])))$
  61 matrix_element_mult:lambda([x,y],expand4(x.y));
  62 /*
  63 matcoeff4(matrix([1,-2*%i-3],[3-2*%i,1]),basematrix1);
  64 gg:matcoeff4(matrix([1,3],[-2,1]),basematrix1);
  65 */