function qm = quatmult(q1,q2)
%  quatmult computes the quaternion multiplication of q1 and q2
%
%  $BEGIN
%  $HEADER quatmult
%
%  quatmult computes the quaternion multiplication of q1 and q2
%
%  $HEADER REFERENCE
%  W. Bean, Operational Use of Quaternions
%  in the Description of Rigid Body Motion,
%  MSC Internal Note #69-FM-237, (1969)
%
%  $HEADER PROGRAMMER
%  John W. Craft, Jr.
%  University of Houston, (2005)
%
%  $HEADER VERSION
%  1.0.0
%
%  Example
%  my_quaternion = quatmult(q1,q2)
%
%  See Also quatmat quatrot
%
%  $END
%


 s3 = (q1(1)*q2(1) - q1(2)*q2(2)- q1(3)*q2(3) -q1(4)*q2(4));

 a3 = (q1(1)*q2(2)  +q1(2)*q2(1) +q1(3)*q2(4) -q1(4)*q2(3));

 b3 = (q1(1)*q2(3) - q1(2)*q2(4) +q1(3)*q2(1) +q1(4)*q2(2));
 
 c3 = (q1(1)*q2(4)  +q1(2)*q2(3) -q1(3)*q2(2) +q1(4)*q2(1));

 qm = [s3 a3 b3 c3];

return

