mulQ:=proc(xx,yy)
a1:=op(1,xx); b1:=op(2,xx); c1:=op(3,xx); d1:=op(4,xx);
a2:=op(1,yy); b2:=op(2,yy); c2:=op(3,yy); d2:=op(4,yy);
pp:= (a1*a2 - b1*b2 - c1*c2 - d1*d2);
qq:= (a1*b2 + b1*a2 + c1*d2 - d1*c2);
rr:= (a1*c2 - b1*d2 + c1*a2 + d1*b2);
ss:= (a1*d2 + b1*c2 - c1*b2 + d1*a2);
[pp,qq,rr,ss];
end;
invQ:=proc(xx)
a1:=op(1,xx); a2:=op(2,xx); a3:=op(3,xx); a4:=op(4,xx);
[a1,-a2,-a3,-a4];
end;
rotQ:=proc(axis1,angle)
cc:=cos(angle/2);
ss:=sin(angle/2);
x1:=op(1,axis1);
y1:=op(2,axis1);
z1:=op(3,axis1);
tt := sqrt(x1*x1+y1*y1+z1*z1);
xx:=x1/tt;
yy:=y1/tt;
zz:=z1/tt;
evalf([cc,ss*xx,ss*yy,ss*zz]);
end;
a:=rotQ([1,2,3],Pi/3.0);
mulQ(a,invQ(a));
unassign('p','q','r','s','x','y','z','w');
unassign('a','b','c','d','x','y','z','w');
X1:=mulQ([p,q,r,s],[0,x,y,z]);
X :=mulQ(X1,[p,-q,-r,-s]);
X1:=mulQ([a,b,c,d],[0,x,y,z]);
X :=mulQ(X1,[a,-b,-c,-d]);
simplify(op(1,X));
simplify(op(2,X));
simplify(op(3,X));
simplify(op(4,X));
JSFH