function q = rotator(axis,angle)
%function q = rotator(axis,angle)
%inputs:  angle, real number
%          axis,   3X1 vector
%output    q       1X4 vector (unit quaternion)
t2 = angle/2;
cs = cos(t2); si = sin(t2);
nv = norm(axis);
  if nv == 0
      printf('need a nonzero axis vector \n')
      return
  end
  axis = axis/nv;
  q = si*axis';
  q = [cs, q];
  
  