;+
; NAME:
;	smei_zld_remove
; PURPOSE:
;	Removes zodiacal light distribution from SMEI orbital sky maps
; CATEGORY:
;	ucsd/camera/idl/zld
; CALLING SEQUENCE:
	PRO smei_zld_remove, wanted_map , $

		destination = destination	, $
		to_smeidb	= to_smeidb		, $

		checkversion= checkversion	, $
		overwrite	= overwrite 	, $
		force		= force			, $
		silent		= silent		, $

		width		= width 		, $
		reverse_loop= reverse_loop	, $

		_extra		= _extra		, $
		keep_gegenschein = keep_gegenschein

; INPUTS:
;	wanted_map		array[n]; type: see href=smei_getfile=
;						Selects skymap files to be processed
;						in combination with other keywords passed
;						using _extra.
; OPTIONAL INPUT PARAMETERS:
;
;	/to_smeidb		if set, write to directory structure as used
;						for SMEI data base (separate subdirs for mode
;						and camera) with 'destination' as the base dir.
;	destination=destination
;					scalar; type: string
;						destination directory for skymaps
;						if /to_smeidb NOT set then maps are written
;							directly to 'destination'.
;							The default in this case is $TUB
;						if /to_smeidb SET then files are written to
;							a directory structure as used for the
;							SMEI data base with 'destination' as the
;							base directory.
;							The default in this case is 'SMEIDB?' (the
;							SMEI database).
;						If destination starts with a "user@host:" identifier then
;						the output is first created locally in $TUB
;						and is then transferred to the remote host using scp;
;						the local files will be deleted after the transfer;
;						note that this will only work if the remote host has
;						been configured for ssh using private/public keys
;						allowing login without passwords.
;	/silent 		if set, suppresses informational messages
;
;	/checkversion	overwrite existing skymaps only if the version
;					number has increased.
;	/overwrite		overwrite existing skymaps
;
;	_extra=_extra	keywords passed to href=smei_getfile=
; OUTPUTS:
;	Files in destination directory
; INCLUDE:
	@compile_opt.pro		; On error, return to caller
; CALLS:
;	InitVar, CheckDir, IsType, smei_filename, ArrayLocation
;	GetFileSpec, gzip_file, smei_camera_gain
;	hide_env, smei_zld_model, gridgen, BadValue, smei_sky, big_eph
;	jpl_body, CvSky, TimeSet, TimeGet, TimeOp, TimeUnit, TimeSystem
; PROCEDURE:
; MODIFICATION HISTORY:
;	SEP-2005, Paul Hick (UCSD/CASS), V1.0
;		Started, but never finished
;	DEC-2007, Paul Hick (UCSD/CASS), v2.0
;		First working version. Now uses /model_2 in smei_zld_model
;	FEB-2008, Paul Hick (UCSD/CASS)
;		Location of files now done using smei_getfile.
;	APR-2008, Paul Hick (UCSD/CASS; ppick@ucsd.edu)
;		Rewrite. Not using smoothsphere to reduce resolution
;		(too slow). Instead the raw maps (equatorial, polar)
;		are now smoothed using the IDL smooth function. Then
;		The sun-centered ecliptic map is obtained using
;		bi-linear interpolation on the smoothed raw maps.
;		Version number will be increased once the next
;		(final?) zld model becomes available.
;	SEP-2008, Paul  Hick (UCDS/CASS; pphick@ucsd.edu)
;		Added code to check star subtraction version and
;		sidereal background version when /checkversion is
;		set.
;-
; Version number for the default zodiacal model

zld_version = smei_zld_model(/get_version)
message, /info, 'zodiacal light distribution, model'+string(zld_version,'(F5.2)')

InitVar, to_smeidb  , /key
InitVar, destination, ([getenv('TUB'),'SMEIDB?'])[to_smeidb]

InitVar, silent 	 , 0
InitVar, checkversion, /key
InitVar, overwrite	 , /key
InitVar, force		 , /key

overwrite OR= force

InitVar, width, 1.0
InitVar, reverse_loop, /key

from_mode = 'equ'
to_mode	  = 'ecl'

; Set up local and remote destination
; The presence of a colon indicates that transfer to a remote machine is requested

to_remote = strpos(destination,':') NE -1

CASE to_remote OF
0: local_destination  = destination
1: BEGIN				; Write to remote directory
	local_destination  = getenv('TUB')	
	remote_destination = destination
END
ENDCASE

; Selection of files. Returns fully-qualified filenames.

given_map = smei_getfile(wanted_map, $
	mode	= from_mode	, $
	count	= count		, $
	_extra	= _extra)

IF count EQ 0 THEN	$
	RETURN


stop_time = TimeSystem(/silent)

IF silent LE 3 THEN message, /info, strcompress(count,/rem)+' '+from_mode+' map(s)'

cvsmei_init, /get, dim=ndim_lo, mode='lores'
ndim_eq = ndim_lo

; Lores sun-centered ecliptic grids

range_lo = [[-180.0,180.0],[-90.0,90.0]]
grid_lo = gridgen(ndim_lo,range=range_lo,/open)

crpix1_lo = -0.5d0-range_lo[0]/(range_lo[1]-range_lo[0])*ndim_lo[0]
crpix2_lo = -0.5d0-range_lo[2]/(range_lo[3]-range_lo[2])*ndim_lo[1]
cdelt1_lo = double(range_lo[1]-range_lo[0])/ndim_lo[0]
cdelt2_lo = cdelt1_lo

crpix1_lo += 1			; Change to 1-based index
crpix2_lo += 1

; Lores equatorial grid

range_eq = [[0.0,360.0],[-90.0,90.0]]
grid_eq = gridgen(ndim_eq,range=range_eq,/open)

crpix1_eq = -0.5d0;-range_eq[0]/(range_eq[1]-range_eq[0])*ndim_eq[0]
crpix2_eq = -0.5d0-range_eq[2]/(range_eq[3]-range_eq[2])*ndim_eq[1]
cdelt1_eq = double(range_eq[1]-range_eq[0])/ndim_eq[0]
cdelt2_eq = cdelt1_eq

crpix1_eq += 1			; Change to 1-based index
crpix2_eq += 1


badval  = BadValue(0.0)
missing = badval		; "Missing" value for smooth fnc

dec_cutoff = 55.0
box_cutoff = 60.0-dec_cutoff

rpd = !dpi/180.0d0

imap_low  =		  0L *(1-reverse_loop)+(count-1)*reverse_loop
imap_high = (count-1)*(1-reverse_loop)+		 0L *reverse_loop
imap_step =		  1L *(1-reverse_loop)+(	-1L)*reverse_loop

FOR imap=imap_low,imap_high,imap_step DO BEGIN

	start_time = stop_time

	skyfile = given_map[imap]
	skyhide = hide_env(skyfile)

	skytime = smei_filename(skyfile,camera=camera,postfix=postfix)

	workname = smei_filename(skytime,camera=camera,mode=to_mode,type='.fts')

	CASE 1 OF
	to_remote: tmp = local_destination			; Local destination is $TUB
	to_smeidb: tmp = smei_filepath(source=local_destination,mode=to_mode,camera=camera,postfix=postfix)
	ELSE	 : tmp = local_destination
	ENDCASE

	IF NOT CheckDir(tmp) THEN continue

	workfile = filepath(root=tmp,workname)

	IF checkversion THEN BEGIN

		; Extract version numbers from 'equ' file

		hdr = headfits(skyfile)

		equ_map_bstar_version = fxpar(hdr,'BSTARVER')
		equ_has_bstar_version = !err EQ 0

		equ_map_bkgnd_version = fxpar(hdr,'BKGNDVER')
		equ_has_bkgnd_version = !err EQ 0

		; Extract version numbers from 'ecl' file

		tmp = (file_search(workfile+'*'))[0]

		IF tmp NE '' THEN BEGIN 	; Skymap exists

			; Extract version numbers from 'ecl' file

			hdr = headfits(tmp)

			map_zld_version = fxpar(hdr,'ZLDVER')
			has_zld_version = !err EQ 0

			map_bstar_version = fxpar(hdr,'BSTARVER')
			has_bstar_version = !err EQ 0

			map_bkgnd_version = fxpar(hdr,'BKGNDVER')
			has_bkgnd_version = !err EQ 0

			process_map = 0

			; If the version of the bright star subtraction
			; in the 'ecl' file is lower than in the 'equ' file
			; the zld subtraction needs to be redone

			IF NOT process_map AND has_bstar_version AND equ_has_bstar_version THEN BEGIN
				process_map = map_bstar_version LT equ_map_bstar_version
				IF process_map THEN 		$
					IF silent LE 2 THEN		$
						message, /info, 'old "star subtraction", '+hide_env(workfile)
			ENDIF

			; If the version of the sidereal background subtraction
			; in the 'ecl' file is lower than in the 'equ' file
			; the zld subtraction needs to be redone

			IF NOT process_map AND has_bkgnd_version AND equ_has_bkgnd_version THEN BEGIN
				process_map = map_bkgnd_version LT equ_map_bkgnd_version
				IF process_map THEN 		$
					IF silent LE 2 THEN		$
						message, /info, 'old "sidereal background", '+hide_env(workfile)
			ENDIF

			; If the zld subtraction version in the 'ecl' is
			; lower than the software version then the zld
			; subtraction needs to be redone

			IF NOT process_map AND has_zld_version THEN BEGIN
				process_map = map_zld_version LT zld_version
				IF silent LE 2 THEN 		$
					IF NOT process_map THEN $
						message, /info, '>= version, '+hide_env(workfile)
			ENDIF

			; CONTINUE not allowed inside CASE block

			IF NOT process_map THEN CONTINUE

		ENDIF

	ENDIF ELSE IF NOT overwrite THEN BEGIN

		tmp = (file_search(workfile+'*'))[0]
		IF tmp NE '' THEN BEGIN 	; Skymap exists
			IF silent LE 2 THEN message, /info, hide_env(tmp)+' exists'
			CONTINUE
		ENDIF

	ENDIF

	IF silent LE 2 THEN message, /info, skyhide


	; ======================
	; Process equatorial map

	smei_sky, skyfile, /exists, /equatraw , sky=skybump, nbin=1, /noplot, exten_no=ext_eq, hdr=hdr, silent=silent+1
	smei_sky, skyfile, /exists, /equatflat, sky=skyflat, nbin=1, /noplot, exten_no=ext_eqflat, silent=silent+1

	tt = TimeSet(smei=hdr.orbit)
	gain = smei_camera_gain(tt, camera=camera)

	; Subtract zodiacal light model
	; (limiting the zodiac calculation to finite bins in skybump is not
	; necessary, but should speed up things since big sections of the
	; single-camera skymaps will not be finite)

	hdr = headfits(skyfile, exten=ext_eq)
	tmp = fxpar(hdr,'ZLD')
	IF !err EQ -1 THEN tmp = 0
	remove_zld = 1-tmp

	IF remove_zld THEN BEGIN
		i = where(finite(skybump))
		IF i[0] NE -1 THEN BEGIN
			tmp = ArrayLocation(i,dim=size(skybump,/dim))
			tmp = cvsmei(from_map=tmp,mode='eqmap',/to_equatorial,/degrees,/silent)
			tmp = smei_zld_model(skytime,tmp,/equatorial,/degrees,/silent,keep_gegenschein=keep_gegenschein)*gain
			skybump[i] -= tmp
			skyflat[i] -= tmp
		ENDIF
	ENDIF

	; Reduce resolution of hires equatorial map by factor 'rebin'

	; We will do a boxcar smoothing with the IDL smooth function.
	; The north and south side of the map overlap with the polar
	; map by at least 10 degrees. This is enough to avoid problems
	; with the edge effects of the smooth function as long as the
	; smoothing box is less than 5 degrees wide.

	; On the left (RA=0) and right (RA=360) side of the map we pad
	; the skymap by 5 degrees by copying a strip from the opposite
	; side of the maps. This is stripped of again after smooth
	; has been called.

	cvsmei_init,/get,mpb=mpb,mode='eqmap',/degrees
	npad = round(box_cutoff/mpb)		; Pad width (should be 50)

	ndim = size(skybump,/dim)
	skybump = [ skybump[ndim[0]-npad:ndim[0]-1,*], skybump, skybump[0:npad-1,*] ]
	skyflat = [ skyflat[ndim[0]-npad:ndim[0]-1,*], skyflat, skyflat[0:npad-1,*] ]

	; Use a width-degree smoothing window (+/- 0.5*width degrees)
	; Note that the smoothing window needs to stay inside the
	; padding zone of npad bins. This puts an upper limit of
	; 2*box_cutoff degrees on the smoothing box.

	nsmooth = (width/mpb+1) < (2*npad+1)

	skybump = smooth(skybump,nsmooth,/nan)
	skyflat = smooth(skyflat,nsmooth,/nan)

	; Remove the padding of npad bins left and right again

	skybump = skybump[npad:ndim[0]+npad-1,*]
	skyflat = skyflat[npad:ndim[0]+npad-1,*]

	; Location Sun in geocentric ecliptic coordinates

	loc_sun = big_eph( tt, $
		body	= jpl_body(/sun  ,/string)	, $
		center	= jpl_body(/earth,/string)	, $
		/to_ecliptic, $
		/to_sphere	, $
		/degrees	, $
		/precess	, $
		/onebody	, $
		/silent 	)


	; The output map will be at the lores resolution in
	; sun-centered ecliptic coordinates.

	; Convert lores grid from sun-centered ecliptic to regular
	; ecliptic coordinates by adding ecliptic longitude of Sun

	grid = grid_lo
	grid[0,*] += loc_sun[0]

	; Convert from ecliptic to equatorial coordinates.
	; These are the coordinates used in the skymaps (equatorial
	; and polar maps with different map projections).

	grid = CvSky(tt,from_ecliptic=grid,/to_equatorial,/degrees,/silent)

	; Set up lores sun-centered full-sky map

	sun_skybump_ecl = replicate(badval, ndim_lo)
	sun_skyflat_ecl = sun_skybump_ecl

	; Elements in grid with abs(dec) <= dec_cutoff are
	; interpolated from the smoothed equatorial map

	i = where( abs(grid[1,*]) LE dec_cutoff, n )
	IF i[0] NE -1 THEN BEGIN
		tmp = cvsmei(from_equatorial=grid[*,i],/degrees,/to_map,mode='eqmap',/silent)
		sun_skybump_ecl[i] = interpolate(skybump, tmp[0,*], tmp[1,*], missing=missing)
		sun_skyflat_ecl[i] = interpolate(skyflat, tmp[0,*], tmp[1,*], missing=missing)
	ENDIF

	; set up lores equatorial full-sky map

	sun_skybump_equ = replicate(badval, ndim_eq)
	sun_skyflat_equ = sun_skybump_equ

	i = where( abs(grid_eq[1,*]) LE dec_cutoff, n )
	IF i[0] NE -1 THEN BEGIN
		tmp = cvsmei(from_equatorial=grid_eq[*,i],/degrees,/to_map,mode='eqmap',/silent)
		sun_skybump_equ[i] = interpolate(skybump, tmp[0,*], tmp[1,*], missing=missing)
		sun_skyflat_equ[i] = interpolate(skyflat, tmp[0,*], tmp[1,*], missing=missing)
	ENDIF

	; =========================
	; Process map of north pole

	smei_sky, skyfile, /exists, /northraw , sky=skybump, nbin=1, /noplot, exten_no=ext_np, hdr=hdr, silent=silent+1
	smei_sky, skyfile, /exists, /northflat, sky=skyflat, nbin=1, /noplot, exten_no=ext_npflat, silent=silent+1

	; Subtract zodiacal light model

	IF remove_zld THEN BEGIN
		i = where(finite(skybump))
		IF i[0] NE -1 THEN BEGIN
			tmp = ArrayLocation(i,dim=size(skybump,/dim))
			tmp = cvsmei(from_map=tmp,mode='npmap',/to_equatorial,/degrees,/silent)
			tmp = smei_zld_model(skytime,tmp,/equatorial,/degrees,/silent,keep_gegenschein=keep_gegenschein)*gain
			skybump[i] -= tmp
			skyflat[i] -= tmp
		ENDIF
	ENDIF

	; Smooth the north polar map.

	skybump = smooth(skybump,nsmooth,/nan)
	skyflat = smooth(skyflat,nsmooth,/nan)


	; Elements in grid with dec > dec_cutoff are
	; interpolated from the smoothed north polar map

	i = where( grid[1,*] GT dec_cutoff, n )
	IF i[0] NE -1 THEN BEGIN
		tmp = cvsmei(from_equatorial=grid[*,i],/degrees,/to_map,mode='npmap',/silent)
		sun_skybump_ecl[i] = interpolate(skybump, tmp[0,*], tmp[1,*], missing=missing)
		sun_skyflat_ecl[i] = interpolate(skyflat, tmp[0,*], tmp[1,*], missing=missing)
	ENDIF

	i = where( grid_eq[1,*] GT dec_cutoff, n )
	IF i[0] NE -1 THEN BEGIN
		tmp = cvsmei(from_equatorial=grid_eq[*,i],/degrees,/to_map,mode='npmap',/silent)
		sun_skybump_equ[i] = interpolate(skybump, tmp[0,*], tmp[1,*], missing=missing)
		sun_skyflat_equ[i] = interpolate(skyflat, tmp[0,*], tmp[1,*], missing=missing)
	ENDIF

	; =========================
	; Process map of south pole

	smei_sky, skyfile, /exists, /southraw , sky=skybump, nbin=1, /noplot, exten_no=ext_sp, hdr=hdr, silent=silent+1
	smei_sky, skyfile, /exists, /southflat, sky=skyflat, nbin=1, /noplot, exten_no=ext_spflat, silent=silent+1

	; Subtract zodiacal light model

	IF remove_zld THEN BEGIN
		i = where(finite(skybump))
		IF i[0] NE -1 THEN BEGIN
			tmp = ArrayLocation(i,dim=size(skybump,/dim))
			tmp = cvsmei(from_map=tmp,mode='spmap',/to_equatorial,/degrees,/silent)
			tmp = smei_zld_model(skytime,tmp,/equatorial,/degrees,/silent,keep_gegenschein=keep_gegenschein)*gain
			skybump[i] -= tmp
			skyflat[i] -= tmp
		ENDIF
	ENDIF

	; Smooth the south polar map.

	skybump = smooth(skybump,nsmooth,/nan)
	skyflat = smooth(skyflat,nsmooth,/nan)

	; Elements in grid with dec < -dec_cutoff are
	; interpolated from the smoothed south polar map

	i = where( grid[1,*] LT -dec_cutoff, n )
	IF i[0] NE -1 THEN BEGIN
		tmp = cvsmei(from_equatorial=grid[*,i],/degrees,/to_map,mode='spmap',/silent)
		sun_skybump_ecl[i] = interpolate(skybump, tmp[0,*], tmp[1,*], missing=missing)
		sun_skyflat_ecl[i] = interpolate(skyflat, tmp[0,*], tmp[1,*], missing=missing)
	ENDIF

	i = where( grid_eq[1,*] LT -dec_cutoff, n )
	IF i[0] NE -1 THEN BEGIN
		tmp = cvsmei(from_equatorial=grid_eq[*,i],/degrees,/to_map,mode='spmap',/silent)
		sun_skybump_equ[i] = interpolate(skybump, tmp[0,*], tmp[1,*], missing=missing)
		sun_skyflat_equ[i] = interpolate(skyflat, tmp[0,*], tmp[1,*], missing=missing)
	ENDIF


	; ====================
	; Process low res maps

	smei_sky, skyfile, /exists, /orbtime, hdr=hdr_time, sky=sky_time, /noplot, exten_no=ext_time, silent=silent+1
	smei_sky, skyfile, /exists, /orbsecs, hdr=hdr_secs, sky=sky_secs, /noplot, exten_no=ext_secs, silent=silent+1
	smei_sky, skyfile, /exists, /hit    , hdr=hdr_hits, sky=sky_hits, /noplot, exten_no=ext_hits, silent=silent+1

	; Convert from RA,dec to array indices in lores map.
	; 'grid' at this point is the lores sun-centered ecliptic grid
	; converted to equatorial (RA/dec).

	grid = cvsmei(from_equatorial=grid, /degrees, /to_map, mode='lores',/silent)

	; Note that this linear interpolation does not work very well on the "edge" of the
	; skymap (where start and end of the orbit meet). In equatorial coordinates the
	; edge is sharp with the time jumping by 1 orbital period from one pixel to the
	; next. The interpolation will cause invalid times to be inserted in the map
	; close to this transition.

	sky_time = interpolate(sky_time, grid[0,*], grid[1,*], missing=missing)
	sky_secs = interpolate(sky_secs, grid[0,*], grid[1,*], missing=missing)
	sky_hits = interpolate(sky_hits, grid[0,*], grid[1,*], missing=missing)

	sky_time = reform(sky_time,ndim_lo,/overwrite)
	sky_secs = reform(sky_secs,ndim_lo,/overwrite)
	sky_hits = reform(sky_hits,ndim_lo,/overwrite)


	; ====================
	; Write the Fits files

	; Mark bad data with BAD_DATA value from Fits header

	hdr = headfits(skyfile, exten=ext_eq)
	bad_data = fxpar(hdr,'BAD_DATA')

	bad = where(1-finite(sun_skybump_ecl))
	IF bad[0] NE -1 THEN sun_skybump_ecl[bad] = bad_data

	bad = where(1-finite(sun_skyflat_ecl))
	IF bad[0] NE -1 THEN sun_skyflat_ecl[bad] = bad_data

	bad = where(1-finite(sky_time))
	IF bad[0] NE -1 THEN sky_time[bad] = bad_data

	bad = where(1-finite(sky_secs))
	IF bad[0] NE -1 THEN sky_secs[bad] = bad_data

	bad = where(1-finite(sky_hits))
	IF bad[0] NE -1 THEN sky_hits[bad] = bad_data

	; Set up Fits header template

	; Make main header. Copy most entries from the main header of the original skymap
	; Fits keyword SMEI_HTM comes right after the mandatory keywords.

	; Main header updates

	fxaddpar, hdr, 'DATE', TimeGet(TimeSystem(/silent),format='YYYY-MN-DD hh:mm:ss',/scalar)
	fxaddpar, hdr, 'MAP' , 'Helioecliptic map HLNG=[-180,180],HLAT=[-90,+90]', ' Map description'

	fxaddpar, hdr, 'NAXIS1', ndim_lo[0]
	fxaddpar, hdr, 'NAXIS2', ndim_lo[1]

	tmp = fxpar(hdr, 'CRPIX1')
	IF !err GE 0 THEN BEGIN
		fxaddpar, hdr, 'CRPIX1', crpix1_lo	, format='(F10.2)', ' Array index (base 1) for HLON=0'
		fxaddpar, hdr, 'CRPIX2', crpix2_lo	, format='(F10.2)', ' Array index (base 1) for HLAT=0'
		fxaddpar, hdr, 'CDELT1', cdelt1_lo	, format='(F10.2)'
		fxaddpar, hdr, 'CDELT2', cdelt2_lo	, format='(F10.2)'
		fxaddpar, hdr, 'CTYPE1', 'HLON', ' Axis type: helioecliptic longitude'
		fxaddpar, hdr, 'CTYPE2', 'HLAT', ' Axis type: helioecliptic latitude'
		fxaddpar, hdr, 'CRVAL1',   0.0 , ' HLON=0.0', format='(F10.2)' 
		fxaddpar, hdr, 'CRVAL2',   0.0 , ' HLAT=0.0', format='(F10.2)'
	ENDIF ELSE BEGIN
		fxaddpar, hdr, 'CXEQ', crpix1_lo	, format='(F10.2)', ' Array index (base 1) for HLON=0'
		fxaddpar, hdr, 'CYEQ', crpix2_lo	, format='(F10.2)', ' Array index (base 1) for DEC=0'
		fxaddpar, hdr, 'D_EQ', 1.0d0/cdelt1_lo, format='(F10.2)'
	ENDELSE

	fxaddpar, hdr, 'NAME'	 , GetFileSpec(workfile,part='name')
	fxaddpar, hdr, 'IMG_TYPE', 'skymap_'+to_mode

	; Main header additions. Add new keywords in reverse order as needed in Fits file.

	IF remove_zld THEN BEGIN
		fxaddpar, hdr, 'ZLD'	, fix(remove_zld), ' Zodiacal light subtracted (0=N,1=Y)?'	, after='BAD_DATA'
		fxaddpar, hdr, 'ZLDVER'	, zld_version 	 , ' Zodiacal light model'	, after='ZLD', format='(F5.2)'
	ENDIF

	fxaddpar, hdr, 'SUN_DIS', loc_sun[2], ' [AU] Sun-Earth distance'		, after='BAD_DATA', format='(F10.6)'
	fxaddpar, hdr, 'SUN_LAT', loc_sun[1], ' [deg] Ecliptic latitude Sun'	, after='BAD_DATA', format='(F10.6)'
	fxaddpar, hdr, 'SUN_LNG', loc_sun[0], ' [deg] Ecliptic longitude Sun'	, after='BAD_DATA', format='(F10.5)'

	; ===========
	; Start writing file

	; Extension 0, Map without zodiacal light HLNG=[-180,180],HLAT=[-90,+90]

	writefits, workfile, sun_skybump_ecl, hdr

	; Extension 1, Map (flat) without zodiacal light HLNG=[-180,180],HLAT=[-90,+90]

	; Extract the part of the header containing the keywords
	; pertaining to the map projection, adding at the start the
	; 'MAP' keyword, and at the end the 'BAD_DATA' keyword
	; This includes keywords
	; MAP,RADESYS,EQUINOX,CTYPE1,CTYPE2,CUNIT1,CUNIT2,CRPIX1,CRPIX2,
	; CDELT1,CDELT2,CRVAL1,CRVAL2,BAD_DATA

	hdr_proj = hdr[(where(strpos(hdr,'MAP') EQ 0))[0]:(where(strpos(hdr,'BAD_DATA') EQ 0))[0]]

	; Set up a basic header; then insert hdr_proj before 'END'

	mkhdr, hdr, IsType(sun_skyflat_ecl), ndim_lo, /image
	hdr = [hdr[0:(where(strpos(hdr,'END   ') EQ 0))[0]-1],hdr_proj]
	fxaddpar, hdr, 'MAP', 'Helioecliptic map HLNG=[-180,180],HLAT=[-90,+90] (flat)', ' Map description'

	writefits, workfile, sun_skyflat_ecl, hdr, /append

	; Copy couple of lowres extensions from original skymap

	; Extension 2, Time map (fraction of orbit) HLNG=[-180,180],HLAT=[-90,+90

	hdr = headfits(skyfile, exten=ext_time)

	fxaddpar, hdr, 'NAXIS1', ndim_lo[0]
	fxaddpar, hdr, 'NAXIS2', ndim_lo[1]

	tmp = fxpar(hdr, 'CRPIX1')
	IF !err GE 0 THEN BEGIN
		fxaddpar, hdr, 'CRPIX1', crpix1_lo	, format='(F10.2)', ' Array index (base 1) for HLON=0'
		fxaddpar, hdr, 'CRPIX2', crpix2_lo	, format='(F10.2)', ' Array index (base 1) for HLAT=0'
		fxaddpar, hdr, 'CDELT1', cdelt1_lo	, format='(F10.2)'
		fxaddpar, hdr, 'CDELT2', cdelt2_lo	, format='(F10.2)'
		fxaddpar, hdr, 'CTYPE1', 'HLON',' Axis type: helioecliptic longitude'
		fxaddpar, hdr, 'CTYPE2', 'HLAT',' Axis type: helioecliptic latitude'
		fxaddpar, hdr, 'CRVAL1',   0.0 ,' HLON=0.0', format='(F10.2)'
		fxaddpar, hdr, 'CRVAL2',   0.0 ,' HLAT=0.0', format='(F10.2)'
	ENDIF ELSE BEGIN
		fxaddpar, hdr, 'CXLO', crpix1_lo	, format='(F10.2)', ' Array index (base 1) for HLON=0'
		fxaddpar, hdr, 'CYLO', crpix2_lo	, format='(F10.2)', ' Array index (base 1) for HLAT=0'
		fxaddpar, hdr, 'D_LO', 1.0d0/cdelt1_lo, format='(F10.2)'
	ENDELSE

	fxaddpar, hdr, 'MAP', 'Time (fraction of orbit) HLNG=[-180,180],HLAT=[-90,+90]'

	writefits, workfile, sky_time, hdr, /append

	; Extension 3, Time map (sec since 1st frame used) HLNG=[-180,180],HLAT=[-90,+90

	hdr = headfits(skyfile, exten=ext_secs)

	fxaddpar, hdr, 'NAXIS1', ndim_lo[0]
	fxaddpar, hdr, 'NAXIS2', ndim_lo[1]

	tmp = fxpar(hdr, 'CRPIX1')
	IF !err GE 0 THEN BEGIN
		fxaddpar, hdr, 'CRPIX1', crpix1_lo	, format='(F10.2)', ' Array index (base 1) for HLON=0'
		fxaddpar, hdr, 'CRPIX2', crpix2_lo	, format='(F10.2)', ' Array index (base 1) for HLAT=0'
		fxaddpar, hdr, 'CDELT1', cdelt1_lo	, format='(F10.2)'
		fxaddpar, hdr, 'CDELT2', cdelt2_lo	, format='(F10.2)'
		fxaddpar, hdr, 'CTYPE1', 'HLON',' Axis type: helioecliptic longitude'
		fxaddpar, hdr, 'CTYPE2', 'HLAT',' Axis type: helioecliptic latitude'
		fxaddpar, hdr, 'CRVAL1',   0.0 ,' HLON=0.0', format='(F10.2)'
		fxaddpar, hdr, 'CRVAL2',   0.0 ,' HLAT=0.0', format='(F10.2)'
	ENDIF ELSE BEGIN
		fxaddpar, hdr, 'CXLO', crpix1_lo	, format='(F10.2)', ' Array index (base 1) for HLON=0'
		fxaddpar, hdr, 'CYLO', crpix2_lo	, format='(F10.2)', ' Array index (base 1) for HLAT=0'
		fxaddpar, hdr, 'D_LO', 1.0d0/cdelt1_lo, format='(F10.2)'
	ENDELSE

	fxaddpar, hdr, 'MAP', 'Time (sec since TORIGIN) HLNG=[-180,180],HLAT=[-90,+90]'

	writefits, workfile, sky_secs, hdr, /append

	; Extension 4, Contributing-pixel count map HLNG=[-180,180],HLAT=[-90,+90

	hdr = headfits(skyfile, exten=ext_hits)

	fxaddpar, hdr, 'NAXIS1', ndim_lo[0]
	fxaddpar, hdr, 'NAXIS2', ndim_lo[1]

	tmp = fxpar(hdr, 'CRPIX1')
	IF !err GE 0 THEN BEGIN
		fxaddpar, hdr, 'CRPIX1', crpix1_lo	, format='(F10.2)', ' Array index (base 1) for HLON=0'
		fxaddpar, hdr, 'CRPIX2', crpix2_lo	, format='(F10.2)', ' Array index (base 1) for HLAT=0'
		fxaddpar, hdr, 'CDELT1', cdelt1_lo	, format='(F10.2)'
		fxaddpar, hdr, 'CDELT2', cdelt2_lo	, format='(F10.2)'
		fxaddpar, hdr, 'CTYPE1', 'HLON',' Axis type: helioecliptic longitude'
		fxaddpar, hdr, 'CTYPE2', 'HLAT',' Axis type: helioecliptic latitude'
		fxaddpar, hdr, 'CRVAL1',   0.0 ,' HLON=0.0', format='(F10.2)'
		fxaddpar, hdr, 'CRVAL2',   0.0 ,' HLAT=0.0', format='(F10.2)'
	ENDIF ELSE BEGIN
		fxaddpar, hdr, 'CXlO', crpix1_lo	, format='(F10.2)', ' Array index (base 1) for HLON=0'
		fxaddpar, hdr, 'CYLO', crpix2_lo	, format='(F10.2)', ' Array index (base 1) for HLAT=0'
		fxaddpar, hdr, 'D_LO', 1.0d0/cdelt1_lo, format='(F10.2)'
	ENDELSE

	fxaddpar, hdr, 'MAP', 'Contributing pixel count HLNG=[-180,180],HLAT=[-90,+90]'

	writefits, workfile, sky_hits, hdr, /append

	; Add lores equatorial maps

	; Extension 5, Map without zodiacal light RA=[0,360],DEC=[-90,+90]

	mkhdr, hdr, IsType(sun_skybump_equ), ndim_eq, /image
	hdr = [hdr[0:(where(strpos(hdr,'END   ') EQ 0))[0]-1],hdr_proj]
	fxaddpar, hdr, 'MAP', 'Equatorial map RA=[0,360],DEC=[-90,+90]', ' Map description'

	fxaddpar, hdr, 'NAXIS1', ndim_eq[0]
	fxaddpar, hdr, 'NAXIS2', ndim_eq[1]

	tmp = fxpar(hdr, 'CRPIX1')
	IF !err GE 0 THEN BEGIN
		fxaddpar, hdr, 'CRPIX1', crpix1_eq	, format='(F10.2)', ' Array index (base 1) for HLON=0'
		fxaddpar, hdr, 'CRPIX2', crpix2_eq	, format='(F10.2)', ' Array index (base 1) for HLAT=0'
		fxaddpar, hdr, 'CDELT1', cdelt1_eq	, format='(F10.2)'
		fxaddpar, hdr, 'CDELT2', cdelt2_eq	, format='(F10.2)'
		fxaddpar, hdr, 'CTYPE1', 'RA' , ' Axis type'
		fxaddpar, hdr, 'CTYPE2', 'DEC', ' Axis type'
		fxaddpar, hdr, 'CRVAL1',   0.0 , ' RA=0.0' , format='(F10.2)' 
		fxaddpar, hdr, 'CRVAL2',   0.0 , ' DEC=0.0', format='(F10.2)'
	ENDIF ELSE BEGIN
		fxaddpar, hdr, 'CXEQ', crpix1_eq	, format='(F10.2)', ' Array index (base 1) for HLON=0'
		fxaddpar, hdr, 'CYEQ', crpix2_eq	, format='(F10.2)', ' Array index (base 1) for DEC=0'
		fxaddpar, hdr, 'D_EQ', 1.0d0/cdelt1_eq, format='(F10.2)'
	ENDELSE

	writefits, workfile, sun_skybump_equ, hdr, /append

	; Extension 6, Map (flat) without zodiacal light RA=[0,360],DEC=[-90,+90]

	fxaddpar, hdr, 'MAP', 'Equatorial map (flat) RA=[0,360],DEC=[-90,+90]', ' Map description'
	writefits, workfile, sun_skyflat_equ, hdr, /append

	; Extension 7, time (fraction of orbit) RA=[0,360],DEC=[-90,+90]

	tmp = readfits(skyfile, hdr, exten_no=ext_time, /silent)
	writefits, workfile, tmp, hdr, /append

	; Extension 8, time (sec since TORIGIN) RA=[0,360],DEC=[-90,+90]

	tmp = readfits(skyfile, hdr, exten_no=ext_secs, /silent)
	writefits, workfile, tmp, hdr, /append

	; Extension 9, contributing pixel count RA=[0,360],DEC=[-90,+90]

	tmp = readfits(skyfile, hdr, exten_no=ext_hits, /silent)
	writefits, workfile, tmp, hdr, /append

	; File complete
	; =============
	; /compress on writefits is ignored with /append present, so gzip here.

	tmp = gzip_file(workfile, /force)
	workfile += '.gz'
	message, /info, hide_env(workfile)+' is '+to_mode+' map'


	IF to_remote THEN BEGIN

		tmp = strpos(remote_destination,':')+1

		CASE to_smeidb OF
		0: tmp = remote_destination
		1: tmp = strmid(remote_destination,0,tmp)	+ $
			smei_filepath(source=strmid(remote_destination,tmp),mode=to_mode,camera=camera,postfix=postfix)
		ENDCASE

		; If destination is a remote directory scp the files
		; and delete the local copies

		tmp = 'scp '+workfile+' '+filepath(root=tmp,workname)
		print, tmp
		spawn, tmp
		tmp = do_file(workfile, /delete)

	ENDIF

	stop_time = TimeSystem(/silent)

	IF silent LE 2 THEN BEGIN
		tmp = TimeOp(/subtract,stop_time,start_time,TimeUnit(/min))
		tmp = string(tmp,format='(F10.2)')
		message, /info, skyhide+' '+strtrim(tmp,2)+' min'
		print
		print
	ENDIF

ENDFOR

RETURN  &  END