FUNCTION smei_sky_getmask, ut, ra, dec, $
		equator	= equator	, $
		degrees	= degrees	, $
		scalar	= scalar	, $
		silent	= silent

	@compile_opt.pro

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

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

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)

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

	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