FUNCTION smei_ccd2sky, rvec, camera=camera, quaternion=quaternion,	$
	qcamera=qcamera, degrees=degrees, inside=inside,				$
	boolean=boolean, distort=distort, rectangular=rectangular

;+
; NAME:
;	smei_ccd2sky
; PURPOSE:
;	Convert from CCD coordinates to positions in the sky
; CATEGORY:
;	SMEI
; CALLING SEQUENCE:
;	R = smei_ccd2sky( rvec, [ camera=camera, quaternion=quaternion, /degrees])
; INPUTS:
;	rvec			array[2,*]; type: float
;						x- and y- coordinates of locations on the CCD
;						These are converted to sky coordinates if inside the
;						field of view. If no points lie inside the field of
;						view then the scalar BadValue(rvec) is returned.
;						If rvec is not specified then all pixels in a 
;						frame are processed.
; OPTIONAL INPUT PARAMETERS:
;	/rectangular	if set then the return value is in rectangular coordinates
;					(with 1st dimension of 3).
;
;	camera=camera 	This keyword can be any one of the following three.
;					array[1]; type: smei_buf_frm structure
;						this determines the camera and the quaternion.
;						The keywords 'quaternion' and 'qcamera' are ignored.
;					array[1]; type: smei_camera_fov structure
;						structure containing information about SMEI field of
;						view; usually the return value of href=smei_camera=
;						with the /get_structure keyword set.
;					scalar; type: integer
;						camera number (1,2 or 3).
;					If not specified then camera=1 is assumed.
;
;					The camera number is used as argument to smei_camera to
;					obtain the fov structure.
;
;	quaternion=quaternion
;					array[4]; type: float
;						Coriolis quaternion, or (if /qcamera is set) the
;						camera quaternion (i.e. the quaternion
;						constructed by href=smei_cam_quaternion=).
;	/qcamera		indicates that 'quaternion' already is the camera
;					quaternion for the camera specified in 'camera'.
;	/degrees		if set, all in- and output angles are in degrees
;						(default: radians). Note that the 'fov' structure
;						must be consistent with the setting of /degrees.
;	/boolean		if set then the output array 'inside' is returned as an
;					2D array of indices. (set this keyword only if rvec
;					is NOT specified, i.e. when processing a whole frame.
;	distort=distort
;					NOT IMPLEMENTED YET
; OUTPUTS:
;	R			array[2,n]; type: float
;					spherical coordinates (longitude, latitude)
;					for all pixels inside the SMEI fov
; OPTIONAL OUTPUT PARAMETERS:
;	inside=inside
;				array[  n]; type: long integer (if /twod NOT set)
;				array[2,n]; type: long integer (if /twod set0
;					one- or two-D CCD pixel indices for all the pixels
;					listed in R
; INCLUDE:
	@compile_opt.pro		; On error, return to caller
; CALLS:
;	ToRadians, ToDegrees, InitVar, smei_camera, IsType, gridgen, AngleRange
;	Inside_Wedge, boost, EulerRotate, SubArray, BadValue, SyncDims
;	smei_radial2theta, smei_cam_quaternion, CvRotation
; PROCEDURE:
; MODIFICATION HISTORY:
;	OCT-1998, Paul Hick (UCSD/CASS)
;	FEB-2001, Paul Hick (UCSD/CASS; pphick@ucsd.edu); added quaternion keyword
;-

rpm = ToRadians(degrees=degrees)
dpm = ToDegrees(degrees=degrees)

InitVar, rectangular, /key
InitVar, boolean	, /key
InitVar, qcamera	, /key
InitVar, distort	, intarr(6)

CASE 1 OF
IsType(camera,/structure): begin
	CASE tag_names(camera, /structure_name) OF
	'SMEI_FRM_HDR': BEGIN
		fov = smei_camera(camera=camera.camera, /get_structure, degrees=degrees)
		quaternion = camera.quaternion
	END
	'SMEI_CAMERA_FOV': fov = camera
	END
END
IsType(camera,/undefined): fov = smei_camera(camera=1, /get_structure, degrees=degrees)
ELSE					 : fov = smei_camera(camera=camera, /get_structure, degrees=degrees)
ENDCASE

; If no input rvec is specified set up a 2D grid of x,y positions across
; the entire CCD, with the origin at the position set by 'fov'
; If rvec is specified (in zero-based pixel coordinates) then shift
; the origin to the 'fov' position.

CASE IsType(rvec, /defined) OF
0: r = gridgen( smei_camera(fov, /get_nsize), orig=smei_camera(fov, /get_center) )
1: r = rvec-smei_camera(fov, /get_center)#replicate(1,n_elements(rvec)/2)
ENDCASE

sz_rvec = size(r)
npart = sz_rvec[1]		; Always 2
nvec  = n_elements(r)/npart

r = reform(r, npart, nvec, /overwrite)
 
; Convert to polar coordinates
; Make sure the phase angle is in [-!pi,+!pi] before calling Inside_Wedge
; (Why??? PPH 02/2003)

r = cv_coord(from_rect=r, /to_polar, degrees=degrees)
r[0,*] = AngleRange(r[0,*], /pi, degrees=degrees)

; Only process positions inside the fov. 'inside' is zero outside and 1
; inside the fov.

inside = Inside_Wedge(smei_camera(fov, /get_fov_limits), r)

;inside = reform(inside, fov.nsize , /overwrite)
;if n_elements(frame) eq 0 then frame = inside
;view, in=frame, /stretch, /bl

; Only retain coordinates inside fov

inside = where(inside, n_inside)

IF n_inside NE 0 THEN BEGIN

	r = r[*,inside] 					; Retain coordinates inside fov only.

	; r now contains polar coordinates on the CCD relative to fov center.
	; Convert to angles in the sky relative to the optical axis.

	r = r-smei_camera(fov,/get_optical_axis)#replicate(1,n_inside)

	; Convert the radial offset to an angle thetay in the sky.
	; smei_radial2theta uses angles in degrees only

	r[0,*] = r[0,*]*dpm
	r[1,*] = smei_radial2theta( r[1,*], r[0,*], tol=1.0E-5 )

	r = sin(r/!radeg)					; Direction cosines in x-y plane
	boost, r, sqrt((1-total(r*r,1)) > 0); Make unit vector in rectangular coordinates

	IF IsType(quaternion, /defined) THEN BEGIN

		; Apply the quaternion to rotate to sky coordinates.

		CASE qcamera OF
		0: dcm = smei_cam_quaternion(quaternion, fov.camera)
		1: dcm = quaternion
		ENDCASE

		dcm[0:2] = -dcm[0:2]		; Invert quaternion
		dcm = CvRotation(from_quaternion=dcm, /to_dcm)
		r = transpose(dcm)#r

	ENDIF

	; Convert to spherical coordinates (phase angle and latitude)
	; Drop the distance (is 1.0)

	IF NOT rectangular THEN	$
		r = (cv_coord(from_rect=r, /to_sphere, degrees=degrees))[0:1,*]
	
ENDIF

npart = 2+rectangular
out = make_array(type=IsType(r), dim=[npart,nvec], value=BadValue(r) )

IF n_inside NE 0 THEN out[*,inside] = r

IF boolean THEN BEGIN
	CASE IsType(rvec, /defined) OF
	0: tmp = bytarr(smei_camera(fov, /get_nsize))
	1: tmp = byte(SubArray(rvec,element=0)*0)
	ENDCASE

	IF n_inside EQ 0 THEN inside = tmp ELSE inside = tmp[inside]
ENDIF

sz_rvec[sz_rvec[0]+2] = sz_rvec[sz_rvec[0]+2]/sz_rvec[1]*npart 
sz_rvec[1] = npart	; From 2 --> 2 or 3

SyncDims, out, size=sz_rvec

RETURN, out  &  END