;+
; NAME:
;	smei_star_box
; PURPOSE:
;	Set up box around star centroid of about the same
;	size as the standard star map.
; CATEGORY:
;	camera/idl/star
; CALLING SEQUENCE:
	FUNCTION smei_star_box, $
		fullmap	, $
		maptype , $
		std_ndim, $
		ndim	, $
		xystar	, $
		ibeg	, $
		iend	, $
		jbeg	, $
		jend	, $
		rr		, $
		ibeg0=ibeg0, $
		jbeg0=jbeg0
; INPUTS:
;	fullmap 			scalar; type: integer
;							1 if array size 'ndim' is for a full skymap
;							0 if not.
;							Matters only for equatorial maps (maptype=0). For
;							full maps, if the box around the stars partly falls
;							outside the map (RA < 0 or RA > 360) then the box
;							can be filled using the opposite side of the map.
;	maptype 			scalar; type: integer
;							-1: south-polar map
;							 0: equatorial map
;							 1: north-polar map
;							Only used to test for equatorial map (for full maps,
;							see description of 'fullmap' argument.
;	std_ndim			array[2]; type: integer
;							size of the standard star map (50x50)
;	ndim				array[2]; type: integer
;							size of the skymap
;	xystar				array[2]; type: double
;							array indices into the skymap (of size 'ndim')
;							of the star centroid
; OUTPUTS:
;	result				array[m], m <= std_ndim[0]; type: integer
;							see PROCEDURE
;	ibeg				scalar; type: integer
;	jbeg				scalar; type: integer
;							horizontal and vertical index values in the
;							returned 'rr' array for the lower-left corner
;	iend				scalar; type: integer
;	jend				scalar; type: integer
;							same for the upper-right corner
;	rr					array[2,n], n <= std_ndim[0]*std_ndim[1]; type: double
;							see PROCEDURE
; OPTIONAL OUTPUT PARAMETERS:
;	ibeg0=ibeg0			scalar; type: integer
;	jbeg0=jbeg0			scalar; type: integer
;							horizontal and vertical index values in 'rr' for
;							the lower-left corner before rows/columns
;							outside the map were truncated.
;							ibeg0=ibeg and jbeg0=jbeg if the lower-left corner
;							of the 'std_ndim' sized box lies inside the map;
;							if not then at least one of ibeg0 or jbeg0 will
;							be an invalid negative index into the skymap.
; INCLUDE:
	@compile_opt.pro		; On error, return to caller
; CALLS:
;	gridgen
; PROCEDURE:
;	rr[2,n] provides horizontal (xx[0,*]) and vertical (xx[1,*]) indices into
;	a skymap of dimensions 'ndim' for a square box of skylocations centered on
;	the star centroid 'xystar'. The size of the box is usually 'std_ndim', but is
;	smaller if part of the box falls outside the skymap.
;
;	result[n] is an array of horizontal indices into the skymap across the box.
;
;	If the whole box of size 'std_ndim' lies inside the skymap then
;	rr[2,n] gives locations for a box of size `std_ndim`; result[n] is of size
;	std_ndim[0]: result = reform(rr[0,std_ndim[0]) (i.e. a monotonic increasing
;	sequence for the horizontal locations of the bottom row of the box).
;
;	If part the box of size 'std_ndim' sticks out 'above' or 'below' of the map,
;	then the rows outside the map are removed and only coordinates of rows
;	inside the map are returned in rr[2,n]; n=nx x ny with ny < std_ndim[1].
;
;	If part of the box of size 'std_ndim' sticks out 'left' or 'right' of the map
;	then, for polar maps (maptype=+/-1) and for incomplete (fullmap=0) equatorial
;	(maptype=0) maps, the columns outside the map are remove, and only coordinates
;	for columns inside the map are returned in rr[2,n]; n = nx x ny with nx < std_ndim[0].
;	In this case it is still true result[nx] = rr[0:nx-1] (i.e. still a monotonic
;	increasing sequence of horizontal locations of the bottom row of the box.
;
;	For a complete (fullmap=1) equatorial (maptype=0) map, instead of truncating
;	columns from the box of size 'std_ndim' that extend below RA=0 or above RA=360,
;	the columns 'outside' the map are filled in by picking up colmns from the
;	opposite side of the map. In this case rr[2,n]; n=nx x ny with nx=std_ndim[0], and
;	result[nx] with nx=std_ndim[0]. rr[2,n] will contain invalid horizontal
;	indices for the columns outside the map. However, the locations in result[nx]
;	have been adjusted to reflect the 'wrap'-operation, and will contain valid
;	horizontal indices only. This 'result' array should be used by caller to index
;	the skymaps.
; MODIFICATION HISTORY:
;	JUL-2007, Paul Hick (UCSD/CASS)
;	FEB-2010, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;		Added documentation
;-

; Pick up coordinates for part of skymap around star as array
; indices into original sky map. Note that the center of the
; std_ndim sized box is set to xystar, the centroid of a star.
; Pick up a square about the same size as the standard star array.

rr = round(gridgen(std_ndim,origin=0.5*(std_ndim-1)-xystar))

; We know that the centroid xystar is inside the map,
; but part of the std_ndim sized box may still be outside.
; Deal with part of the box outside map.

; The box now is std_ndim[0] x std_ndim[1] pixels

yy = reform(rr[1,*])			; Row indices

jbeg  = yy[0]					; Smallest yy
jend  = yy[n_elements(yy)-1]	; Largest yy
jbeg0 = jbeg
	
IF jbeg LT 0 OR jend GE ndim[1] THEN BEGIN

	; Some of the rows are outside the map.
	; In the vertical direction for both polar and equatorial maps
	; exclude the part of the box outside the map.
	; This excludes entire rows of std_ndim[0] pixels from the box.

	rr   = rr[*,where(0 LE yy AND yy LT ndim[1])]
	jbeg >= 0 			; Smallest yy
	jend <= ndim[1]-1	; Largest yy

ENDIF

; The box now is std_ndim[0] x (jend-jbeg+1) pixels
; where jend-jbeg+1 could be < std_ndim[1]

xx = reform(rr[0,*])			; Column indices

ibeg  = xx[0]					; Smallest xx
iend  = xx[n_elements(xx)-1]	; Largest xx
ibeg0 = ibeg

; ipix is used to access the skymap 'final_sky'
; ipix is a monotonic sequence of x-values unless the star bridges
; RA=0=360, in which case ipix is adjusted to wrap around to
; the opposite side of the skymap.

ipix = xx[0:std_ndim[0]-1] 		; Full range of xx values (from bottom row of box)

IF ibeg LT 0 OR iend GE ndim[0] THEN BEGIN

	; Some of the columns are outside the map.

	CASE maptype EQ 0 AND fullmap OF

	0: BEGIN					; North or south polar map, or
								; a partial equatorial map

		; For a polar map, and for an incomplete equatorial map
		; (that may not cover 360 deg in RA) exclude part of box
		; outside map 
		; This will exclude entire columns of (jend-jbeg+1) pixels.

		rr = rr[*,where(0 LE xx AND xx LT ndim[0])]
		ibeg >= 0 				; Smallest xx
		iend <= ndim[0]-1		; Largest xx

		ipix = reform(rr[0,0:iend-ibeg])
								; Reduced range of xx values (from bottom row of box)
	END

	1: BEGIN					; Equatorial map

		; For a full equatorial map (covering 360 deg in RA)
		; with the RA outside the map (< 0 or > 360) fill up the
		; box by wrapping around 360 degrees and grabbing indices
		; from the opposite RA side of the map.

		; Only array 'ipix' (used to access the sky map arrays
		; 'final_sky' and 'stars_sky') is adjusted.

		tmp = where(ipix LT 0)
		IF tmp[0] NE -1 THEN ipix[tmp] += ndim[0]

		tmp = where(ipix GE ndim[0])
		IF tmp[0] NE -1 THEN ipix[tmp] -= ndim[0]

	END

	ENDCASE

ENDIF

RETURN, ipix  &  END