//  Copyright (c) CNES  2008
//
//  This software is part of CelestLab, a CNES toolbox for Scilab
//
//  This software is governed by the CeCILL  license under French law and
//  abiding by the rules of distribution of free software.  You can  use,
//  modify and/ or redistribute the software under the terms of the CeCILL
//  license as circulated by CEA, CNRS and INRIA at the following URL
//  'http://www.cecill.info'.

function [q] = CL_rot_angles2quat(naxes,angles)
// Rotation angles to quaternion
//
// Calling Sequence
// q = CL_rot_angles2quat(naxes,angles)
//
// Description
// <itemizedlist><listitem>
// <p>Computes the quaternion <b>q</b> that defines the rotation resulting
// from the combination of successive elementary rotations, each elementary rotation 
// being described by an axis number (1=x-axis, 2=y-axis, 3=z-axis) and an angle. </p>
// </listitem>
// </itemizedlist>
//
// Parameters
// naxes : Axis numbers: 1=x-axis, 2=y-axis or 3=z-axis (1xP or Px1)
// angles : Rotation angles around respective axes [rad] (PxN) 
// q : quaternion (that defines the rotation) (4xN).
//
// Authors
// CNES - DCT/SB
//
// See also
// CL_rot_quat2matrix
// CL_rot_matrix2quat
//
// Examples
// // 40deg rotation around X followed by 60deg around Y: 
// q = CL_rot_angles2quat([1,2], [40;60] * %pi/180)
// 
// // Same considering individual rotations: 
// q1 = CL_rot_angles2quat(1, 40 * %pi/180); 
// q2 = CL_rot_angles2quat(2, 60 * %pi/180); 
// q = q1*q2 
// 


// Declarations:


// Code:

// validity checking
Nrot = length(naxes);
N = size(angles,2); 

// check: number of axes = number of angles
if (size(angles,1) <> Nrot)
  CL__error("Invalid arguments sizes"); 
end

// return [] if arguments are empty
if (Nrot == 0)
  q = CL_rot_defQuat([]); 
  return; 
end

// check that axis numbers are equal to 1, 2 or 3
I = find(naxes <> 1 & naxes <> 2 & naxes <> 3); 
if (~isempty(I))
  CL__error("Invalid axes number");
end

Id = eye(3,3); 

q = CL_rot_axAng2quat(Id(:,naxes(1)), angles(1,:));  
for k = 2 : Nrot
  q = q * CL_rot_axAng2quat(Id(:,naxes(k)), angles(k,:));
end

endfunction