function w = quatmul(u,v)
%inputs: u,v quaternions to be multiplied   (1x4 arrays)
%output:  w = u^v their quaternion product  (1x4 array)
%   if    u = [u0, U], U=[u1,u2,u3]
%         v = [v0, V], V=[v1,v2,v3]
% then
%        W = [ u0*v0-U.V, u0*U+v0*V+UXV]
% where:   * is scalar multiplication
%          . is the dot product of U and V (=u1*v1+u2*v2+u3*v3)
%          X is the cross product of U and V (intrinsic)
%   UXV= [u2*v3-u3*v2,u3*v1-u1*v3,u1*v2-u2*v1]
% humor matlab's idiocyncracies by adding 1 to all indices
   U  = u(2:4);   V = v(2:4); 
   u0 = u(1);    v0 = v(1);
   w0 = u0*v0-dot(U,V);
   W  = u0*V+v0*U + cross1(U,V);
   w  = [w0,W];
   

function p = cross1(u,v)
%  p,a,b are 1X3 arrays
p = [u(2)*v(3)-u(3)*v(2),u(3)*v(1)-u(1)*v(3),u(1)*v(2)-u(2)*v(1)]; 
% Note that cross is an intrinsic function in matlab! (help cross)