;+
; NAME:
;	smei_sky_getmask
; PURPOSE:
;	Get SMEI mask at specified time and sky location(s)
; CATEGORY:
;	camera/idl/sky
; CALLING SEQUENCE:
	FUNCTION smei_sky_getmask, ut, ra, dec, $
		equator	= equator	, $
		degrees	= degrees	, $
		scalar	= scalar	, $
		silent	= silent
; INPUTS:
;	ut		array[1]; type: time structure or valid time string
;	ra		array[n]; type: float
;				right ascensions or ecliptic longitudes
;	dec		array[m]; type: float
;				equatorial declinations or ecliptic latitudes
; OPTIONAL INPUT PARAMETERS:
;	/equator	if set, ra,dec are intepreted as equatorial coordinates
;				(default: ecliptic coordinates)
;	/degrees	if set, all angles are in degrees (default: radians)
;	/scalar		if set, a single position is returned as a scalar
;				instead of a 1-element array
;	silent=silent	scalar; type: integer
;				larger numbers suppress more messages.
; OUTPUTS:
;	R		array[n.m]; type: float
;				mask values
; INCLUDE:
	@compile_opt.pro		; On error, return to caller
; CALLS:
;	InitVar, TimeUnit, TimeGet, smei_filepath, smei_filename, hide_env
;	IsType, gridgen, CvSky, AngleRange, cvsmei
; PROCEDURE:
; MODIFICATION HISTORY:
;	APR-2009, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;-

InitVar, equator, /key
InitVar, degrees, /key
InitVar, scalar , /key
InitVar, silent , 0

mode = 'msk'
uday = TimeUnit(/day)
day  = TimeGet(ut, /bot, uday)

; Contruct the name of the mask file for specified day and read daily mask

mask_file = filepath(root=smei_filepath(mode=mode),smei_filename(day,camera=0,mode=mode,upto=uday))
IF silent LE 0 THEN message, /info, hide_env(mask_file)
mask = readfits(mask_file,/silent)

IF IsType(mask,/array) THEN BEGIN

	rr = gridgen(/merge, ra, dec)

	; Convert to J2000 equatorial coordinates

	IF NOT equator THEN rr = CvSky(ut, from_ecliptic=rr, /to_equator, degrees=degrees, silent=silent)

	; Convert from J2000 equatorial coordinates to low-res array indices

	rr[0,*] = AngleRange(rr[0,*],degrees=degrees)
	rr = cvsmei(from_equatorial=rr, /to_map, mode='lores', degrees=degrees, silent=silent)

	ndim = size(mask,/dim)

	; The IDL function interpolate does not work outside the array limits
	; (R[0,*] < 0 or > ndim[0]-1, R[1,*] < 0 or > ndim[1]-1)
	; However, the points in R can have values of up to 0.5 outside this
	; range. So add one row/column along the edge of the array.

	; Add 1st row at bottom and last row at top after shifting 180 degree

	mask = [shift(mask[*,0],ndim[0]/2), reform(mask,ndim[0]*ndim[1]), shift(mask[*,ndim[1]-1],ndim[0]/2)]
	mask = reform(mask,ndim[0],ndim[1]+2)

	; Add last column before 1st column and 1st column after last column

	mask = [mask[ndim[0]-1,*],mask,mask[0,*]]

	; Now interpolate. Note that 1 needs to be added to the interpolation points

	mask = interpolate( mask, rr[0,*]+1, rr[1,*]+1, missing=0 )

	mask = reform(mask,n_elements(ra),n_elements(dec),/overwrite)
	mask = reform(mask,/overwrite)

	IF scalar THEN IF n_elements(mask) EQ 1 THEN mask = mask[0]

ENDIF

RETURN, mask  &  END