;+
; NAME:
;	smei_mkzldresidue
; PURPOSE:
;	Stores Andy's zld correction map (1 static and 52 weekly
;	corrections) in a single Fits file
; CATEGORY:
;	camera/idl/zld
; CALLING SEQUENCE:
	PRO smei_mkzldresidue
; INPUTS:
;	(none)
; OUTPUTS:
;	c3zldresidue.fts.gz
; OPTIONAL OUTPUT PARAMETERS:
; INCLUDE:
	@compile_opt.pro		; On error, return to caller
; CALLS:
;	grd_read, GetFileSpec, IsType, hide_env
; PROCEDURE:
;	The Fits headers should be WCS-compliant.
;	The coordinate system used is helioecliptic J2000.
;	Each of Andys maps is 360 x 360 with 0.5 degree resolution.
;	covering [-90,+90] in helioecliptic longitude. This is expanded
;	with zeroes to [-180,+180] making a standard helioecliptic
;	low-res map.
; MODIFICATION HISTORY:
;	APR-2009, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;-

destin = getenv('SMEIDB')
source = filepath(root=destin,'zld')

fts_residue = filepath(root=destin, 'c3zldresidue.fts')
fts_name = GetFileSpec(fts_residue,part='name',/strict)
mode = 'zld'

ff = 'c3avg_2003-8rolledoff.grd'
source_map = filepath(root=source,ff)
message, /info, hide_env(source_map)
IF NOT grd_read(source_map,zldhalf,/silent) THEN	$
	RETURN

zldhalf = reverse(zldhalf)			; 360 x 360

zldmap = fltarr(720,360)			; 720 x 360
zldmap[179+1:179+360,*] = zldhalf	; Insert in middle

n = size(zldmap,/dim)

; Now, create FITS headers for the maps.

message, /info, 'create main FITS header'

mkhdr, hdr, IsType(zldmap), n, /extend

fxaddpar, hdr, 'NAME'	,	fts_name, ' File name'
fxaddpar, hdr, 'IMG_TYPE',   'skymap_'+mode, ' Image type'
fxaddpar, hdr, 'MAP' 	, 'Helioecliptic zld residue HLON=[-90,+90],HLAT=[-90,+90]', ' Map description'
fxaddpar, hdr, 'RADESYS', 	   'FK5', ' IAU 1984 system'
fxaddpar, hdr, 'EQUINOX',     2000.0, ' J2000 coordinates'				, format='f8.1'
fxaddpar, hdr, 'CTYPE1'	, 	  'HLON', ' Axis type'
fxaddpar, hdr, 'CTYPE2'	,	  'HLAT', ' Axis type'
fxaddpar, hdr, 'CUNIT1'	,	   'deg', ' Angles in degrees'
fxaddpar, hdr, 'CUNIT2'	,	   'deg', ' Angles in degrees'
fxaddpar, hdr, 'CRPIX1'	,n[0]/2+0.50, ' Array index (base 1) for HLON=0', format='f8.2'
fxaddpar, hdr, 'CRPIX2'	,n[1]/2+0.50, ' Array index (base 1) for HLAT=0', format='f8.2'
fxaddpar, hdr, 'CDELT1'	,	    0.50, ' Degrees per bin'				, format='f8.2'
fxaddpar, hdr, 'CDELT2'	,	    0.50, ' Degrees per bin'				, format='f8.2'
fxaddpar, hdr, 'CRVAL1'	,		0.00, ' HLON=0.0'						, format='f8.2'
fxaddpar, hdr, 'CRVAL2' ,		0.00, ' HLAT=0.0'						, format='f8.2'
fxaddpar, hdr, 'BAD_DATA',		0.00, ' Missing data flag'				, format='f8.2'

;view, in=zldmap,/stretch

message, /info, 'write '+hide_env(fts_residue)
writefits, fts_residue, zldmap, hdr

hdr0 = hdr

message, /info, 'create FITS header for extensions'

mkhdr, hdr, IsType(zldmap), size(zldmap,/dim), /image

iend = (where(strpos(hdr ,'END     ') EQ 0))[0]
imap = (where(strpos(hdr0,'MAP     ') EQ 0))[0]
ibad = (where(strpos(hdr0,'BAD_DATA') EQ 0))[0]

hdr  = [hdr[0:iend-1],hdr0[imap:ibad],hdr[iend]]

FOR iweek=1,52 DO BEGIN

	cweek = string(iweek,format='(I2.2)')
	ff = 'cmap'+cweek+'_2003through2008.grd'

	; Abort if the grd file cannot be read

	source_map = filepath(root=source,ff)
	message, /info, hide_env(source_map)
	IF NOT grd_read(source_map,zldhalf,/silent) THEN BEGIN
		spawn, 'rm -vf '+fts_residue
		RETURN
	ENDIF

	zldhalf = reverse(zldhalf)

	zldmap = fltarr(720,360)
	zldmap[179+1:179+360,*] = zldhalf

	fxaddpar , hdr, 'MAP', 'Helioecliptic zld residue (week '+cweek+') HLON=[-90,+90],HLAT=[-90,+90]', ' Map description'
	writefits, fts_residue, zldmap, hdr, /append

ENDFOR

spawn, 'gzip -vf '+fts_residue

RETURN  &  END