FUNCTION CvT3d, degrees=Degrees, translate=Translate, scale=Scale,	$
	rotate=Rotate, xyexch= XYexch, xzexch = XZexch, yzexch = YZexch,	$
	vector=Vector, matrix=Matrix

;	PERSPECTIVE = pers, OBLIQUE = oblique

;+
; NAME:
;	CvT3d
; PURPOSE:
;	Implement three-dimensional transformation of coordinate systems.
;	!!! t3d provides a transformation matrix for transforming vectors in a
;		fixed coordinate frame. This routine provides a matrix for
;		transforming the coordinate frame, while keeping the vectors constant.
;	Accumulates one or more sequences of translation,
;	scaling, rotation, (perspective, and oblique) transformations
;	and returns the resulting 4x4 transformation matrix.
; CATEGORY:
;	Number strangling
; CALLING SEQUENCE:
;	pt = CvT3D(translate=T, scale=S, rotate=R, /degrees)
; OPTIONAL INPUTS:
;	/degrees	if set, all angles are assumed to be in degrees (default: radians)
;	matrix=Matrix		float array[4,4]
;				Used as starting transformation matrix
;				if absent the identity matrix is used as starting point.
;	vector=Vector		float array[3,*]
;				x,y,z coordinates in original coordinate frame.
;				If present, the transformation matrix is applied and the
;				result is returned, instead of the matrix itself.
;
;	All inputs to T3D are in the form of keywords.  Any, all, or none of
;	the following keywords can be present in a call to T3D.
;
;	The transformation specified by each keyword is performed in the
;	order of their descriptions below (e.g., if both TRANSLATE and
;	SCALE are specified, the translation is done first):
;
;	Translate	3-element vector of the translations in the X, Y, and Z directions.
;	Scale		3-element vector of scale factors for the X, Y, and Z axes.
;	Rotate		n-element vector of the rotations, about the X, Y, and Z axes.
;				The elements in Rotate are interpreted in sets of 3 for rotations
;				around X, Y, and Z axes.
;	XYexch		exchange the X and Y axes.
;	XZexch		exchange the X and Z axes.
;	YZexch		exchange the Y and Z axes.
;
; Not implemented
;	PERSPECTIVE	Perspective transformation.  This parameter is a scalar (p)
;		that indicates the Z distance of the center of the projection.
;		Objects are projected into the XY plane at Z=0, and the "eye"
;		is at point [0,0,p].
;	OBLIQUE:	A two-element vector of oblique projection parameters.
;		Points are projected onto the XY plane at Z=0 as follows:
;			x' = x + z(d COS(a)), and y' = y + z(d SIN(a)).
;		where OBLIQUE[0] = d, and OBLIQUE[1] = a.
; OUTPUTS:
;	If keyword vector is not specified, then the 4x4 transformation matrix is returned.
;	If keyword vertor is specified, then the transformation matrix is applied and
;		the result is returned as a float array[3,*]
; INCLUDE:
	@compile_opt.pro		; On error, return to caller
; CALLS:
;	ToRadians
; RESTRICTIONS:
;	This routine implements general rotations about the three axes.
; PROCEDURE:
; >	See also IDL procedure t3d.pro
; >	If a set of vectors is given in the form x[3,*] (containing
;	x,y,z-coordinates in the original coordinate frame), then
;	T3DMatrix(..)#x is an array[3,*] with coordinates in the new
;	coordinate frame.
; >	Using the /xyexch, /yzexch or /zxexch keywords turns a
;	right-handed coordinate frame into a left-handed frame. This impacts
;	the interpretation of the rotation angles (positive rotations are
;	counter-clockwise in a right-handed coordinate frame), and is better
;	avoided. These exchange keywords are only useful as final
;	transformations to swap the corresponding vector coordinates.
; MODIFICATION HISTORY:
;	DMS, Nov, 1987.
;	DMS, June 1990.	Fixed bug that didn't scale or translate
;			matrices with perspective properly.
;	DMS, July, 1996.  Added MATRIX keyword.
;	FEB-1998, Paul Hick (UCSD/CASS, pphick@ucsd.edu); based on T3D.pro
;-

rpd = ToRadians(degrees=Degrees)

id = fltarr(4,4)			; Make identity matrix
id[[0,1,2,3],[0,1,2,3]] = 1.0
							; Use input of start from identity
IF n_elements(Matrix) EQ 16 THEN a = Matrix ELSE a = id

; Difference with t3d:
; The translation is put in the bottom row, instead of the right column,
;		and are inserted with a minus sign.
; Scaling factors are entered as 1./Scale, instead of Scale
; New transformations are added to the left, instead of the right

IF n_elements(Translate) NE 0 THEN BEGIN
	ri = id
	ri[ [0,1,2], 3 ] = -Translate
	a = ri#a							; Apply translation
ENDIF

IF n_elements(Scale) NE 0 THEN BEGIN
	ri = id
	ri[ [0,1,2], [0,1,2] ] = 1./Scale
	a = ri#a							; Apply scale
ENDIF

IF n_elements(Rotate) NE 0 THEN BEGIN	; Rotation
	ri = id[0:2,0:2]					; Use 3 by 3 identity
	r  = ri  							; Accumulate in r
	sx = sin(Rotate*rpd)
	cx = cos(Rotate*rpd)

	FOR i=0,n_elements(Rotate)-1 DO BEGIN
		IF Rotate[i] NE 0 THEN BEGIN
			rr = ri						; 3x3 identity
			CASE i mod 3 OF
			0: BEGIN					; X axis
				rr[1,1] =  cx[i]
				rr[1,2] =  sx[i]
				rr[2,1] = -sx[i]
				rr[2,2] =  cx[i]
			END
			1: BEGIN					; Y axis
				rr[0,0] =  cx[i]
				rr[0,2] = -sx[i]
				rr[2,0] =  sx[i]
				rr[2,2] =  cx[i]
			END
			2: BEGIN					; Z angle
				rr[0,0] =  cx[i]
				rr[0,1] =  sx[i]
				rr[1,0] = -sx[i]
				rr[1,1] =  cx[i]
			END
			ENDCASE
			r = rr#r
		ENDIF
	ENDFOR

	rr = fltarr(4,4)					; Set up 4x4 matrix
	rr[0,0] = r
	rr[3,3] = 1.0

	a = rr#a							; Apply cumulative rot transforms
ENDIF

;IF n_elements(pers) NE 0 THEN BEGIN	; Perspective
;	r = id
;	r[2,3] = -1./pers
;	a = a # r
;ENDIF

;IF n_elements(oblique) NE 0 THEN BEGIN	; Oblique projection
;	r = id
;	r[2,2] = 0.0
;	r[2,0] = oblique[0]*cos(oblique[1]*rpd)
;	r[2,1] = oblique[0]*sin(oblique[1]*rpd)
;	a = a # r
;ENDIF

IF keyword_set(xyexch) THEN exch = [0,1]
IF keyword_set(xzexch) THEN exch = [0,2]
IF keyword_set(yzexch) THEN exch = [1,2]
IF n_elements(exch) NE 0 THEN BEGIN	; Exchange axes.
	t = a[exch[0],*]
	a[exch[0],0] = a[exch[1],*]
	a[exch[1],0] = t
ENDIF

CASE n_elements(vector) OF
0   : return, a
ELSE: return, (a#[vector, replicate(1.,1,n_elements(vector)/3)])[0:2,*]
ENDCASE

END