function CombineRotations, R1, R2, R3, R4, R5, quaternion=quaternion, euler_angles=euler_angles, dcm=dcm, geometric=geometric, degrees=degrees

@compile_opt.pro		; On error, return to caller

;+
; NAME:
;	CombineRotations
; PURPOSE:
;	Combine two rotations, i.e. find the rotation resulting from
;	first applying one (R1), then another (R2) rotation.
; CATEGORY:
;	gen/idl/toolbox/math
; CALLING SEQUENCE:
;	R = CombineRotations(R1, R2 [,R3,R4,R5], /euler_angles, degrees=degrees)
; INPUTS:
;	R1				array[3,n], array[4,n], array[3,3,n]; type: float
;						first rotation
;	R2				array[3,m], array[4,m], array[3,3,m]; type: float
;						second rotation
;
;					R1 and/or R2 can be a single rotation with the other
;					multiple rotations, i.e. n = 1 and m != 1; or n != 1
;					and m = 1. The single rotation will be combined with
;					each of the multiple rotations.
;
;	R3,R4,R5		Additional rotations. Up to 5 rotations can be combined
;					in a single call.
; OPTIONAL INPUT PARAMETERS:
;	/quaternion 	input and output rotations are in quaternions;
;	/euler_angles	in- and output rotations are in Euler angles triplets
;	/dcm			in- and output rotations are in direction cosine matrices
;	/geometric		in- and putput rotations are in geometic form
;						(unit vector, plua rotation around unit vector).
;
;	/degrees		if set then all angles are in degrees (default: radians)
; OUTPUTS:
;	R				array[3,n > m]; type: same as input
;						combined rotation in same form as input
; CALLS:
;	InitVar, IsType, SyncArgs, CvRotation
; SEE ALSO:
;	EulerRotate
; PROCEDURE:
;	See definition of quaternions in: Foley, van Dam, Feiner and Hughes
;	'Computer graphics, principles and practice', 2nd edition.
;	Addison-Wesley, 1992, page 1063
; MODIFICATION HISTORY:
;	MAR-2003, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;-

nrot = n_params()

if nrot gt 2 then begin
	if nrot gt 5 then message, 'Cannot handle more than 5 rotations at a time'

	R = CvRotation(from_quaternion=[0,0,0,1], to_quaternion=quaternion, to_euler=euler_angles, to_dcm=dcm, to_geometric=geometric)

	switch 1 of
	IsType(R5, /defined): R = CombineRotations(R5, R, quaternion=quaternion, euler_angles=euler_angles, dcm=dcm, geometric=geometric, degrees=degrees)
	IsType(R4, /defined): R = CombineRotations(R4, R, quaternion=quaternion, euler_angles=euler_angles, dcm=dcm, geometric=geometric, degrees=degrees)
	IsType(R3, /defined): R = CombineRotations(R3, R, quaternion=quaternion, euler_angles=euler_angles, dcm=dcm, geometric=geometric, degrees=degrees)
	IsType(R2, /defined): R = CombineRotations(R2, R, quaternion=quaternion, euler_angles=euler_angles, dcm=dcm, geometric=geometric, degrees=degrees)
	IsType(R1, /defined): R = CombineRotations(R1, R, quaternion=quaternion, euler_angles=euler_angles, dcm=dcm, geometric=geometric, degrees=degrees)
	endswitch

	return, R
endif

InitVar, quaternion		, /key
InitVar, euler_angles	, /key
InitVar, dcm			, /key
InitVar, geometric		, /key

if quaternion+euler_angles+dcm+geometric eq 0 then euler_angles = 1

case 1 of
quaternion: begin
	q1 = R1
	q2 = R2
end
euler_angles: begin
	q1 = CvRotation(from_euler=R1, /to_quaternion, degrees=degrees)
	q2 = CvRotation(from_euler=R2, /to_quaternion, degrees=degrees)
end
dcm: begin
	q1 = CvRotation(from_dcm=R1, /to_quaternion)
	q2 = CvRotation(from_dcm=R2, /to_quaternion)
end
geometric: begin
	q1 = CvRotation(from_geometric=R1, /to_quaternion, degrees=degrees)
	q2 = CvRotation(from_geometric=R2, /to_quaternion, degrees=degrees)
end
endcase

SyncArgs, q1, q2

dims = size(q1, /dim)
nrot = n_elements(q1)/dims[0]		; dims[0] is 4

q1 = reform(q1, [dims[0],nrot], /overwrite)
q2 = reform(q2, [dims[0],nrot], /overwrite)


									; Scalar in last position: q = [qx,qy,qz,q0]
R = [	q1[3,*]*q2[0,*]+q1[0,*]*q2[3,*]+q1[1,*]*q2[2,*]-q1[2,*]*q2[1,*] ,	$
		q1[3,*]*q2[1,*]+q1[1,*]*q2[3,*]+q1[2,*]*q2[0,*]-q1[0,*]*q2[2,*] ,	$
		q1[3,*]*q2[2,*]+q1[2,*]*q2[3,*]+q1[0,*]*q2[1,*]-q1[1,*]*q2[0,*] ,	$
		q1[3,*]*q2[3,*]-q1[0,*]*q2[0,*]-q1[1,*]*q2[1,*]-q1[2,*]*q2[2,*] ]
									; Scalar in first position: q = [q0,qx,qy,qz]
;R = [	q1[0,*]*q2[0,*]-q1[1,*]*q2[1,*]-q1[2,*]*q2[2,*]-q1[3,*]*q2[3,*],	$
;		q1[0,*]*q2[1,*]+q1[1,*]*q2[0,*]+q1[2,*]*q2[3,*]-q1[3,*]*q2[2,*],	$
;		q1[0,*]*q2[2,*]+q1[2,*]*q2[0,*]+q1[3,*]*q2[1,*]-q1[1,*]*q2[3,*],	$
;		q1[0,*]*q2[3,*]+q1[3,*]*q2[0,*]+q1[1,*]*q2[2,*]-q1[2,*]*q2[1,*] 	]

R = reform(R, dims, /overwrite)

case 1 of
quaternion	:
euler_angles: R = CvRotation(from_quaternion=R, /to_euler	 , degrees=degrees)
dcm 		: R = CvRotation(from_quaternion=R, /to_dcm 	 )
geometric	: R = CvRotation(from_quaternion=R, /to_geometric, degrees=degrees)
endcase

return, R  &  end