;+
; NAME:
;	smei_mkmask
; PURPOSE:
;	Produces a daily "sky mask" indicating which parts of the sky were seen
;	by SMEI on a given day (1=always seen; 0=never seen; between 0 and 1:
;	seen in some but not all skymaps).
; CALLING SEQUENCE:
	PRO smei_mkmask, range, $
		nbin		= nbin		, $
		to_smeidb	= to_smeidb	, $
		norbit		= norbit	, $
		remote		= remote	, $
		silent		= silent
; INPUTS:
;	range			array[2]; type: time structure or valid time strings
;						default: ['2006_033','2008_318']
;						time range (start and end day)
; OPTIONAL INPUTS:
;	/to_smeidb		if set then write Fits file to standard
;						directory $SMEIDB/msk
;	norbit=norbit	scalar; type: integer; default: 15
;						number of consecutive orbits to be used for
;						the creation of a mask (the default 15 covers
;						about one day).
;	nbin=nbin		scalar; type: integer; default: 5
;						Sets the resolution of the mask
;						Must be a factor of 1800.
; INCLUDE:
	@compile_opt.pro
; CALLS:
;	InitVar, TimeSet, gridgen, TimeUnit, CheckDir, IsType
;	smei_orbit_get, smei_sky, TimeGet, hide_env
; PROCEDURE:
;	The size of the mask array is [3600,1800]/nbin. For nbin=5 (default)
;	this gives a 720x360 mask with a resolution of 0.5 degrees,
;	covering [0,360] in RA and [-90,+90] in declination.
;
;	The mask value for each skybin is a value between 0 and 1.
;	1 means that the skybin has a value in each of the morbit orbits
;	tested (where morbit <= norbit, is the number of skymaps found);
;	0 means that the skybin did not contain anything in all skymaps.
; MODIFICATION HISTORY:
;	MAR-2009, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;-

InitVar, to_smeidb, /key	; Write to $SMEIDB/msk
InitVar, norbit, 15			; Number of orbits to be used for mask creation
InitVar, nbin, 5			; Makes low-res (720 x 360) map
InitVar, silent, 0

InitVar, range, ['2006_033','2008_318']
tt = TimeSet(range)
IF n_elements(tt) EQ 1 THEN tt = replicate(tt,2)

tt = gridgen(unit=TimeUnit(/day), range=tt)	; Array of times (0 UT)
nt = n_elements(tt)			; Number of times

CASE to_smeidb OF			; Set destination for Fits files
0: InitVar, destination, getenv('TUB')
1: destination = smei_filepath(mode='msk')
ENDCASE

IF IsType(remote,/defined) THEN BEGIN
	remote_destination = remote+':'+destination
	destination = getenv('TUB')
ENDIF

IF NOT CheckDir(destination) THEN	$	; Make sure destination exists
	RETURN

version = 0.00

nx = 3600/nbin				; Resolution of mask
ny = 1800/nbin
mask = fltarr(nx,ny)		; Mask array: 0-360 in RA; -90,90 in Dec

; Set up the main Fits header (should be WCS-compliant).

IF silent LE 0 THEN message, /info, 'Creating main FITS header'
mkhdr, hdr, IsType(mask), size(mask,/dim);, /extend
fxaddpar, hdr, 'NAME'    ,         '', ' File name'
fxaddpar, hdr, 'IMG_TYPE',   'skymap_msk', ' Image type'
fxaddpar, hdr, 'TIME'    ,         '', ' Ref time for gain changes'
fxaddpar, hdr, 'NORBIT'	 ,          0, ' Number of contributing orbits'
fxaddpar, hdr, 'MAP'     , 'Equatorial mask RA=[0,360],DEC=[-90,+90]', ' Map description'
fxaddpar, hdr, 'MASK_VER',    version, ' Version number'				 , format='f8.2'
fxaddpar, hdr, 'RADESYS' ,      'FK5', ' IAU 1984 system'
fxaddpar, hdr, 'EQUINOX' ,     2000.0, ' J2000 coordinates'              , format='f8.1'
fxaddpar, hdr, 'CTYPE1'  ,       'RA', ' Axis type'
fxaddpar, hdr, 'CTYPE2'  ,      'DEC', ' Axis type'
fxaddpar, hdr, 'CUNIT1'  ,      'deg', ' Angles in degrees'
fxaddpar, hdr, 'CUNIT2'  ,      'deg', ' Angles in degrees'
fxaddpar, hdr, 'CRPIX1'  ,       0.50, ' Array index (base 1) for RA=0'  , format='f8.2'
fxaddpar, hdr, 'CRPIX2'  , (ny+1)/2.0, ' Array index (base 1) for DEC=0' , format='f8.2'
fxaddpar, hdr, 'CDELT1'  ,   360.0/nx, ' Degrees per bin'                , format='f8.2'
fxaddpar, hdr, 'CDELT2'  ,   180.0/ny, ' Degrees per bin'                , format='f8.2'
fxaddpar, hdr, 'CRVAL1'  ,       0.00, ' RA=0.0'                         , format='f8.2'
fxaddpar, hdr, 'CRVAL2'  ,       0.00, ' DEC=0.0'                        , format='f8.2'
fxaddpar, hdr, 'BAD_DATA',       0.00, ' Missing data flag'              , format='f8.2'

FOR i=0L,nt-1 DO BEGIN						; Loop over all days

	orbit = TimeGet(/smei,tt[i])			; Orbit structure
	orbit = smei_orbit_get(orbit,/number)	; Orbit that contains start of day

	mask *= 0.0								; Zero out the mask array
	morbit = 0
	FOR j=orbit,orbit+norbit-1 DO BEGIN
		smei_sky, /noplot, j, skymap=skymap, mode='sky', /exact, nbin=nbin, silent=silent
		IF IsType(skymap,/array) THEN BEGIN
			mask += finite(skymap)
			morbit++						; Count contributing orbits
		ENDIF
	ENDFOR

	IF morbit GT 0 THEN BEGIN				; At least one orbit contributing
		mask /= morbit

		ut_time = TimeGet(/_ydoy,tt[i],upto=TimeUnit(/day),/scalar)
		ut_file = 'c0msk_'+ut_time

		fxaddpar, hdr, 'NAME'  , ut_file
		fxaddpar, hdr, 'TIME'  , ut_time
		fxaddpar, hdr, 'NORBIT', morbit

		ut_file = filepath(root=destination,ut_file+'.fts')
		writefits, ut_file, mask, hdr, /compress
		ut_file += '.gz'
		IF silent LE 1 THEN message, /info, 'write '+hide_env(ut_file)

		IF IsType(remote_destination,/defined) THEN	$
			spawn, 'scp '+ut_file+' '+remote_destination+'; rm'+([' ',' -v '])[silent LE 0]+ut_file
	ENDIF

ENDFOR

RETURN  &  END