;+
; NAME:
;	smei_star_remove
; PURPOSE:
;	Removes bright stars from SMEI orbital sky maps
; CATEGORY:
;	ucsd/camera/idl/star
; CALLING SEQUENCE:
	PRO smei_star_remove, wanted_map , $

		destination = destination	, $
		to_smeidb	= to_smeidb		, $

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

		keepbkgnd	= keepbkgnd 	, $
		removezld 	= removezld		, $

		cleanedge	= cleanedge 	, $

		version 	= version		, $
		get_version	= get_version	, $
		reverse_loop= reverse_loop	, $

		_extra		= _extra
; INPUTS:
;	wanted_map		array[n]; type: string, time structure or integer
;						passed to href=smei_getfile=
;						Selects skymap files to be processed
; 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 	scalar; type: string; default: $TUB
;					array[2]; type: string
;						destination directory for skymaps and time series file
;						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 two directories are specified the first is used
;						for the skymap and the second for the time series
;						If destination includes a "user@host:" identifier then
;						the output is created locally in $TUB
;						and is then transferred to the remote host using scp;
;						the local files are 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.
;	/overwrite		forces overwriting of existing skymaps and timeseries files
;					(by default no files are overwritten)
;	/checkversion	if present, files will be overwritten only if the current
;					star subtraction version number and the current sidereal
;					background version is higher then the version in the Fits
;					header of the skymap.
;					(if /keepbkgnd is set then the background version is not
;					checked)
;
;	silent=silent 	if set, suppresses informational messages
;	upto_mags=upto_mags
;		 			scalar; type: float; default: 6.0
;						maximum SMEI magnitude; only stars brighter than this
;						are removed. (passed to href=smei_star_fit)
;	/keepbkgnd		if set, the sidereal background (Milky Way, faint stars,
;						etc.) is not subtracted from the skymaps after
;						the bright stars are subtracted.
;
;	/cleanedge		if set, a correction is made for straylight from sky
;						just outside the fov into the outer areas (direction
;						cosine angles larger than 20 deg) of the camera fov.
;	_extra=_extra	The _extra keyword is passed to href=smei_star_fit=,
;						and href=smei_getfile=, so check those for
;						additional input keywords.
;
;	/get_version	retrieve software version number
;						(returned in keyword "version")
; OUTPUTS:
;	Files in destination directory
; OPTIONAL OUTPUTS:
;	version=version Only defined if /get_version is set
;					array[2]; type: float
;						version[0]: version of bright star subtraction
;						version[1]: version of sidereal background map
; INCLUDE:
	@compile_opt.pro		; On error, return to caller
; CALLS:
;	InitVar, CheckDir, IsType, destroyvar
;	smei_filename, smei_filepath, smei_getfile, smei_star_fit
;	smei_camera_gain, smei_star_writepnt, smei_sky_cleanedge_fov
;	GetFileSpec, TimeSystem, TimeUnit, gzip_file
;	hide_env
; EXAMPLE:
;	smei_star_remove, 1500
;		Subtract stars for orbit 1500 (specified as integer) for camera 1
;	smei_star_remove, smei_coriolis([1500,1501]), camera=1
;		Subtract stars for orbits 1500 and 1501 (specified as time
;		structure) for camera 1
;	smei_star_remove, '2003/04/16 18:25:57'
;		Subtract stars for orbit 1500 specified as time string
;	smei_star_remove, '*.fts.gz'
;		Use dialog_pickfile to select single skymap
;	smei_star_remove, 'c1sky_2006_045_021031.fts.gz'
;		Subtract stars for specified map
; PROCEDURE:
;	For each skymap a new skymap (Fits file) with the star-subtracted
;	maps is created, and an ASCII file with information about the
;	fitted stars. Both files are created in the directories
;	specified in keyword 'destination'.
;
;	The current star subtraction is build around the following
;	set of keywords:
;
;	smei_star_remove, f, /cleanedge, /use_weights, /auto_wing
;
;	This implies a wind_radius value of 1.4 degrees.
;	A good way to redo a section of maps in a time range 'trange' in the SMEI
;	skymap data base on $SMEISKY0 is the command:
;	
;	smei_star_remove, trange,cam=3,/cleanedge,/use_weights,/auto_wing,	$
;		/checkversion,silent=3,/to_smeidb
;
; MODIFICATION HISTORY:
;	NOV-2005, Paul Hick (UCSD/CASS)
;	JUN-2006, Paul Hick (UCSD/CASS)
;		Modified to accept arrays as input in wanted_map argument.
;	JUN-2006, Paul Hick (UCSD/CASS)
;		Added Fits keyword VERSION to sky map
;		Added NAME, CREATED and VERSION keywords to time series file
;		Added keyword /overwrite
;	JUL-2006, Jordan Vaughan, Paul Hick (UCSD/CASS)	V2.0
;		Added subtraction of sidereal background map (can be disabled
;		by keyword /keepbkgnd)
;
;		Added two extensions (copied from original skymap):
;		-extension 6: fraction of orbit (extension 7 from original skymap)
;		-extension 7: time in sec since first frame used
;				(extension 14 from original skymap)
;
;		Keyword VERSION is now called BSTARVER (to distinguish it from
;		keyword BKGNDVER (the version of the sidereal background map)
;		Keyword MILKYWAY is now called BKGND
;		If the sidereal backgound is subtracted keyword BKGND is set
;		to 1, and two more keywords are added:
;		BKGNDVER is the version of the sidereal background map
;		BKGAIN   is the gain correction (multiplicative factor) applied
;		to the sidereal background map before subtracting it.
;
;		Fixed problem with header for extension 4 (map of removed stars
;		in equatorial map.
;	AUG-2006, Paul Hick (UCSD/CASS)	V3.0
;		Modifications related to separate fitting of background and PSF.
;		Added couple of entries to smei_star_fit structure.
;		Time series file is now written by procedure smei_star_writepnt
;	AUG-2006, Paul Hick (UCSD/CASS)	V3.0
;		Padded records in time series file with 5 spaces to make them
;		exactly the same length as the titles record.
;	SEP-2006, Paul Hick (UCSD/CASS) V3.1
;		Added /use_filled_psf keyword
;		Added star_info keyword to smei_star_fit and smei_star_writepnt
;	NOV-2006, Paul Hick (UCSD/CASS) V4.0
;		Rearranged layout of time series file
;		Added extension (lowres map copied from original skymap)
;		-extension 8: number of pixels per skybin
;				(extension 16 from original skymap)
;	MAR-2007, Paul Hick (UCSD/CASS) V4.1
;		Added extension (lowres map copied from original skymap)
;		-extension 8: direction cosine angle
;				(extension 15 from original skymap)
;		-extension 9: old extension 8 (pixels per skybin)
;	APR-2007, Paul Hick (UCSD/CASS) V4.2
;		Added cleanedge keyword to allow removal of straylight
;		from sky just outside the camera fov.
;	MAY-2007, Paul Hick (UCSD/CASS) V4.3
;		Modified fovangle field and added dfovangle field
;		to star_fit structure
;	MAY-2007, Paul Hick (UCSD/CASS) V4.31
;		If keywords /cleanedge is set and /keepbkgnd NOT set,
;		then the sidereal background is added again AFTER
;		the star subtraction, and PRIOR to the removal of
;		the scattered light along the edge of the FOV. After
;		the scattered light correction the sidereal background
;		is taken out again.
;		Added keyword /to_smeidb to write directly into SMEI data base
;		Added keyword remote to do transfers to remote host
;	JUN-2007, Paul Hick (UCSD/CASS)
;		Added keyword /from_smeidb
;	JUL-2007, Paul Hick (UCSD/CASS)
;		Bug fix. removebkgnd was written into Fits file as byte
;		instead of short int.
;	JUL-2007, Paul Hick (UCSD/CASS), V5.01
;		Added three extensions to output Fits file.
;		These contain 3 skymaps (equatorial,north and south pole)
;		with the stars subtracted and the core PSF replaced by
;		the background model.
;		Modified use of area defined with wing_radius. Bins in this
;		area that are part of the core PSF of a brighter neighbour
;		star (which has been subtracted already) are excluded from
;		being used in the lsq fit. They are still used when the
;		star model is subtracted. 
;		Removed keywords wing_radius and use_filled_psf
;		(now passed to smei_star_fit using _extra keyword).
;	NOV-2007, Paul Hick (UCSD/CASS), V5.02
;		Fixes bug in calculation of camera gain for cameras 1 and 3
;		smei_camera_gain.
;	DEC-2007, Paul Hick (UCSD/CASS), V5.03
;		Added /inside_range keyword to be able to process
;		files in a specified time range (specified as a 2-element
;		time array in 'wanted_map')
;		Changed keyword CREATED in time series file to DATE.
;		Shortened header for equatorial maps in ext 3 and 6
;		(TX1MOON, TX2MOON, TX3MOON removed).
;	DEC-2007, Paul Hick (UCSD/CASS), V5.04
;		Added two lines to header for timeseries files
;		(clean edge and remove bkgnd keywords)
;	DEC-2007, Paul Hick (UCSD/CASS), V5.05
;		Added some more keyword to timeseries header:
;		BSTAR,BSTARVER,BKGNDVER,BKGAIN,CLNEDGE
;		(these are also present in the Fits file with the
;		star-subtracted skymaps).
;	DEC-2007, Paul Hick (UCSD/CASS), V5.06
;		Added keywords ONEORBIT and ORBIT to timeseries header.
;	JAN-2008, Paul Hick (UCSD/CASS), V5.07
;		Fixed some format issues in pnt file.
;		Star name was truncated at 11 char, should be 12.
;		Planet names were print right-adjusted, should be
;		left-adjusted
;	FEB-2008, Paul Hick (UCSD/CASS)
;		Location of files now done using smei_getfile.
;		Remove keyword remote (functionally is now provided by
;		keyword destination).
;	FEB-2008, Paul Hick (UCSD/CASS), V5.10
;		Added star_fit.remove entry.
;		Bug fix in calculation of psfwidth (sometimes gave NaN).
;	FEB-2008, Paul Hick (UCSD/CASS), V5.11
;		Added USNO asteroids and 'stars_not_subtracted' catalogue to
;		list of point sources (these are fitted, but not subtracted).
;	MAR-2008, Paul Hick (UCSD/CASS), V5.12
;		Added keyword SMEI_SKY, SMEI_CAL and SMEI_ORB to pnt header
;	MAR-2008, Paul Hick (UCSD/CASS), V5.13
;		Fixed problem with two DATE keywords in header of output
;		Fits file.
;	MAR-2008, Paul Hick (UCSD/CASS), V5.2
;		Bug fix in edge-clean procedure. The bug may have been
;		introduced with V5.01 (JUL-2007) with the inclusion of the
;		'flat' star-subtracted maps.
;	APR-2008, Paul Hick (UCSD/CASS), V5.3
;		Both sidereal background and zodiacal dust brightness
;		are now subtracted prior to fitting stars.
;		The weighting in href=smei_star_lsqfit= was changed to 
;		from 1/b to 1/(50+b), where b is the residual brightness
;		(after background and zodiac subtraction).
;		Added keyword removezld.
;		Added header keywords ZLD, ZLDVERSION to output files.
;		Renamed header keyword gain from BKGAIN to SKYGAIN
;	APR-2008, Paul Hick (UCSD/CASS), V5.31
;		Replaced sidereal background map with new version.
;	JUN-2008, Paul Hick (UCSD/CASS)
;		Added PSF angle to output Fits file.
;	JUL-2008, Paul Hick (UCSD/CASS), V5.32
;		Changed gain decrease for camera 1 from 1%/yr to 0.5%/yr
;		(in smei_camera_gain)
;	OCT-2008, Paul Hick (UCSD/CASS; pphick@ucsd.edu), V5.33
;		Fixed bug in calculation of observation time (used stime as
;		time origin instead of extracted the origin from the header
;		for the time map (keyword TORIGIN).
;	DEC-2008, Paul Hick (UCSD/CASS; pphick@ucsd.edu), V6.00
;		New sidereal background map (version 3)
;		New zodiacal dust cloud model (model 4.01)
;		Gain change for camera 1 set back to -1% per year.
;-

bstar_version = 6.00 				; Version nr for bright star subtraction

InitVar, to_smeidb  , /key 			; Write to $SMEISKY0
InitVar, destination, ([getenv('TUB'),'SMEIDB?'])[to_smeidb]

InitVar, silent 	 , 0
InitVar, overwrite	 , /key
InitVar, checkversion, /key
InitVar, keepbkgnd	 , /key
InitVar, removezld	 , /key
InitVar, cleanedge	 , /key

InitVar, get_version , /key
InitVar, reverse_loop, /key

; Make sure the sidereal background FITS file exists
; If it exists, read the equatorial, north and south polar background maps

tmp = filepath(root=smei_filepath(mode='sid'),subdir='v3','c0sid*.fts.gz')
bkgnd_map = (file_search(tmp))[0]
have_bkgnd = bkgnd_map NE ''		; Background maps exist

IF get_version THEN BEGIN
	version = [bstar_version,-1.00]
	IF have_bkgnd THEN	$		; Sidereal background version
		version[1] = fxpar(headfits(bkgnd_map,/silent),'BKGNDVER')
	RETURN
ENDIF

from_mode = 'sky'
to_mode   = [ 'equ', 'pnt']
to_type   = ['.fts','.txt']

; Need two destinations:
; for skymap and time series file

InitVar, destination , count=1, [destination,destination]
to_remote = strpos(destination,':') NE -1

; Set up local and remote destination

ndestination = n_elements(destination)

FOR i=0,ndestination-1 DO BEGIN
	CASE to_remote[i] OF
	0: BEGIN
		boost, local_destination , destination[i]
		boost, remote_destination, ''
	END
	1: BEGIN
		boost, local_destination , getenv('TUB')
		boost, remote_destination, destination[i]
	END
	ENDCASE
ENDFOR

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

IF count EQ 0 THEN	$
	RETURN


CASE have_bkgnd OF

0: message, /info, 'no sidereal background, '+hide_env(tmp)

1: BEGIN
	message, /info, 'sidereal background, '+hide_env(bkgnd_map)

	eq_bkgnd0 = readfits(bkgnd_map, hdr, exten_no=0, /silent); Equatorial map

	bkgnd_version = fxpar(hdr, 'BKGNDVER') 					; Version number of sidereal background map

	eq_cra  = fxpar(hdr,'CRPIX1')-1
	eq_cdec = fxpar(hdr,'CRPIX2')-1
	eq_bdeg = 1.0/fxpar(hdr,'CDELT2')

	np_bkgnd0 = readfits(bkgnd_map, hdr, exten_no=1, /silent); North pole map
	np_cra  = fxpar(hdr,'CRPIX1')-1
	np_cdec = fxpar(hdr,'CRPIX2')-1
	np_bdeg = 1.0/fxpar(hdr,'CDELT2')

	sp_bkgnd0 = readfits(bkgnd_map, hdr, exten_no=2, /silent); South pole map
	sp_cra  = fxpar(hdr,'CRPIX1')-1
	sp_cdec = fxpar(hdr,'CRPIX2')-1
	sp_bdeg = 1.0/fxpar(hdr,'CDELT2')
END

ENDCASE

removebkgnd = have_bkgnd AND (1-keepbkgnd)

IF NOT removebkgnd THEN		$
	IF silent LE 2 THEN		$
		message, /info,'not subtracting sidereal background'

stop_time = TimeSystem(/silent)

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

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)

	destroyvar, workname, workfile

	bad_dir = 0

	FOR i=0,ndestination-1 DO BEGIN

		boost, workname, smei_filename(skytime,camera=camera,mode=to_mode[i],type=to_type[i])

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

		IF NOT CheckDir(tmp) THEN bad_dir += 1
			
		boost, workfile, filepath(root=tmp,workname[i])

	ENDFOR

	IF bad_dir NE 0 THEN continue

	IF checkversion THEN BEGIN

		tmp = (file_search(workfile[0]+'*'))[0]
		IF tmp NE '' THEN BEGIN 	; Skymap exists
			hdr = headfits(tmp)
			map_bstar_version = fxpar(hdr,'BSTARVER')
			has_bstar_version = !err EQ 0

			; CONTINUE not allowed inside CASE block

			IF have_bkgnd THEN BEGIN
				map_bkgnd_version = fxpar(hdr,'BKGNDVER')
				has_bkgnd_version = !err EQ 0

				IF (has_bstar_version AND map_bstar_version GE bstar_version) AND $
				   (has_bkgnd_version AND map_bkgnd_version GE bkgnd_version) THEN BEGIN
					IF silent LE 2 THEN message, /info, '>= version, '+hide_env(workfile[0])
					CONTINUE
				ENDIF
			ENDIF ELSE BEGIN
				IF has_bstar_version AND map_bstar_version GE bstar_version THEN BEGIN
					IF silent LE 2 THEN message, /info, '>= version, '+hide_env(workfile[0])
					CONTINUE
				ENDIF
			ENDELSE
		ENDIF

	ENDIF ELSE IF NOT overwrite THEN BEGIN

		tmp = (file_search(workfile[0]+'*'))[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

	; Read maps and headers; then remove stars

	smei_sky,skyfile,/exists,/equatraw,sky=eq_raw,hdr=eq_struct,/noplot,silent=silent+1,nbin=1, exten_no=eq_exten
	eq_exists = size(eq_raw,/n_dim) NE 0
	eq_full   = n_elements(eq_raw) EQ 3600L*1200L
	IF eq_exists THEN eq_good = where(finite(eq_raw)) ELSE eq_good = -1

	smei_sky,skyfile,/exists,/northraw,sky=np_raw,hdr=np_struct,/noplot,silent=silent+1,nbin=1, exten_no=np_exten
	np_exists = size(np_raw,/n_dim) NE 0
	np_full   = n_elements(np_raw) EQ 800L*800L
	IF np_exists THEN np_good = where(finite(np_raw)) ELSE np_good = -1

	smei_sky,skyfile,/exists,/southraw,sky=sp_raw,hdr=sp_struct,/noplot,silent=silent+1,nbin=1, exten_no=sp_exten
	sp_exists = size(sp_raw,/n_dim) NE 0
	sp_full   = n_elements(sp_raw) EQ 800L*800L
	IF sp_exists THEN sp_good = where(finite(sp_raw)) ELSE sp_good = -1

	; The sidereal background is removed before subtracting stars.
	; Presumably this makes the background near the galactic plane
	; flatter (removing first and higher spatial derivatives), and thus
	; improving the star fitting.

	skygain = smei_camera_gain(skytime, camera=camera)

;view, in=rebin(eq_raw < 500,900,300),/stretch

	; Calculate the sidereal background and the zodiacal light
	; distribution for all good bins

	IF NOT (eq_full AND np_full AND sp_full) THEN message, 'not yet implemented'

	; Set the weights to be used in the linear fit in smei_star_lsqfit
	; from the skymaps prior to subtracting sidereal background and
	; zodiacal light distribution. This makes the weights (and hence the
	; result of the fit) largely independent of the background.

	; The weights will only be used (in smei_star_lsqfit) if /use_weights is set.
	; (For now we are depending on the default in smei_star_lsqfit of 1/(50+brightness)

	;eq_weights = 1.0/eq_raw
	;np_weights = 1.0/np_raw
	;sp_weights = 1.0/sp_raw

	IF eq_good[0] NE -1 THEN BEGIN
		IF have_bkgnd THEN eq_bkgnd = eq_bkgnd0[eq_good]*skygain ELSE eq_bkgnd = 0
		tmp = ArrayLocation(eq_good,dim=size(eq_raw,/dim))
		tmp = cvsmei(from_map=tmp,mode='eqmap',/to_equatorial,/degrees,/silent)
		eq_zld = smei_zld_model(skytime,tmp,/equatorial,/degrees,/silent,version=zld_version)*skygain
	ENDIF

	IF np_good[0] NE -1 THEN BEGIN
		IF have_bkgnd THEN np_bkgnd = np_bkgnd0[np_good]*skygain ELSE np_bkgnd = 0
		tmp = ArrayLocation(np_good,dim=size(np_raw,/dim))
		tmp = cvsmei(from_map=tmp,mode='npmap',/to_equatorial,/degrees,/silent)
		np_zld = smei_zld_model(skytime,tmp,/equatorial,/degrees,/silent,version=zld_version)*skygain
	ENDIF

	IF sp_good[0] NE -1 THEN BEGIN
		IF have_bkgnd THEN sp_bkgnd = sp_bkgnd0[sp_good]*skygain ELSE sp_bkgnd = 0
		tmp = ArrayLocation(sp_good,dim=size(sp_raw,/dim))
		tmp = cvsmei(from_map=tmp,mode='spmap',/to_equatorial,/degrees,/silent)
		sp_zld = smei_zld_model(skytime,tmp,/equatorial,/degrees,/silent,version=zld_version)*skygain
	ENDIF

	; Remove sidereal background and zodiacal light prior to star subtraction

	IF silent LE 2 THEN	$
		message, /info, 'remove sidereal background and zodiacal light prior to star subtraction'

	IF eq_good[0] NE -1 THEN eq_raw[eq_good] -= eq_bkgnd+eq_zld	; Equator
	IF np_good[0] NE -1 THEN np_raw[np_good] -= np_bkgnd+np_zld	; North pole
	IF sp_good[0] NE -1 THEN sp_raw[sp_good] -= sp_bkgnd+sp_zld	; South pole

;view, in=rebin(eq_raw < 500,900,300),/new,/stretch

	; Fit and subtract the stars: *_raw -> *_sky

	smei_star_fit, hdr=eq_struct, sky=eq_raw, final_sky=eq_sky	, $
		flat_sky=eq_flat, star_sky=eq_stars, star_fit=eq_fit	, $
		count=eq_count, /degrees, badfit=0.0, silent=silent 	, $
		_extra=_extra, star_info=star_info, weights=eq_weights
	print
	print

	smei_star_fit, hdr=np_struct, sky=np_raw, final_sky=np_sky	, $
		flat_sky=np_flat, star_sky=np_stars, star_fit=np_fit	, $
		count=np_count, /degrees, badfit=0.0, silent=silent 	, $
		_extra=_extra, /northpole, weights=np_weights
	print
	print

	smei_star_fit, hdr=sp_struct, sky=sp_raw, final_sky=sp_sky	, $
		flat_sky=sp_flat, star_sky=sp_stars, star_fit=sp_fit	, $
		count=sp_count, /degrees, badfit=0.0, silent=silent 	, $
		_extra=_extra, /southpole, weights=sp_weights
	print
	print

	IF silent LE 2 THEN BEGIN
		tmp = [eq_count,np_count,sp_count]
		message, /info, strjoin(strcompress(tmp,/rem),'+')+'='	+ $
			strcompress(round(total(tmp)),/rem)+' stars, '+skyhide
	ENDIF

	IF cleanedge AND NOT (eq_full AND np_full AND sp_full) THEN BEGIN
		cleanedge = 0
		message, /info, '/cleanedge ignored: not a full-skymap, '+skyhide
	ENDIF

;view, in=rebin(eq_sky > (-25) < 25,900,300),/new,/stretch
;view, in=rebin(eq_flat > (-25) < 25,900,300),/new,/stretch
;return

	; At this point the sidereal background, the zodiacal light
	; and all the stars have been removed.

	CASE cleanedge OF

	0: BEGIN

		; If necessary, add sidereal background and/or
		; zodiacal dust brightness back in.

		IF eq_good[0] NE -1 THEN BEGIN
			tmp = (1-removebkgnd)*eq_bkgnd+(1-removezld)*eq_zld
			eq_sky [eq_good] += tmp	; Equator
			eq_flat[eq_good] += tmp	; Equator
		ENDIF

		IF np_good[0] NE -1 THEN BEGIN
			tmp = (1-removebkgnd)*np_bkgnd+(1-removezld)*np_zld
			np_sky [np_good] += tmp	; North pole
			np_flat[np_good] += tmp	; North pole
		ENDIF

		IF sp_good[0] NE -1 THEN BEGIN
			tmp = (1-removebkgnd)*sp_bkgnd+(1-removezld)*sp_zld
			sp_sky [sp_good] += tmp	; South pole
			sp_flat[sp_good] += tmp	; South pole
		ENDIF

	END

	1: BEGIN

		IF silent LE 2 THEN 	$
			message, /info, 'scattered light at edge of FOV'

		; Add the sidereal background and zodiacal light back in
		; (they contribute to the scattered light removed by smei_sky_cleanedge).

		IF eq_good[0] NE -1 THEN BEGIN
			tmp = eq_bkgnd+eq_zld
			eq_sky [eq_good] += tmp		; Equator
			eq_flat[eq_good] += tmp		; Equator
		ENDIF

		IF np_good[0] NE -1 THEN BEGIN
			tmp = np_bkgnd+np_zld
			np_sky [np_good] += tmp		; North pole
			np_flat[np_good] += tmp		; North pole
		ENDIF

		IF sp_good[0] NE -1 THEN BEGIN
			tmp = sp_bkgnd+sp_zld
			sp_sky [sp_good] += tmp		; South pole
			sp_flat[sp_good] += tmp		; South pole
		ENDIF

		; Get the low res data needed for the scattered light removal

		smei_sky, skyfile, /exists, /fovx   , sky=fovx_map, /noplot, silent=silent+1, /degrees
		smei_sky, skyfile, /exists, /orbsecs, sky=time_map, /noplot, silent=silent+1

		; Remove scattered light from just outside fov
		; Concatenate eq_sky and eq_flat.
		; eq_sky is used to calculate the scattered light
		; The same is subtracted from both eq_sky and eq_flat
		; (same for np and sp maps)

		eq_sky = reform([[eq_sky],[eq_flat]], [size(eq_sky,/dim),2])
		np_sky = reform([[np_sky],[np_flat]], [size(np_sky,/dim),2])
		sp_sky = reform([[sp_sky],[sp_flat]], [size(sp_sky,/dim),2])

		smei_sky_cleanedge_fov	, $
			eq_sky	= eq_sky	, $
			np_sky	= np_sky	, $
			sp_sky	= sp_sky	, $

			time_map= time_map	, $
			fovx_map= fovx_map	, $

			/degrees			, $
			silent	= silent

		eq_flat = eq_sky[*,*,1]
		eq_sky  = eq_sky[*,*,0]

		np_flat = np_sky[*,*,1]
		np_sky  = np_sky[*,*,0]

		sp_flat = sp_sky[*,*,1]
		sp_sky  = sp_sky[*,*,0]

		; If necessary, remove sidereal background and/or
		; zodiacal dust brightness back in.

		IF eq_good[0] NE -1 THEN BEGIN
			tmp = removebkgnd*eq_bkgnd+removezld*eq_zld
			eq_sky [eq_good] -= tmp	; Equator
			eq_flat[eq_good] -= tmp	; Equator
		ENDIF

		IF np_good[0] NE -1 THEN BEGIN
			tmp = removebkgnd*np_bkgnd+removezld*np_zld
			np_sky [np_good] -= tmp	; North pole
			np_flat[np_good] -= tmp	; North pole
		ENDIF

		IF sp_good[0] NE -1 THEN BEGIN
			tmp = removebkgnd*sp_bkgnd+removezld*sp_zld
			sp_sky [sp_good] -= tmp	; South pole
			sp_flat[sp_good] -= tmp	; South pole
		ENDIF

	END

	ENDCASE

;view, in=rebin(eq_sky < 500,900,300),/new,/stretch
;return

	; ===================================================
	; Start writing the Fits file for the subtracted maps.
	; hdr is the primary header from the original skymap

	; Mark bad data with BAD_DATA value from Fits header

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

	sky_version = fxpar(hdr,'SMEI_SKY')
	IF !err LT 0 THEN sky_version = 0.00
	cal_version = fxpar(hdr,'SMEI_CAL')
	IF !err LT 0 THEN cal_version = 0.00
	orb_version = fxpar(hdr,'SMEI_ORB')
	IF !err LT 0 THEN orb_version = 0.00

	IF eq_exists THEN BEGIN
		bad = where(1-finite(eq_sky))
		IF bad[0] NE -1 THEN eq_sky  [bad] = bad_data

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

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

	IF np_exists THEN BEGIN
		bad = where(1-finite(np_sky))
		IF bad[0] NE -1 THEN np_sky  [bad] = bad_data

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

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

	IF sp_exists THEN BEGIN
		bad = where(1-finite(sp_sky))
		IF bad[0] NE -1 THEN sp_sky  [bad] = bad_data

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

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

	; Set up Fits header template

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

	mkhdr, main_hdr, IsType(eq_sky), size(eq_sky,/dim), /extend
	tmp = (where(strpos(hdr,'SMEI_SKY') EQ 0))[0]
	IF tmp EQ -1 THEN tmp = (where(strpos(hdr,'SMEI_HTM') EQ 0))[0]

	; This cuts of the IDL Fits header after the EXTEND keyword,
	; dropping the DATE keyword and a long COMMENT line.
	; (there is another DATE keyword copied from the primary header
	; of the input skymap, which is updated below).

	main_hdr = [main_hdr[0:(where(strpos(main_hdr,'DATE    =') EQ 0))[0]-1],hdr[tmp:*]]

	; Main header updates

	fxaddpar, main_hdr, 'NAME'	  , GetFileSpec(workfile[0],part='name')
	fxaddpar, main_hdr, 'IMG_TYPE', 'skymap_'+to_mode[0]
	fxaddpar, main_hdr, 'DATE'	  , TimeGet(TimeSystem(/silent),format='YYYY-MN-DD hh:mm',/scalar)

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

	fxaddpar, main_hdr, 'SKYGAIN' , skygain			, ' Gain factor', after='BAD_DATA'

	IF removezld THEN	$
		fxaddpar, main_hdr, 'ZLDVER'  , zld_version	, ' Zodiacal light model'  , after='BAD_DATA'

	; Don't write byte into Fits file!

	fxaddpar, main_hdr, 'ZLD' 		, fix(removezld), ' Zodiacal light removed (0=N,1=Y)'  , after='BAD_DATA'

	IF removebkgnd THEN	$
		fxaddpar, main_hdr, 'BKGNDVER', bkgnd_version, ' Sidereal background version'	  , after='BAD_DATA', format='(F5.2)'

	; Don't write byte into Fits file!

	fxaddpar, main_hdr, 'BKGND' 	, fix(removebkgnd)	, ' Sidereal background removed (0=N,1=Y)' , after='BAD_DATA'

	fxaddpar, main_hdr, 'BSTARVER'	, bstar_version , ' Star subtraction version'		  , after='BAD_DATA', format='(F5.2)'
	fxaddpar, main_hdr, 'BSTAR' 	, 1				, ' Bright stars subtracted (0=N,1=Y)', after='BAD_DATA'

	; Basic header for equatorial maps

	mkhdr, eq_hdr, IsType(eq_sky), size(eq_sky,/dim), /image
	eq_hdr = [eq_hdr[0:(where(strpos(eq_hdr,'END   ') EQ 0))[0]-1],main_hdr[(where(strpos(main_hdr,'MAP') EQ 0))[0]:(where(strpos(main_hdr,'BAD_DATA') EQ 0))[0]]]

	; Basic header for maps of north pole

	hdr = headfits(skyfile, silent=silent, exten=np_exten)
	mkhdr, np_hdr, IsType(np_sky), size(np_sky,/dim), /image
	np_hdr = [np_hdr[0:(where(strpos(np_hdr,'END   ') EQ 0))[0]-1],hdr[(where(strpos(hdr,'MAP') EQ 0))[0]:*]]

	; Basic header for maps of south pole

	hdr = headfits(skyfile, silent=silent, exten=sp_exten)
	mkhdr, sp_hdr, IsType(sp_sky), size(sp_sky,/dim), /image
	sp_hdr = [sp_hdr[0:(where(strpos(sp_hdr,'END   ') EQ 0))[0]-1],hdr[(where(strpos(hdr,'MAP') EQ 0))[0]:*]]

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

	; Extension 0, equatorial map with stars removed

	hdr = main_hdr
	fxaddpar, hdr, 'MAP', 'Equatorial map RA=[0,360],DEC=[-60,+60]', ' Map description'
	fxaddpar, hdr, 'CLNEDGE', (['F','T'])[cleanedge], ' Clean edge'
	writefits, workfile[0], eq_sky, hdr

	; Extension 1, map of north pole with stars removed

	hdr = np_hdr
	fxaddpar, hdr, 'MAP', 'Map of north pole DEC=[+50,+90]', ' Map description'
	writefits, workfile[0], np_sky, hdr, /append

	; Extension 2, map of south pole with stars removed

	hdr = sp_hdr
	fxaddpar, hdr, 'MAP', 'Map of south pole DEC=[-50,-90]', ' Map description'
	writefits, workfile[0], sp_sky, hdr, /append

	; Extension 3, equatorial map with stars removed with background in PSF area

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

	; Extension 4, map of north pole with stars removed with background in PSF area

	hdr = np_hdr
	fxaddpar, hdr, 'MAP', 'Map of north pole DEC=[+50,+90] (flat)', ' Map description'
	writefits, workfile[0], np_flat, hdr, /append

	; Extension 5, map of south pole with stars removed with background in PSF area

	hdr = sp_hdr
	fxaddpar, hdr, 'MAP', 'Map of south pole DEC=[-50,-90] (flat)', ' Map description'
	writefits, workfile[0], sp_flat, hdr, /append

	; Extension 6, equatorial map of removed stars

	hdr = eq_hdr
	fxaddpar, hdr, 'MAP'   , 'Stars removed from equator RA=[0,360],DEC=[-60,+60]', ' Map description'
	fxaddpar, hdr, 'NSTARS', round(eq_count), ' Number of stars removed'
	writefits, workfile[0], eq_stars, hdr, /append

	; Extension 7, map of stars removed from north pole

	hdr = np_hdr
	fxaddpar, hdr, 'MAP'   , 'Stars removed from north pole DEC=[+50,+90]', ' Map description'
	fxaddpar, hdr, 'NSTARS', round(np_count), ' Number of stars removed'
	writefits, workfile[0], np_stars, hdr, /append

	; Extension 8, map of stars removed from south pole

	hdr = sp_hdr
	fxaddpar, hdr, 'MAP'   , 'Stars removed from south pole DEC=[-50,-90]', ' Map description'
	fxaddpar, hdr, 'NSTARS', round(sp_count), ' Number of stars removed'
	writefits, workfile[0], sp_stars, hdr, /append

	; Copy two time extensions from original skymap

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

	tmp = readfits(skyfile, hdr, exten_no=5+2, /silent)
	writefits, workfile[0], tmp, hdr, /append

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

	tmp = readfits(skyfile, hdr, exten_no=5+9, /silent)
	writefits, workfile[0], tmp, hdr, /append

	; Extension 11, Direction cosine angle-x RA=[0,360],DEC=[-90,+90]

	tmp = readfits(skyfile, hdr, exten_no=5+10, /silent)
	writefits, workfile[0], tmp, hdr, /append

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

	tmp = readfits(skyfile, hdr, exten_no=5+11, /silent)
	writefits, workfile[0], tmp, hdr, /append

	; Extension 13, PSF position angle from north RA=[0,360],DEC=[-90,+90]

	tmp = readfits(skyfile, hdr, exten_no=5+3, /silent)
	writefits, workfile[0], tmp, hdr, /append

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

	tmp = gzip_file(workfile[0], /force)
	workfile[0] += '.gz'
	IF silent LE 3 THEN message, /info, hide_env(workfile[0])+' is '+to_mode[0]+' map'

	dec_cutoff = 55.0			; Same units as RA,dec in fit_* structure (degrees)

	; Write the .pnt text file.
	; First write the header only

	IF IsType(eq_fit,/defined) OR IsType(np_fit,/defined) OR IsType(sp_fit,/defined) THEN BEGIN

		IF silent LE 3 THEN $
			message, /info, hide_env(workfile[1])+' is '+to_mode[1]+' file'

		boost, star_info,		$
			['SMEI_SKY: '+string(sky_version  ,format='(F5.2)')	, $
			 'SMEI_CAL: '+string(cal_version  ,format='(F5.2)')	, $
			 'SMEI_ORB: '+string(orb_version  ,format='(F5.2)')	, $
			 'BSTAR   : yes'									, $
			 'BSTARVER: '+string(bstar_version,format='(F5.2)')	]

		boost, star_info, 'BKGND   : '+(['no','yes'])[removebkgnd]
		IF removebkgnd THEN		$
			boost, star_info, 'BKGNDVER: '+string(bkgnd_version,format='(F5.2)')

		boost, star_info, 'ZLD     : '+(['no','yes'])[removezld]
		IF removezld THEN		$
			boost, star_info, 'ZLDVER: '+string(zld_version,format='(F5.2)')

		boost, star_info, 'SKYGAIN : '+strcompress(skygain,/remove_all)

		boost, star_info, 'CLNEDGE : '+(['no','yes'])[cleanedge]

		oneorbit = fxpar(main_hdr, 'ONEORBIT')
		boost, star_info, 'ONEORBIT: '+(['no','yes'])[oneorbit]

		IF oneorbit THEN BEGIN
			orbitnum = fxpar(main_hdr,'ORBIT')
			boost, star_info, 'ORBIT   : '+strcompress(orbitnum, /rem)
		ENDIF

		smei_star_writepnt, workfile[1] , $
			dec_cutoff	= dec_cutoff, $
			camera		= camera	, $
			star_info	= star_info	, $
			/degrees

	ENDIF ELSE 	$
		IF silent LE 3 THEN message, /info, 'no time series file created'

	; Write the structures (nothing happens if the structure doesn't exist)

	smei_star_writepnt, workfile[1], eq_fit, dec_cutoff=dec_cutoff, /append, /degrees
	smei_star_writepnt, workfile[1], np_fit, dec_cutoff=dec_cutoff, /append, /degrees
	smei_star_writepnt, workfile[1], sp_fit, dec_cutoff=dec_cutoff, /append, /degrees

	FOR i=0,ndestination-1 DO BEGIN

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

		IF to_remote[i] THEN BEGIN

			tmp = strpos(remote_destination[i],':')+1

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

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

		ENDIF

	ENDFOR

	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