;+
; NAME:
;	smei_mksidereal
; PURPOSE:
;	Combine three grid files:
;		eq_allskymeed.grd.gz
;		np_allskymeed.grd.gz, and
;		sp_allskymeed.grd.gz
;	into one compressed Fits file:
;		$SMEIDB/sid/c0sid_2002_365_000000.fts.gz.
;	The resulting FITS file represents sidereal background information
;	useful for star removal algorithms.
; CALLING SEQUENCE:
	PRO smei_mksidereal, version
; INPUTS:
;	None
; OPTIONAL INPUT PARAMETERS:
;	None
; INCLUDE:
	@compile_opt.pro
; CALLS:
;	grd_read, mkhdr, fxaddpar, writefits, gzip_file, TimeGet
;	TimeSet, TimeUnit, hide_env, IsType, smei_filepath
; RESTRICTIONS:
;	The input grid file paths and output FITS file path are hardcoded
;	into the program.
; PROCEDURE:
;	The three grid files mentioned above are assumed to be valid grid
;	files. Additional FITS header information is added to each map as
;	they are written to the output FITS file.
;	A special version number header keyword, BKGNDVER, is added to the
;	main FITS header.
; MODIFICATION HISTORY:
;	JUL-2006, Jordan Vaughan (UCSD)
;	DEC-2007, Paul Hick (UCSD/CASS)
;		Introduced WCS-compliant keywords to describe coordinate system.
;	APR-2008, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;		Added new skymap version (2.0; not quite final yet).
;-        

InitVar, version, 1.0
version = float(version)

mode = 'sid'

tgain   = TimeGet(/_ydoy,TimeSet(yr=2003,doy=0),upto=TimeUnit(/sec),/scalar)
name    = 'c0'+mode+'_'+tgain

; We will find the background skymaps in the directory $SMEIDB/sidereal.

smeidb = smei_filepath(mode=mode)

rawmapdir = filepath(root=smeidb, 'raw')
skymapdir = filepath(root=smeidb, 'v'+strcompress(round(version),/rem))
IF NOT checkdir(skymapdir) THEN spawn, 'mkdir -v '+skymapdir

skymapfts = filepath(root=skymapdir,name+'.fts')

; Make sure each grid file exists. If they exist, read them.
; If not, display an error message and exit.

CASE version OF
1: marker = 'allskymeed'
2: marker = 'allskymeed2'
3: marker = 'allskymeed_avg12'
ENDCASE

tmp = filepath(root=rawmapdir, 'eq_'+marker+'.grd.gz')

message, /info, hide_env(tmp)
IF NOT grd_read(tmp, eqmap, /silent) THEN BEGIN
	message, /info, 'ERROR: File does not exist or could not be read!'
	RETURN
ENDIF

tmp = filepath(root=rawmapdir, 'np_'+marker+'.grd.gz')
message, /info, hide_env(tmp)
IF NOT grd_read(tmp, npmap, /silent) THEN BEGIN
	message, /info, 'ERROR: File does not exist or could not be read!'
	RETURN
ENDIF

tmp = filepath(root=rawmapdir, 'sp_'+marker+'.grd.gz')
message, /info, hide_env(tmp)
IF NOT grd_read(tmp, spmap, /silent) THEN BEGIN
	message, /info, 'ERROR: File does not exist or could not be read!'
	RETURN
ENDIF

; Do rotations on the maps to get the correct orientations.

eqmap = rotate(eqmap, 5)			; x -> -x
npmap = rotate(npmap, 7)			; y -> -y
spmap = rotate(spmap, 2)			; x -> -x and y -> -y

; Now, create FITS headers for the maps.
; For eqmap
message, /info, 'Creating main FITS header'
mkhdr, mainheader, IsType(eqmap), size(eqmap,/dim), /extend
fxaddpar, mainheader, 'NAME'	,	    name, ' File name'
fxaddpar, mainheader, 'IMG_TYPE',   'skymap_'+mode, ' Image type'
fxaddpar, mainheader, 'BKGNDVER',	 version, ' Sidereal background version number', format='f8.2'
fxaddpar, mainheader, 'TIME'	,	 tgain 	, ' Ref time for gain changes'
fxaddpar, mainheader, 'STIME'	,	 tgain 	, ' Ref time for gain changes'
fxaddpar, mainheader, 'ETIME'	,	 tgain 	, ' Ref time for gain changes'
fxaddpar, mainheader, 'MAP' 	, 'Equatorial map RA=[0,360],DEC=[-60,+60]', ' Map description'
fxaddpar, mainheader, 'RADESYS'	, 	   'FK5', ' IAU 1984 system'
fxaddpar, mainheader, 'EQUINOX'	,     2000.0, ' J2000 coordinates'				, format='f8.1'
fxaddpar, mainheader, 'CTYPE1'	, 		'RA', ' Axis type'
fxaddpar, mainheader, 'CTYPE2'	,	   'DEC', ' Axis type'
fxaddpar, mainheader, 'CUNIT1'	,	   'deg', ' Angles in degrees'
fxaddpar, mainheader, 'CUNIT2'	,	   'deg', ' Angles in degrees'
fxaddpar, mainheader, 'CRPIX1'	,	    0.50, ' Array index (base 1) for RA=0'	, format='f8.2'
fxaddpar, mainheader, 'CRPIX2'	,	  600.50, ' Array index (base 1) for DEC=0' , format='f8.2'
fxaddpar, mainheader, 'CDELT1'	,	    0.10, ' Degrees per bin'				, format='f8.2'
fxaddpar, mainheader, 'CDELT2'	,	    0.10, ' Degrees per bin'				, format='f8.2'
fxaddpar, mainheader, 'CRVAL1'	,		0.00, ' RA=0.0'							, format='f8.2'
fxaddpar, mainheader, 'CRVAL2'  ,		0.00, ' DEC=0.0'						, format='f8.2'
fxaddpar, mainheader, 'BAD_DATA',       0.00, ' Missing data flag'				, format='f8.2'

; For npmap
message, /info, 'Creating extension 1 FITS header'
mkhdr, npheader, IsType(npmap), size(npmap,/dim), /image
fxaddpar, npheader, 'MAP'		, 'Map of north pole DEC=[+50,+90]', ' Map description'
fxaddpar, npheader, 'RADESYS'	,      'FK5', ' IAU 1984 system'
fxaddpar, npheader, 'EQUINOX'	,     2000.0, ' J2000 coordinates'					, format='f8.1'
fxaddpar, npheader, 'CTYPE1'	, 'RA---ARC', ' Axis type: equidistant azimuthal'
fxaddpar, npheader, 'CTYPE2'	, 'DEC--ARC', ' Axis type: equidistant azimuthal'
fxaddpar, npheader, 'CUNIT1'	,	   'deg', ' Angles in degrees'
fxaddpar, npheader, 'CUNIT2'	,	   'deg', ' Angles in degrees' 
fxaddpar, npheader, 'CRPIX1'	,     400.50, ' Array index (base 1) for north pole', format='f8.2'
fxaddpar, npheader, 'CRPIX2'	,     400.50, ' Array index (base 1) for north pole', format='f8.2'
fxaddpar, npheader, 'CDELT1'	,       0.10, ' Degrees per bin'					, format='f8.2'
fxaddpar, npheader, 'CDELT2'	,       0.10, ' Degrees per bin'					, format='f8.2'
fxaddpar, npheader, 'CRVAL1'	,      90.00, ' RA of north pole'					, format='f8.2'
fxaddpar, npheader, 'CRVAL2'	,      90.00, ' DEC of north pole'					, format='f8.2'
fxaddpar, npheader, 'LONPOLE'	,       0.00, ' Native longitude of celestial pole'	, format='f8.2'
fxaddpar, npheader, 'BAD_DATA'	,       0.00, ' Missing data flag'					, format='f8.2'

; For spmap
message, /info, 'Creating extension 2 FITS header'
mkhdr, spheader, IsType(spmap), size(spmap,/dim), /image
fxaddpar, spheader, 'MAP'		, 'Map of south pole DEC=[-50,-90]', ' Map description'
fxaddpar, spheader, 'RADESYS'	,      'FK5', ' IAU 1984 system'
fxaddpar, spheader, 'EQUINOX'	,     2000.0, ' J2000 coordinates'					, format='f8.1'
fxaddpar, spheader, 'CTYPE1'	, 'RA---ARC', ' Axis type: equidistant azimuthal'
fxaddpar, spheader, 'CTYPE2'	, 'DEC--ARC', ' Axis type: equidistant azimuthal'
fxaddpar, spheader, 'CUNIT1'	,	   'deg', ' Angles in degrees'
fxaddpar, spheader, 'CUNIT2'	,	   'deg', ' Angles in degrees' 
fxaddpar, spheader, 'CRPIX1'	,     400.50, ' Array index (base 1) for south pole', format='f8.2'
fxaddpar, spheader, 'CRPIX2'	,     400.50, ' Array index (base 1) for south pole', format='f8.2'
fxaddpar, spheader, 'CDELT1'	,      -0.10, ' Degrees per bin'					, format='f8.2'
fxaddpar, spheader, 'CDELT2'	,       0.10, ' Degrees per bin'					, format='f8.2'
fxaddpar, spheader, 'CRVAL1'	,     -90.00, ' RA of north pole'					, format='f8.2'
fxaddpar, spheader, 'CRVAL2'	,     -90.00, ' DEC of north pole'					, format='f8.2'
fxaddpar, spheader, 'LONPOLE'	,       0.00, ' Native longitude of celestial pole'	, format='f8.2'
fxaddpar, spheader, 'BAD_DATA'	,       0.00, ' Missing data flag'					, format='f8.2'

; Now, write the data to the combined FITS file.
message, /info, 'write '+hide_env(skymapfts)

writefits, skymapfts, eqmap, mainheader
writefits, skymapfts, npmap, npheader, /append
writefits, skymapfts, spmap, spheader, /append

tempvar = gzip_file(skymapfts, /force)		; GZIP the FITS file.

; Success!
RETURN  &  END