FUNCTION smei_cam2ccd, runit, camera=camera, degrees=degrees, distort=distort,	$
	inside=inside, boolean=boolean,	scalefov=scalefov,	$
	min_zunit=min_zunit

;+
; NAME:
;	smei_sky2ccd
; PURPOSE:
;	Convert from positions in the sky to CCD coordinates
; CATEGORY:
;	camera/idl
; CALLING SEQUENCE:
;	R = smei_sky2ccd( runit, [ camera=camera, /degrees] )
; INPUTS:
;	runit			array[3,n]; type: float
;						positions in the sky as unit vectors in the UCSD camera
;						frames (as output by smei_frm_sky2ccd)
;						These are converted to CCD coordinates if inside the FOV
;						If no points lie inside the
;						field of view then the scalar BadValue(rvec) is returned.
; OPTIONAL INPUT PARAMETERS:
;	camera=camera	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.
;
;	/degrees	if set, all in- and output angles are in degrees
;					(default: radians). Note that the 'camera' structure
;					must be consistent with the setting of /degrees.
;	/boolean	if set then the output array 'inside' is a byte array of zeros
;					and ones indicating whether the corresponding vectory is
;					outside/inside the field of view (instead of a list of 
;					indices from the IDL 'where' function identifying vectors
;					inside the field of view).
;	scalefov=scalefov
;				scalar; type: float; default=3.0 (covers essentially the whole
;					hemisphere centered on the camera optical axis).
;					Only locations closer then fov.size[0]/2 (= 30 degrees) to the
;					pole potentially lie inside the fov. The multiplier scalefov
;					extends the search to at most 89.5 degrees (essentially the
;					whole hemisphere centered on the pole. See PROCEDURE.
;	distort=distort
;				NOT IMPLEMENTED YET
; OUTPUTS:
;	R			array[2,n]; type: float
;					x and y pixel coordinates inside the SMEI fov
; OPTIONAL OUTPUT PARAMETERS:
;	inside=inside	array; type: long int (/boolean NOT set) or byte (/boolean set).
;						if /boolean is not set then 'inside' is a list of indices of
;						vectors inside field of view
;						if /boolean is set then 'inside' is a byte array with 0 for
;						vectors inside or 1  for vectors inside the field of view.
; INCLUDE:
	@compile_opt.pro		; On error, return to caller
; CALLS:
;	ToRadians, ToDegrees, InitVar, smei_camera, IsType
;	Inside_Wedge, BadValue, SyncDims, smei_theta2radial
; PROCEDURE:
;	The determination of sky locations inside the field of view is done in
;	two cuts. The first cut retains only points within scalefov*30 deg from
;	the optical axis. The default scalefov=3.0 so all points within 90 deg
;	are retained (actually 89.5 degrees).
;
;	The second cut retains only points strictly inside the field of view as
;	defined in the 'camera' structure. Only CCD coordinates for these points are
;	returned and are marked in the 'inside' array.
; MODIFICATION HISTORY:
;	FEB-2003, Paul Hick (UCSD/CASS)
;	FEB-2005, Paul Hick (UCSD/CASS)
;	    Added keyword runit
;	MAR-2005, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;		Split smei_sky2ccd into two parts: smei_sky2cam and smei_cam2ccd.
;-

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

InitVar, camera 	, 1
InitVar, boolean	, /key
InitVar, scalefov	, 3.0
InitVar, distort	, fltarr(6)

fov = smei_camera(camera=camera, /get_structure, degrees=degrees)

sz_2 = size(runit)
nvec = sz_2[sz_2[0]+2]/3
sz_2[1] = 2
sz_2[sz_2[0]+2] = nvec*2

r = reform(runit,3,nvec)

; Only locations closer then fov.size[0]/2 (= 30 degrees) to the pole
; potentially lie inside the fov. The multiplier scalefov extends the search to
; at most 89.5 degrees (essentially the whole hemisphere centered on the pole.
; CCD coordinates are returned only for points strictly inside the field of view.
; (in the corners of the fov the polar angle is slightly more than 30
; degrees, so maybe we need to add something to the fov extent used here??>

min_zunit = smei_camera(fov, /get_fov_size)
min_zunit = scalefov*(min_zunit[0,1]-min_zunit[0,0])/2
min_zunit = min_zunit < (89.9/dpm)
min_zunit = cos(min_zunit*rpm)
inside = where ( reform(r[2,*]) GE min_zunit, n_inside )

IF n_inside NE 0 THEN BEGIN

	r = r[0:1,inside]							; Drop the z-component.

	; Convert to theta-x, theta-y angles from Andy's memo (in degrees to set up call
	; to smei_theta2radial). Calculate the radial offset to the optical axis

	r = asin(r)*!radeg
	r[1,*] = smei_theta2radial(r[1,*], r[0,*]); Unit is pixels
	r[0,*] = r[0,*]/dpm 					; Convert angle to mystery units

	; Add optical axis parameters to get polar coordinates on the CCD

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

	; Final check for inside/outside fov. Only these are converted
	; to rectangular pixel coordinates.

	tmp = Inside_Wedge(smei_camera(fov,/get_fov_limits), r)
	tmp = where(tmp, n_inside)

	IF n_inside NE 0 THEN BEGIN
		inside = inside[tmp]

		; Convert points inside fov from polar to rectangular
		; and add the pixel coordinates of the fov center.

		r = smei_camera(fov, /get_center)#replicate(1,n_inside)+	$
			cv_coord(from_polar=r[*,tmp], /to_rect, degrees=degrees)
	ENDIF

ENDIF

out = make_array(type=IsType(r), dim=[2,nvec], value=BadValue(r) )
IF n_inside NE 0 THEN out[*,inside] = r
SyncDims, out, size=sz_2

IF boolean THEN BEGIN
	IF nvec EQ 1 THEN tmp = 0B ELSE tmp = bytarr(sz_2[2:sz_2[0]])
	IF n_inside NE 0 THEN tmp[inside] = 1B
	inside = tmp
ENDIF

RETURN, out  &  END