;+
; NAME:
;	smei_mkstarmask
; PURPOSE:
;	Creates the star mask file used in the creation of
;	skymaps by href=smei_skyd=
; CATEGORY:
;	camera/idl/star
; CALLING SEQUENCE:
	PRO smei_mkstarmask, notest=notest
; INCLUDE:
	@compile_opt.pro		; On error, return to caller
; PROCEDURE:
;	Uses grd files in $SSWDB_SMEI as input, and creates
;	starmask.fts.gz in the same directory.
;	starmask.fts needs to be moved to $SMEIDB
;	where the Fortran routine href=smei_get_starmask= looks for it.
; MODIFICATION HISTORY:
;	DEC-2008, Paul Hick (UCSD/CASS)
;	JUN-2011, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;		Bugfix: maps were stored in wrong coordinate system. The new
;		version converst from Andy's convention used in the grd file
;		to the SMEI skymap convention (the same is done for the
;		sidereal background maps in smei_mksidereal.pro)
;		Bugfix: resolution (bins/degrees) was specified wrong in Fits header.
;		Made Fits header WCS compliant.
;		Added Fits header keyword for star mask version number (now 2.0)
;-

InitVar, notest, /key

root = filepath(root=getenv('SSWDB_SMEI'),'starmask')

eq_mask = filepath(root=root,'EqStarMask_440.grd')
np_mask = filepath(root=root,'NpStarMask_440.grd')
sp_mask = filepath(root=root,'SpStarMask_440.grd')

eq_i = grd_read(eq_mask, eq_map, /norow)
np_i = grd_read(np_mask, np_map, /norow)
sp_i = grd_read(sp_mask, sp_map, /norow)

eq_map = byte(eq_map)
np_map = byte(np_map)
sp_map = byte(sp_map)

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

eq_map = rotate(eq_map, 5)		; x -> -x
np_map = rotate(np_map, 7)		; y -> -y
sp_map = rotate(sp_map, 2)		; x -> -x and y -> -y

sz = size(eq_map,/dim)
cdeq = 360.0d0/sz[0]
;d_eq = 1/cdeq
cxeq = 0.5d0
cyeq = double(sz[1]/2)+0.5d0

print, cxeq, cyeq, cdeq

sz = size(np_map,/dim)
cdpl = 0.1d0
;d_pl = 1/cdpl
cxpl = double(sz[0]/2)+0.5d0
cypl = double(sz[1]/2)+0.5d0

print, cxpl, cypl, cdpl

fts_mask = filepath(root=root,'starmask'+(['_test',''])[notest]+'.fts')

mkhdr, hdr, eq_map, /extend
fxaddpar, hdr, 'SMASKVER', 2.00, 'Star mask version number'
fxaddpar, hdr, 'MAP','Binary star mask Equatorial map RA=[0,360],DEC=[-60,+60]', 'Map description'
;fxaddpar, hdr, 'CX',cxeq, 'Array index (base 1) for RA=0'
;fxaddpar, hdr, 'CY',cyeq, 'Array index (base 1) for DEC=0'
;fxaddpar, hdr, 'BD',d_eq, 'Number of bins per degree'

fxaddpar, hdr, 'RADESYS', 'FK5'  , 'IAU 1984 system'
fxaddpar, hdr, 'EQUINOX', 2000.0 , 'J2000 coordinates'
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' , cxeq   , 'Array index (base 1) for RA=0'
fxaddpar, hdr, 'CRPIX2' , cyeq   , 'Array index (base 1) for DEC=0'
fxaddpar, hdr, 'CDELT1' , cdeq   , 'Degrees per bin'
fxaddpar, hdr, 'CDELT2' , cdeq   , 'Degrees per bin'
fxaddpar, hdr, 'CRVAL1' , 0.00   , 'RA=0.0'
fxaddpar, hdr, 'CRVAL2' , 0.00   , 'DEC=0.0'

writefits, fts_mask, eq_map, hdr

mkhdr, hdr, np_map, /image
fxaddpar, hdr, 'MAP','Binary star mask map of north pole DEC=[+50,+90]', 'Map description'
;fxaddpar, hdr, 'CX', cxpl, 'Array index (base 1) for north pole'
;fxaddpar, hdr, 'CY', cypl, 'Array index (base 1) for north pole'
;fxaddpar, hdr, 'BD', d_pl, 'Number of bins per degree'

fxaddpar, hdr, 'RADESYS' , 'FK5'     , 'IAU 1984 system'
fxaddpar, hdr, 'EQUINOX' , 2000.0    , 'J2000 coordinates'
fxaddpar, hdr, 'CTYPE1'  , 'RA---ARC', 'Axis type: equidistant azimuthal'
fxaddpar, hdr, 'CTYPE2'  , 'DEC--ARC', 'Axis type: equidistant azimuthal'
fxaddpar, hdr, 'CUNIT1'  , 'deg'     , 'Angles in degrees'
fxaddpar, hdr, 'CUNIT2'  , 'deg'     , 'Angles in degrees'
fxaddpar, hdr, 'CRPIX1'  , cxpl      , 'Array index (base 1) for north pole'
fxaddpar, hdr, 'CRPIX2'  , cypl      , 'Array index (base 1) for north pole'
fxaddpar, hdr, 'CDELT1'  , cdpl      , 'Degrees per bin'
fxaddpar, hdr, 'CDELT2'  , cdpl      , 'Degrees per bin'
fxaddpar, hdr, 'CRVAL1'  , 90.00     , 'RA of north pole'
fxaddpar, hdr, 'CRVAL2'  , 90.00     , 'DEC of north pole'
fxaddpar, hdr, 'LONPOLE' ,  0.00     , 'Native longitude of celestial pole'

writefits, fts_mask, np_map, hdr, /append

mkhdr, hdr, sp_map, /image
fxaddpar, hdr, 'MAP','Binary star mask map of south pole DEC=[-50,-90]', 'Map description'
;fxaddpar, hdr, 'CX', cxpl, 'Array index (base 1) for south pole'
;fxaddpar, hdr, 'CY', cypl, 'Array index (base 1) for south pole'
;fxaddpar, hdr, 'BD', d_pl, 'Number of bins per degree'

fxaddpar, hdr, 'RADESYS', 'FK5'     , 'IAU 1984 system'
fxaddpar, hdr, 'EQUINOX', 2000.0    , 'J2000 coordinates'
fxaddpar, hdr, 'CTYPE1' , 'RA---ARC', 'Axis type: equidistant azimuthal'
fxaddpar, hdr, 'CTYPE2' , 'DEC--ARC', 'Axis type: equidistant azimuthal'
fxaddpar, hdr, 'CUNIT1' , 'deg'     , 'Angles in degrees'
fxaddpar, hdr, 'CUNIT2' , 'deg'     , 'Angles in degrees'
fxaddpar, hdr, 'CRPIX1' ,  cxpl     , 'Array index (base 1) for south pole'
fxaddpar, hdr, 'CRPIX2' ,  cypl     , 'Array index (base 1) for south pole'
fxaddpar, hdr, 'CDELT1' , -cdpl     , 'Degrees per bin'
fxaddpar, hdr, 'CDELT2' ,  cdpl     , 'Degrees per bin'
fxaddpar, hdr, 'CRVAL1' , -90.00    , 'RA of south pole'
fxaddpar, hdr, 'CRVAL2' , -90.00    , 'DEC of south pole'
fxaddpar, hdr, 'LONPOLE',   0.00    , 'Native longitude of celestial pole'

writefits, fts_mask, sp_map, hdr, /append

spawn, 'gzip -vf '+fts_mask

RETURN  &  END