function [a,b,c,d] = quatern(w,x,y,z)
# quatern: construct or extract a quatern
#  w = a*i + b*j + c*k + d from given data.
#
# calling formats:
# [a,b,c,d]   = quatern(w)		-or-
# [vv,theta] = quatern(w)
# w           = quatern(a,b,c,d)
# w           = quatern(vv,theta)

switch(nargin)
case(1),					# extract information
  if(!(is_vec(w) & length(w) == 4) )
    error("input vector must be of length 4)")
  endif
  # extract data
  switch(nargout)
  case(4),
    a = w(1);
    b = w(2);
    c = w(3);
    d = w(4);
  case(2),
    if(abs(norm(w)-1) > 1e-12)
      warning(sprintf("quatern: ||w||=%e, setting=1 for vv, theta",norm(w)))
      w = w/norm(w);
    endif
    [a,b,c,d] = quatern(w);
    theta = acos(d)*2;
    if(abs(theta) > pi)
      theta = theta - sign(theta)*pi;
    endif
    sin_th_2 = norm([a, b, c]);

    if(sin_th_2 != 0)
      vv = [a,b,c]/sin_th_2;
    else
      vv = [a, b, c];
    endif
    a = vv;
    b = theta;
  otherwise,
    usage("[a,b,c,d] = quatern(w) or [vv,theta] = quatern(w)");
  endswitch

case(2),
  if(nargout != 1)
    usage("w = quaterion(vv,theta)");
  endif
  vv = w;
  theta = x;

  if(!is_vec(vv) | length(vv) != 3)
    error("vv must be a length three vector")
  elseif(!is_scal(theta))
    error("theta must be a scalar");
  elseif(norm(vv) == 0)
    error("quatern: vv is zero.") 
  elseif(abs(norm(vv)-1) > 1e-12)
    warning("quatern: ||vv|| != 1, normalizing")
    vv = vv/norm(vv);
  endif

  if(abs(theta) > 2*pi)
    warning("quatern: |theta| > 2 pi, normalizing")
    theta = rem(theta,2*pi);
  endif
  vv = vv*sin(theta/2);
  d = cos(theta/2);
  a = quatern(vv(1), vv(2), vv(3), d); 

case(4),
  if(nargout != 1)
    usage("w = quaterion(a,b,c,d)");
  endif
  if ( !(is_scal(w) & is_scal(x) & is_scal(y) & is_scal(z)) )
    error("input values must be scalars.")
  endif
  a = [w, x, y, z];

otherwise,
  error("usage: [a,b,c,d] = quatern(w), a = quatern(w,x,y,z)")

endswitch

endfunction
