;+
; NAME:
;	smei_star_box
; PURPOSE:
;	Set up box around star centroid of about the same
;	size as the standard star map.
; CATEGORY:
; 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
;	maptype 			scalar; type: integer
;	std_ndim			array[2]; type: integer
;	ndim				array[2]; type: integer
;	xystar				array[2]; type: double
; OUTPUTS:
;	result				array[m]; type: integer
;	ibeg				scalar; type: integer
;	iend				scalar; type: integer
;	jbeg				scalar; type: integer
;	jend				scalar; type: integer
;	rr					array[2,n]; type: double
; OPTIONAL OUTPUT PARAMETERS:
;	ibeg0=ibeg0			scalar; type: integer
;	jbeg0=jbeg0			scalar; type: integer
; INCLUDE:
	@compile_opt.pro		; On error, return to caller
; CALLS:
;	gridgen
; PROCEDURE:
; MODIFICATION HISTORY:
;	JUL-2007, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;-

; 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 = jbeg > 0 			; Smallest yy
	jend = 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 sky map 'final_sky'
; ipix is a monotonic sequence of x-values unless the star bridges
; RA=0=360, in which case the ipix is adjusted to wrap around to
; the opposite side of the sky map.

ipix = xx[0:std_ndim[0]-1] 		; Full range of xx values

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

		; For polar maps, 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 = ibeg > 0 		; Smallest xx
		iend = iend < (ndim[0]-1); Largest xx

		ipix = reform(rr[0,0:iend-ibeg])

	END

	1: BEGIN					; Equatorial map

		; For equatorial maps with the RA outside the map wrap
		; around 360 degrees to fill up the box.
		; 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