;+
; NAME:
;	cvsmei
; PURPOSE:
;	Convert positions between various SMEI-related coordinate frames
; CATEGORY:
;	camera/idl/toolbox
; CALLING SEQUENCE:
	FUNCTION cvsmei, $
		camera		= camera	 , $
		quaternion	= quaternion , $

		mode		= mode		 , $
		from_mode	= from_mode  , $
		to_mode 	= to_mode	 , $

		from_ccd	= from_ccd	 , $
		from_camera = from_camera, $
		from_equatorial    = from_equatorial   , $
		from_helioecliptic = from_helioecliptic, $
		from_map	= from_map	 , $

		time		= time		 , $

		rectangular = rectangular, $
		degrees 	= degrees	 , $
		force_unit	= force_unit , $

		to_ccd		= to_ccd 	 , $
		to_camera	= to_camera	 , $
		to_equatorial = to_equatorial, $
		to_helioecliptic = to_helioecliptic, $
		to_map		= to_map	 , $

		inside		= inside	 , $
		boolean 	= boolean	 , $
		ndim		= ndim		 , $
		silent		= silent	 , $
		count		= count
; INPUTS:
; OPTIONAL INPUT PARAMETERS:
;   from_map        array[2,n]; type: numerical
;                   scalar; type: string
;
;						If set to string 'MAP' then the position of all
;						pixels in the skymap for the specified mode are used.
;
;                       Pixel locations in a skymap, specified as a zero-base
;                       array index (i.e. lower-left pixel is [0,0].
;
;                       The mode is "lores" (low-resolution skymap),
;                       "eqmap" (high-resolution equatorial map),
;                       "npmap" (hi-res map of north pole) or "spmap"
;                       (hi-res map of south pole). This is set using
;                       keywords "mode" or "from_mode"
;
;	from_equatorial	array[2,n] or array[3,n]
;						positions in the sky as spherical or rectangular
;						equatorial coordinates.
;
;					/rect NOT SET: spherical coordinates
;						rvec[0,*] is RA;
;						rvec[1,*] is declination
;						rvec[3,*] (usually the distance)
;
;					/rect SET: rectangular coordinates
;						rvec[0,*] is x-coordinate (toward vernal equinox)
;						rvec[1,*] is y-coordinate
;						rvec[3,*] is z-coordinate (toward equatorial north)
;
;	from_helioecliptic
;					array[2,n] or array[3,n]
;						positions in the sky as spherical or rectangular
;						helioecliptic coordinates
;
;					/rect NOT SET: spherical coordinates
;						rvec[0,*] is sun-centered ecliptic longitude;
;						rvec[1,*] is ecliptic longitude
;						rvec[3,*] (usually the distance)
;
;					/rect SET: rectangular coordinates
;						rvec[0,*] is x-coordinate (toward sun)
;						rvec[1,*] is y-coordinate
;						rvec[3,*] is z-coordinate (toward ecliptic north)
;
;	from_camera 	array[3,n]
;						positions in the sky as unit vectors in the
;						UCSD camera frame.
;
;	from_ccd		array[2,n]; type: double
;					scalar; type: string
;
;						If set to string 'CCD' then the position of all
;						pixels in the image for the specified mode are used.
;
;						pixel locations on the CCD, specified as a zero-based
;						array index (i.e. lower-left pixel is [0,0])
;						The unit is the pixel size for mode 0.
;						Override this by explicitly setting the 'mode' keyword
;						(so set mode=0 for engineering mode pixels).
;
;	/to_map			returns array[2,n]; type: double
;						output is in pixel coodinates in a sky map
;						selected with the "mode" of "to_mode" keywords
;	/to_equatorial	returns array[3,n] or array[2,n]; type: double
;						output is in equatorial sky coordinates (rectangular
;						array[3,n ] if /rectangular is SET.
;						array[2,n] (RA and dec) if /rectangular NOT set.
;	/to_helioecliptic
;			 		returns array[3,n] or array[2,n]; type: double
;						output is in helioecliptic sky coordinates (rectangular
;						array[3,n ] if /rectangular is SET.
;						array[2,n] (sun-centered ecliptic longitude and latitude
;						if /rectangular NOT set.
;	/to_camera		returns array[3,n]; type: double
;						output is a unitvector in the UCSD camera frame
;						(with z-axis along camera optical axis)
;	/to_ccd 		returns array[2,n]; type: double
;						output is in pixel coordinates on the CCD
;						The units are mode-dependent.
;
;	/rectangular	only used if 'from_equatorial' or /to_equatorial is set
;						if set then 'from_equatorial' is in rectangular
;						coordinates. The 1st dimension MUST be 3; if it is
;						not program execution halts.
;						If /to_equatorial is set then the output is in
;						rectangular coordinates.
;
;	/degrees		only used if 'from_equatorial' or /to_equatorial is set
;						if set, all angles are in degrees (default: radians).
;
;	/force_unit 	only used if 'from_equatorial' or 'from_camera' are set
;						forces conversion of vectors to unit vectors
;						For 'from_equatorial' vectors converted to camera
;						coordinates must be unit vectors. By default only
;						vectors deviating from unit vectors by more than the
;						machine precision are explicitly scaled to unit vectors.
;
;	camera=camera	scalar; type: integer; default: 1
;						camera number (1,2 or 3)
;						This determines the fixed S/C-to-camera quaternion for
;						conversions between camera and equatorial frames;
;						and CCD-related parameters (fov size, etc.) for
;						conversions between camera and ccd coordinates
;	quaternion=quaternion
;					array[4]; type: double
;						Coriolis S/C quaternion
;						Only needed for conversions between camera frame
;						and equatorial frame.
;
;	from_mode=from_mode
;	to_mode=to_mode
;	mode=mode		scalar; type: integer or string; default: 0
;
;					If /to_ccd and/or /from_ccd are set:
;
;						mode number (0,1 or 2)
;
;						The 'from_mode' and 'to_mode' keyword are necessary only
;						when both 'from_ccd' and /to_ccd are used (i.e. when
;						converting CCD coordinates between different modes.
;						'from_mode' takes priority over 'mode' if from_ccd is set.
;						'to_mode' takes priority over 'mode' if /to_ccd is set.
;
;                       If /to_map and/or /from_map are set:
;
;                       mode string: "lores", "eqmap", "npmap" or "spmap".
;
;						The 'from_mode' and 'to_mode' keyword are necessary only
;						when both 'from_map' and /to_map are used (i.e. when
;						converting between different types of skymap.
;
;	/boolean		only used for conversions between camera frame and
;					CCD coordinates.
;						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).
; OUTPUTS:
;	R				array[3,n]; type: float
;						output locations as specified by /to_equatorial, /to_camera
;						or /to_ccd keyword.
; OPTIONAL OUTPUT PARAMETER:
;	inside=inside	array[n]; type: long (/boolean NOT set) or byte (/boolean set).
;						Will only exist for conversion between CCD and sky coordinates
;						(either camera or equatorial)..
;						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
; COMMON BLOCKS:
	common cvsmei_save, map_modes, str_modes, cxy_modes, dpb_modes, dim_modes, $
		from_maptype_save, to_maptype_save
; CALLS:
;	ToRadians, InitVar, smei_camera, IsType, SyncDims, boost, smei_cam_quaternion
;	CvRotation
; PROCEDURE:
;	There are four groups of coordinate systems involved:
;	- CCD coordinates (for three different modes)
;	- The native camera coordinate frame
;	- The J2000 equatorial coordinate frame
;	- Skymap coordinates (for lores maps, and three different
;	  hires maps: equator, north and south pole)
;
;	Conversions are done as follows:
;
;		CCD <---> Camera <---> J2000 Equatorial <---> Skymaps
;
;	Phase 1:
;	if CCD  		if /to_ccd then convert to CCD (m0)
;					if not /to_ccd then convert to Camera
;	if Camera		do nothing
;	if J2000		do nothing
;	if Skymaps		convert to J2000
;
;	After phase 1 we have Camera or J2000 coordinates
;	unless /to_ccd was set (then we have CCD m0 coordinates)
;
;	Phase 2:
;	if Camera 		if /to_map then convert to J2000
;	if J2000		if /to_ccd then convert to Camera
;
;	We now have coordinates that are at most one
;	conversion from the final frame.
;
;	Phase 3:
;	if /to_ccd			convert from CCD (m0) or Camera
;	if /to_camera		already in Camera or convert from J2000
;	if /to_equatorial	already in J2000 or convert from Camera
;	if /to_map			convert from J2000

;	Conversions between camera coordinates and equatorial
;	coordinates require a Coriolis quaternion and a camera
;	quaternion. This is provided by explicitly specifying a
;	camera number and Coriolis quaternion through keywords
;	'camera' and 'quaternion'
; MODIFICATION HISTORY:
;	MAR-2005, Paul Hick (UCSD/CASS)
;		Extracted from smei_sky2ccd.
;	MAR-2008, Paul Hick (UCSD/CASS)
;		Bug fix in conversion from hires equat map to polar
;		maps. Force output for RA into [0,360]
;	APR-2008, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;		Fairly substantial rewrite to incorporate conversions
;		from and to helioecliptic maps.
;-

InitVar, to_ccd 	, /key
InitVar, to_camera	, /key
InitVar, to_equatorial	 , /key
InitVar, to_helioecliptic, /key
InitVar, to_map		, /key

InitVar, rectangular, /key
InitVar, boolean	, /key
InitVar, force_unit , /key

InitVar, silent 	, 0

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

; Initialize scaling parameters to defaults used in full-sky SMEI skymaps

IF IsType(map_modes,/undefined) THEN cvsmei_init
mpb_modes = dpb_modes/dpm

InitVar, camera , 1
IF IsType(camera,/string) THEN camera = round(flt_string(camera))

; Keep track of coordinate system

is_map = 0
is_equatorial = 0
is_helioecliptic = 0
is_camera = 0
is_engccd = 0

CASE IsType(from_map, /defined) OF

0: BEGIN

	InitVar, mode, 0

	from_type = -1
	from_map_helioecliptic = 0
	from_map_equatorial    = 0

END

1: BEGIN

	InitVar, mode	  , from_maptype_save
	InitVar, from_mode, mode

	from_type = strpos(map_modes,strlowcase(from_mode))
	from_type = (where(from_type NE -1))[0]
	IF from_type EQ -1 THEN from_type = 1 ; Default: eqmap

	from_map_helioecliptic = from_type EQ 5
	from_map_equatorial	 = 1-from_map_helioecliptic

END

ENDCASE

; If transforming to skymap distinguish
; between equatorial and helioecliptic maps

CASE to_map OF

0: BEGIN
	InitVar, mode, 0

	to_type = -1
	to_map_helioecliptic = 0
	to_map_equatorial	 = 0
END

1: BEGIN

	InitVar, mode	, to_maptype_save
	InitVar, to_mode, mode

	to_type = strpos(map_modes,strlowcase(to_mode))
	to_type = (where(to_type NE -1))[0]
	IF to_type EQ -1 THEN to_type = 1	; Default: eqmap

	to_map_helioecliptic = to_type EQ 5
	to_map_equatorial	 = 1-to_map_helioecliptic

END

ENDCASE

CASE 1 OF

; Convert CCD coordinates to native camera coordinates, unless
; /to_ccd is set, in which case convert to eng mode (mode 0)
; CCD coordinates.

IsType(from_ccd,/defined): BEGIN

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

	InitVar, from_mode, mode
	IF IsType(from_mode,/string) THEN from_mode = round(flt_string(from_mode))

	; If from_ccd is a string set up a 2D grid of x,y positions across
	; the entire CCD for mode 'from_mode'

	IF IsType(from_ccd, /string) THEN $
		from_ccd = gridgen(smei_camera(fov, /get_nsize)/2^from_mode)

	sz = size(from_ccd)
	nd = (size(from_ccd,/dim))[0]		; # components (2)
	nv = n_elements(from_ccd)/nd        ; # verctors

	xy = reform(from_ccd,nd,nv)

	IF from_mode NE 0 THEN	$			; Convert to engineering mode (mode 0)
		xy = -0.5d0+(xy+0.5d0)*(2.0d0^from_mode)

	CASE to_ccd OF

	0: BEGIN

		xy -= smei_camera(fov,/get_center)#replicate(1,nv)
 	
		; Convert to polar coordinates
		; Make sure the phase angle is in [-!pi,+!pi] before calling Inside_Wedge

		xy = cv_coord(from_rect=xy, /to_polar, degrees=degrees)
		xy[0,*] = AngleRange(xy[0,*], /pi, degrees=degrees) ; NOT REDUNDANT

		r = make_array(type=IsType(xy), dim=[3,nv], value=BadValue(xy))

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

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

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

		i_inside = where(inside, n_inside)

		IF n_inside NE 0 THEN BEGIN

			xy = xy[*,i_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.

			xy -= 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

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

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

			r[*,i_inside] = xy

		ENDIF

		IF NOT boolean THEN inside = i_inside

		is_camera = 1

		; r now contains unit vectors in native camera coordinates

	END

	1: BEGIN

		r = xy
		is_engccd = 1

		; r now contains engineering mode (mode 0) CCD pixel coordinates

	END

	ENDCASE

	str_from = 'CCD (m'+strcompress(from_mode,/rem)+')'

END

; Keep native camera coordinates

IsType(from_camera,/defined): BEGIN

	sz = size(from_camera)
	nd = (size(from_camera,/dim))[0]	; # components (3)
	nv = n_elements(from_camera)/nd		; # vectors

	r = reform(from_camera,nd,nv)

	CASE force_unit OF					; Make sure we have unit vectors
	0: BEGIN
		tmp = total(r*r,1)
		i   = where(abs(tmp-1) GT 3*(machar(double=IsType(r,/double))).eps)
		IF i[0] NE -1 THEN r[*,i] /= [1,1,1]#[sqrt(tmp[i])]
	END
	1: r /= [1,1,1]#[sqrt(total(r*r,1))]
	ENDCASE

	is_camera = 1

	str_from = 'native camera frame'

	; r now contains unit vectors for native camera coordinates

END

; Keep J2000 equatorial sky coordinates

IsType(from_equatorial,/defined): BEGIN

	sz =  size(from_equatorial)
	nd = (size(from_equatorial,/dim))[0]; # components (2 or 3)
	nv =  n_elements(from_equatorial)/nd; # vectors

	IF rectangular AND nd NE 3 then message, 'from_equatorial must be dim [3,n]'

	r = reform(from_equatorial,nd,nv)

	is_equatorial = 1
	is_rectangular = rectangular

	str_from = 'equatorial sky'

	; r is now in J2000 equatorial coordinates
	; (spherical or rectangular as set by is_rectangular)

END

IsType(from_helioecliptic,/defined): BEGIN

	sz =  size(from_helioecliptic)
	nd = (size(from_helioecliptic,/dim))[0]; # components (2 or 3)
	nv =  n_elements(from_helioecliptic)/nd; # vectors

	IF rectangular AND nd NE 3 THEN $
		message, 'from_helioecliptic must be dim [3,n]'

	r = reform(from_helioecliptic,nd,nv)

	is_helioecliptic = 1
	is_rectangular = rectangular

	str_from = 'helioecliptic sky'

END

; Convert map coordinates to J2000 equatorial sky coordinates,
; except when converting between hires equatorial and lores
; maps (these are just linear scalings, so no need to go
; to sky coordinates).

IsType(from_map, /defined): BEGIN

	IF IsType(from_map,/string) THEN ndim = from_map

	CASE IsType(ndim,/defined) OF
	0: r = from_map
	1: BEGIN
		IF IsType(ndim,/string) THEN ndim = dim_modes[*,from_type]
		r = gridgen(ndim)
	END
	ENDCASE

	sz =  size(r)
	nd = (size(r,/dim))[0]			; # components (2)
	nv =  n_elements(r)/nd			; # vectors

	r = reform(r,nd,nv,/overwrite)

	; Stay in map when converting between hires
	; equatorial (from_type=1) and lores (from_type=0,
	; always equatorial) map coordinates

	stay_in_map = to_map AND $
		((from_type LE 1 AND to_type LE 1) OR	$
		 (from_type EQ 5 AND to_type EQ 5))

	CASE stay_in_map OF

	0: BEGIN

		r -= cxy_modes[*,from_type]#replicate(1,nv)	; Shift origin to [cx,cy]

		; Convert to sky coordinates:
		; equatorial (from_type NE 5) or helioecliptic (from_type EQ 5)

		mpb = mpb_modes[from_type]

		CASE from_type OF
		0: r *= mpb			; lores 
		1: r *= mpb			; eqmap
		2: BEGIN			; npmap
			r = [atan(r[1,*],r[0,*])/rpm, reform( 90.0d0/dpm-sqrt(total(r*r,1))*mpb,1,nv)]
			r[0,*] = AngleRange(r[0,*],degrees=degrees)
		END
		3: BEGIN			; spmap
			r = [atan(r[1,*],r[0,*])/rpm, reform(-90.0d0/dpm+sqrt(total(r*r,1))*mpb,1,nv)]
			r[0,*] = AngleRange(r[0,*],degrees=degrees)
		END
		4: r *= mpb			; eqfull
		5: r *= mpb			; eclsun
		ENDCASE

		is_rectangular = 0

		IF from_type EQ 5 THEN is_helioecliptic = 1 ELSE is_equatorial = 1

	END

	1: BEGIN

		is_map = 1

	END

	ENDCASE

	str_from = str_modes[from_type]+' map'

	; r is now in spherical J2000 equatorial coordinates

END

ENDCASE

sz_3 = sz
sz_3[1] = 3 				; Three components
sz_3[sz_3[0]+2] = sz_3[sz_3[0]+2]/nd*3

sz_2 = sz
sz_2[1] = 2 				; Two components
sz_2[sz_2[0]+2] = sz_2[sz_2[0]+2]/nd*2


; r is now a set of unit vectors in the camera frame EXCEPT WHEN
; - converting from_ccd -> to_ccd, in which case we now
;   have pixel coordinates in engineering mode)
; - converting from_map, in which case r is now J2000 RA/dec (is_equatorial=1)

IF is_helioecliptic AND (to_camera OR to_ccd OR to_map_equatorial) THEN BEGIN

	IF NOT IsTime(time) THEN message, 'need time for this conversion'

	; what about precession here?

	r = CvSky(time,from_centered_ecliptic=r,/to_equatorial	, $
		degrees=degrees,rectangular=is_rectangular)

	is_helioecliptic = 0
	is_equatorial = 1

ENDIF

IF is_equatorial AND to_ccd THEN BEGIN

	; Convert to native camera coordinates

	IF NOT is_rectangular THEN BEGIN ; Unit vector in spherical coordinates
		IF (size(r,/dim))[0] EQ 2 THEN boost, r, replicate(1, nv)
		r = cv_coord(from_sphere=r, /to_rect, degrees=degrees)
		is_rectangular = 1
	ENDIF							; Convert to rectangular coordinates

	CASE force_unit OF				; Make sure we have unit vectors
	0: BEGIN
		tmp = total(r*r,1)
		i   = where(abs(tmp-1) GT 2*(machar(double=IsType(r,/double))).eps)
		IF i[0] NE -1 THEN r[*,i] /= [1,1,1]#[sqrt(tmp[i])]
	END
	1: r /= [1,1,1]#[sqrt(total(r*r,1))]
	ENDCASE

	; Rotate to the UCSD camera frame

	InitVar, quaternion, [1.0d0,1.0d0,1.0d0,0.0d0]/sqrt(3.0d0)
	dcm = smei_cam_quaternion(quaternion,camera)
	dcm = CvRotation(from_quaternion=dcm, /to_dcm)
	r = transpose(dcm)#r

	is_equatorial = 0
	is_camera = 1

	; r now contains unit vectors for native camera coordinates

ENDIF

IF is_camera AND (to_map OR to_helioecliptic) THEN BEGIN

	; Convert to J2000 equatorial coordinates

	InitVar, quaternion, [1.0d0,1.0d0,1.0d0,0.0d0]/sqrt(3.0d0)
	dcm = smei_cam_quaternion(quaternion,camera)
	dcm[0:2] = -dcm[0:2]		; Invert quaternion
	dcm = CvRotation(from_quaternion=dcm, /to_dcm)
	r = transpose(dcm)#r

	is_camera = 0
	is_equatorial = 1
	is_rectangular = 1

	; r is now in rectangular J2000 equatorial coordinates

ENDIF

IF is_equatorial AND to_map_helioecliptic THEN BEGIN

	IF NOT IsTime(time) THEN message, 'need time for this conversion'

	; what about precession here?
	; range of heliolongitude should be -180,180

	r = CvSky(time,from_equatorial=r,/to_centered_ecliptic, $
		degrees=degrees,rectangular=rectangular)

	is_equatorial = 0
	is_helioecliptic = 1
	
ENDIF


CASE 1 OF

to_map_equatorial: BEGIN

	CASE 1 OF

	is_engccd		: message, 'oops, is_engccd'
	is_camera		: message, 'oops, is_camera'
	is_helioecliptic: message, 'oops, is_helioecliptic'

	is_equatorial: BEGIN

		; Convert to spherical coordinates

		IF is_rectangular THEN BEGIN
			r = cv_coord(from_rect=r, /to_sphere, degrees=degrees)
			is_rectangular = 0
		ENDIF

		; If present, drop the distance (3rd component)

		IF (size(r,/dim))[0] EQ 3 THEN r = r[0:1,*]

		mpb = mpb_modes[to_type]

		CASE to_type OF
		0: r /= mpb			; lores
		1: r /= mpb			; eqmap
		2: r = (90.0/dpm-[r[1,*],r[1,*]])*[cos(r[0,*]*rpm),sin(r[0,*]*rpm)]/mpb ; npmap
		3: r = (90.0/dpm+[r[1,*],r[1,*]])*[cos(r[0,*]*rpm),sin(r[0,*]*rpm)]/mpb ; spmap
		4: r /= mpb			; eqfull
		ENDCASE

		r += cxy_modes[*,to_type]#replicate(1,nv)

		SyncDims, r, size=sz_2

		; r is now in equatorial skymap coordinates

	END

	is_map: BEGIN

		; Only happens when converting from lores sky map
		; to hires equatorial map or v.v

		IF to_type NE from_type THEN BEGIN
			r -= cxy_modes[*,from_type]#replicate(1,nv)
			r *= mpb_modes[from_type]/mpb_modes[to_type]
			r += cxy_modes[*,to_type]#replicate(1,nv)
		ENDIF

	END

	ENDCASE

	count = nv
	str_to = str_modes[to_type]+' map'

END

to_map_helioecliptic: BEGIN

	CASE 1 OF

	is_engccd		: message, 'oops, is_engccd'
	is_camera		: message, 'oops, is_camera'

	is_helioecliptic: BEGIN

		; Convert to spherical coordinates

		IF is_rectangular THEN BEGIN
			r = cv_coord(from_rect=r, /to_sphere, degrees=degrees)
			is_rectangular = 0
		ENDIF

		; If present, drop the distance (3rd component)

		IF (size(r,/dim))[0] EQ 3 THEN r = r[0:1,*]

		mpb = mpb_modes[to_type]

		r /= mpb			; eclsun
		r += cxy_modes[*,to_type]#replicate(1,nv)

		SyncDims, r, size=sz_2

		; r is now in helioecliptic skymap coordinates

	END

	is_equatorial: message, 'oops, is_equatorial'

	is_map: BEGIN

		; Only happens when converting from from_type=5 to
		; to_type=5, i.e. when doing nothing.
		; So this is never executed, I think.

		IF to_type NE from_type THEN BEGIN
			r -= cxy_modes[*,from_type]#replicate(1,nv)
			r *= mpb_modes[from_type]/mpb_modes[to_type]
			r += cxy_modes[*,to_type]#replicate(1,nv)
		ENDIF

	END

	ENDCASE

	count = nv
	str_to = str_modes[to_type]+' map'

END

to_helioecliptic: BEGIN

	CASE 1 OF

	is_engccd		: message, 'oops, is_engccd'
	is_camera		: message, 'oops, is_camera'
	is_helioecliptic:

	is_equatorial	: BEGIN

		IF NOT IsTime(time) THEN message, 'need time for this conversion'

		; what about precession here?
		; range of heliolongitude should be -180,180

		r = CvSky(time,from_equatorial=r,/to_centered_ecliptic, $
			degrees=degrees,rectangular=is_rectangular)

		is_helioecliptic = 1

	END

	is_map			: message, 'oops, is_map'

	ENDCASE

	IF rectangular NE is_rectangular THEN BEGIN

		CASE rectangular OF
		0: BEGIN
			r = (cv_coord(from_rect=r, /to_sphere, degrees=degrees))[0:1,*]
			SyncDims, r, size=sz_2
		END
		1: BEGIN
			IF (size(r,/dim))[0] EQ 2 THEN boost, r, replicate(1, nv)
			r = cv_coord(from_sphere=r, /to_rect, degrees=degrees)
			SyncDims, r, size=sz_3
		END
		ENDCASE

	ENDIF ELSE BEGIN

		CASE rectangular OF
		0: SyncDims, r, size=sz_2
		1: SyncDims, r, size=sz_3
		ENDCASE

	ENDELSE

	count  = nv
	str_to = 'helioecliptic sky'

END

to_equatorial: BEGIN

	CASE 1 OF
	is_engccd: message, 'oops, is_engccd'

	is_camera: BEGIN

		InitVar, quaternion, [1.0d0,1.0d0,1.0d0,0.0d0]/sqrt(3.0d0)
		dcm = smei_cam_quaternion(quaternion,camera)
		dcm[0:2] = -dcm[0:2]		; Invert quaternion
		dcm = CvRotation(from_quaternion=dcm, /to_dcm)
		r = transpose(dcm)#r

		is_camera = 0
		is_equatorial = 1
		is_rectangular = 1

	END

	is_helioecliptic: BEGIN
		IF NOT IsTime(time) THEN message, 'need time for this conversion'

		; what about precession here?

		r = CvSky(time,from_centered_ecliptic=r,/to_equatorial	, $
			degrees=degrees,rectangular=rectangular)

		is_helioecliptic = 0
		is_equatorial = 1
	END

	is_equatorial:
	is_map		 : message, 'oops, is_map'

	ENDCASE

	IF rectangular NE is_rectangular THEN BEGIN

		CASE rectangular OF
		0: BEGIN
			r = (cv_coord(from_rect=r, /to_sphere, degrees=degrees))[0:1,*]
			r[0,*] = AngleRange(r[0,*],degrees=degrees) ; RA in [0,360]
			SyncDims, r, size=sz_2
		END
		1: BEGIN
			IF (size(r,/dim))[0] EQ 2 THEN boost, r, replicate(1, nv)
			r = cv_coord(from_sphere=r, /to_rect, degrees=degrees)
			SyncDims, r, size=sz_3
		END
		ENDCASE

	ENDIF ELSE BEGIN

		CASE rectangular OF
		0: SyncDims, r, size=sz_2
		1: SyncDims, r, size=sz_3
		ENDCASE

	ENDELSE

	count  = nv
	str_to = 'equatorial sky'

	; r is now in J2000 equatorial coordinates (spherical or rectangular
	; depending on setting of 'rectangular')

END

to_camera: BEGIN

	CASE 1 OF

	is_engccd		: message, 'oops, is_engccd'
	is_camera		:
	is_helioecliptic: message, 'oops, is_helioecliptic'

	is_equatorial: BEGIN

		IF NOT is_rectangular THEN BEGIN ; Unit vector in spherical coordinates
			IF (size(r,/dim))[0] EQ 2 THEN boost, r, replicate(1, nv)
			r = cv_coord(from_sphere=r, /to_rect, degrees=degrees)
			is_rectangular = 1
		ENDIF							; Convert to rectangular coordinates

		CASE force_unit OF				; Make sure we have unit vectors
		0: BEGIN
			tmp = total(r*r,1)
			i   = where(abs(tmp-1) GT 2*(machar(double=IsType(r,/double))).eps)
			IF i[0] NE -1 THEN r[*,i] /= [1,1,1]#[sqrt(tmp[i])]
		END
		1: r /= [1,1,1]#[sqrt(total(r*r,1))]
		ENDCASE

		; Rotate to the UCSD camera frame

		InitVar, quaternion, [1.0d0,1.0d0,1.0d0,0.0d0]/sqrt(3.0d0)
		dcm = smei_cam_quaternion(quaternion,camera)
		dcm = CvRotation(from_quaternion=dcm, /to_dcm)
		r = transpose(dcm)#r

		is_camera = 1

	END

	is_map: message, 'oops, is_map'
	
	ENDCASE

	SyncDims, r, size=sz_3

	count  = nv
	str_to = 'native camera frame'

	; r now contains unit vectors in native camera coordinates

END

to_ccd: BEGIN

	count = nv

	CASE 1 OF

	is_engccd:

	is_camera: BEGIN

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

		; Only locations closer then fov.size[0]/2 (= 30 degrees) to the z-axis
		; potentially lie inside the fov.
		; 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??>

		tmp = smei_camera(fov, /get_fov_size)
		tmp = 1.1*(tmp[0,1]-tmp[0,0])/2
		tmp = tmp < (89.9/dpm)				; Safety belt
		tmp = cos(tmp*rpm)

		i_inside = where(reform(r[2,*]) GE tmp, n_inside)

		xy = r[0:1,*]						; Drop the z-component
		r  = make_array(size=size(xy), value=BadValue(xy)); Output array

		IF n_inside NE 0 THEN BEGIN 		; First cut; points above horizon

			xy = xy[*,i_inside]

			; 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

			xy = asin(xy)*(dpm/rpm)
			xy[1,*] = smei_theta2radial(xy[1,*], xy[0,*]); Unit is engineering pixels
			xy[0,*] = xy[0,*]/dpm 			; Convert angle to mystery units

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

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

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

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

			CASE n_inside EQ 0 OF 			; Second cut; inside fov
			0: BEGIN
				i_inside = i_inside[tmp]
				xy = xy[*,tmp]

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

				r[*,i_inside] = smei_camera(fov, /get_center)#replicate(1,n_inside)+	$
					cv_coord(from_polar=xy, /to_rect, degrees=degrees)
			END
			1: i_inside = -1
			ENDCASE

		ENDIF

		CASE boolean OF
		0: inside = i_inside
		1: BEGIN
			inside = SubArray(r,value=0B)
			IF n_inside NE 0 then inside[i_inside] = 1B
		END
		ENDCASE

		count = n_inside

	END

	is_helioecliptic: message, 'oops, is helioecliptic'
	is_equatorial	: message, 'oops, is_equatorial'
	is_map			: message, 'oops, is_map'

	ENDCASE

	SyncDims, r, size=sz_2

	; Convert from eng mode to requested mode

	InitVar, to_mode, mode
	IF IsType(to_mode,/string) THEN to_mode = round(flt_string(to_mode))
	IF to_mode NE 0 THEN r = -0.5d0+(r+0.5d0)/(2.0d0^to_mode)

	str_to = 'CCD (m'+strcompress(to_mode,/rem)+'}'

	; r is now in CCD pixel coordinates for specified mode

END

ENDCASE

IF silent LE 0 THEN $
	message, /info, who_am_i(/caller)+': from "'+str_from+'" to "'+str_to+'"'

RETURN, r  &  END