;+
; NAME:
;	smei_findpnt
; PURPOSE:
;	Originally created to look for point sources in the star-subtracted
;	equ files to find unknown objects (i.e. objects not included in the
;	SMEI star catalogues).
;
;	To run this program access is needed to the 'equ' files, If keyword
;	/use_catalogue is set then also access to the 'sky' files is required.
; CATEGORY:
;	smei/ucsd/camera/idl/star
; CALLING SEQUENCE:
	PRO smei_findpnt, wanted_map 	, $
		camera		= camera		, $
		from_mode	= from_mode 	, $
		to_mode 	= to_mode		, $
		destination	= destination	, $
		fix_centroid= fix_centroid	, $
		cvmin		= cvmin			, $
		use_catalogue= use_catalogue, $
		npeak		= npeak			, $
		silent		= silent		, $
		_extra		= _extra
; INPUTS:
;	wanted_map		array[n]; type: string, time structure or integer
;						passed to href=smei_getfile=.
;						Selects skymap files to be processed.
;						Usually this is a pair of times (start time and
;						stop time) to process all maps inside a period
;						of time.
; OPTIONAL INPUT PARAMETERS:
;	camera=camera	scalar; type: integer; default: 2
;						camera id (1,2,3). The default camera 2 has the most
;						sky real estate.
;	from_mode=from_mode
;					scalar; type: string; default: 'equ'
;						selects the type of skymap to be used.
;						Can also be used for 'sky' to work from the
;						unsubtracted maps, but that doesn't make much sense
;						since these maps still contain all the catalogue stars.
;	to_mode=to_mode
;					scalar; type: string; default: 'bol'
;						This is only used in the naming convention for the
;						output 'pnt' files: the names will look like:
;						c2bol_YYYY_DOY_hhmmss.pnt
;						where YYYY_DOY_hhmmss is the start time of a SMEI
;						orbit, taken from the associated 'equ' file
;	destination=destination
;				 	scalar; type: string; default: $TUB
;						destination directory of output file(s)
;	fix_centroid=fix_centroid
;					if set, this turns on the 'fix_centroid' option in the
;					smei_star_fit, i.e. the centroid position of the new
;					object is iteratively improved. This gives better results
;					at the expense of a substantial increase in processing time.
;
;	cvmin=cvmin		scalar; type: float; default: 0.7
;						Only new objects with a PSF correlation coefficient
;						greater than cvmin will be written to the output files
;
;	/use_catalogue	if set, an attempt is made to match unknown objects against
;					objects in the SMEI catalogue. This will pick up 'unknown'
;					objects at the location of a SMEI catalogue object, usually
;					indicating that the catalogue star was imperfectly subtracted.
;					This keyword can increase processing time significantly
;					especially if combined with a low cvmin value.
;
;	npeak=npeak		scalar; type: integer; default: none, i.e. no limit
;					passed to href=FindPeaks=.
;					Can be used to limit the number of objects on the initial
;					list of local maxima in the skymap, that are further
;					investagated by smei_star_fit.
;					Probably best avoided (see PROCEDURE)
;
;	silent=silent	scalar; type: integer; default: 0
;						higher values suppress more informational messages
; OUTPUTS:
;	Files in destination directory
; INCLUDE:
	@compile_opt.pro		; On error, return to caller
; CALLS:
;	InitVar, smei_filename, smei_filepath, smei_getfile, hide_env
;	BadValue, smei_sky, TimeSet, sphere_distance, smei_star_fit
; PROCEDURE:
;	The equ files have all SMEI catalogue stars (presumably all stars
;	brighter than 6th magnitude), the sidereal background and the
;	zodiacal dust cloud subtracted.
;
;	The search procedure consists of the following steps:
;
;	1. Smoothing of skymaps.
;
;	First the map is smoothed by subtracting a running mean for an 8 degree wide box.
;	This takes out variations over spatial scales much larger than the PSF. This step
;	mainly serves to create local maxima that are more easily sorted, i.e. it makes
;	it easier to process maxima starting with the highest and working down from there.
;
;	The second step is to apply a smoothing using a running mean of 1 degree, the full
;	width of the SMEI PSF. This reduces the noise in the maps, and avoids picking
;	up every noise spike as a separate local maxima.
;
;	2. Prepare a list of candidate objects.
;
;	In this step the procedure FindPeaks is used to make a list of all the local
;	maxima in the skymap, listing the highest ones first.
;	The number of local maxima can be controlled with a number of keywords
;	(see href=FindPeaks= for details). This should not really be necessary:
;	step 1 has proven to be the simplest way to control the output from
;	FindPeaks.
;
;	3. Fit all the candidate maxima to the PSF using href=smei_star_fit=.
;
;	This step uses the default star fitting procedure for SMEI to identify
;	those objects on the candidate list that look like point sources (i.e.
;	look like the SMEI PSF).
;
;	The main keyword to control this step is /fix_centroid.
;	The default mode is NOT to use this keyword. In that case the PSF is fit
;	once at the location of each candidate maximum. If /fix_centroid is ON,
;	then an iterative procedure is used to find the location in the skymap
;	near the candidate maximum that provides the best fit (the highest
;	correlation between the object in the skymap and the standard SMEI PSF).
;
;	The advantage of NOT setting /fix_centroid is a significant reduction
;	in processing time (close to a factor of 10). The disadvantage is that
;	the quality of the PSF fitting is worse, which might lead to good objects
;	going unnoticed because the correlation is much worse.
;
;	The result of step 3 is that all the candidate objects from step 2
;	can be classified by sorting on the correlation coefficient of
;	of the PSF fit. The difference between using (or not using) /fix_centroid
;	is that the candidates will be sorted differently.
;
;	4. Retain the best fitting objects only
;
;	At this stage the keyword cvmin is used to only retain objects that
;	resulted in a correlation coefficient higher than cvmin (0.7 by default).
;	The total list of candidates for one skymap can easily be on the order
;	of 1000 or so. A significant fraction of the correlation coefficient
;	(often 70-80%) of the candidates will be quite low, so it makes sense
;	to reduce the search space by focusing on the higher correlations.
;
;	Note that at this stage the effect of using (or not using) fix_centroid
;	in step 3 becomes evident: since the correlation coefficients are different
;	(lower when fix_centroid is not used), a different set of candidate objects
;	is retained. Note that if cvmin=-1 is used, all objects are retained.
;	This results in significantly larger output files though.
;
;	[Rebekah, could you summarize your tests with /fix_centroid and cvmin here.
;	This is what I remember]
;
;	For cvmin=0.7 a testrun indicates that approximately the same amount of
;	objects remain for runs with /fix_centroid ON and OFF, with about 2/3 of
;	the objects the same in both lists.
;
;	4. Match against catalogue stars
;
;	This step is optional. By default it is not done, but it can be switched
;	on with keyword /use_catalogue.
;
;	The locations of the candidate maxima from step 4 (after trimming the list
;	using cvmin), is compared against the locations in the SMEI star
;	catalogues. If a catalogue star is found whitin 0.5 degrees of a
;	candidate object, most likely the object found is the results of
;	a badly-subtracted catalogue star. In this case the candidate
;	object will be listed in the output file under the name of the candidate
;	object. Note that this is not a positive indentification; it just indicates
;	the proximity of a catalogue star.
;
;	Using a low value of cvmin in combination with /use_catalogue keyword is
;	not advisable: checking a long list of candidate objects against the SMEI
;	star catalogue will drastically increase processing time.
;
;	5. Write output file
;
;	The output files contains the results of the fits for all candidate
;	objects, sorted by correlation (the best one first).
;
;	Note that for objects located in the overlapping area between the polar
;	maps and the equatorial map (approximately between 50 and 60 degrees
;	declination north and south there will be two records in the output file:
;	one for the object seen in the equatorial map, the other in one of the
;	polar maps. This needs to be taken into account when analyzing the
;	content of the input file (e.g. when counting the number of hits for
;	a given object in a sequences of skymaps).
;
;	The names of the objects indicated in which map the object was detected,
;	EQ (for the equatorial map), NP and SP for north and south polar map,
;	respectively; followed by an integer that counts the number of objects.
;	Note that the same object will have	a different name in output files
;	for different skymaps. Entries will need to be matched by RA and dec.
;
; MODIFICATION HISTORY:
;	FEB-2008, Paul Hick (UCSD/CASS)
;	NOV-2011, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;		Thorough review of code; added documentation.
;		Most significant change to the code was the removal of the zodiacal
;		dust cloud subtraction. The equ file now contain maps with the zld
;		already subtracted, so this is not needed anymore.
;       Reduced dependence on 'sky' files. These are needed now only if
;       keyword /use_catalogue is used.
;-

InitVar, camera		, 2
InitVar, from_mode	, 'equ'
InitVar, to_mode	, 'bol'
InitVar, destination, getenv('TUB')
InitVar, fix_centroid, /key
InitVar, cvmin		, 0.7
InitVar, use_catalogue, /key
InitVar, silent 	, 0

; Pick up list of SMEI skymaps.

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

; No maps found? Then stop.

IF count EQ 0 THEN	$
	RETURN

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

badval  = BadValue(0.0)
missing = badval

dpr 		= 180.0/!pi
large_smooth= 81
small_smooth= 11
mindist		= 10
relative	= 1
fraction	= 0.20

flat = ([1.0,0.02])[relative]

FOR imap=0,count-1 DO BEGIN			; Loop over all skymaps

	T0 = TimeSystem(/silent)

	equfile = given_map[imap]
	equhide = hide_env(equfile)

	equtime = smei_filename(equfile,camera=camera,postfix=postfix)

	; 'skyfile' is used only to fit catalogue stars to be matched against the unknown
	; objects found in 'equfile'

	skyfile = smei_filepath(equtime,camera=camera,postfix=postfix,mode='sky',/full)
	eclfile = smei_filepath(equtime,camera=camera,postfix=postfix,mode='ecl',/full)

	bolfile = filepath(root=destination,smei_filename(equtime,camera=camera,mode=to_mode,type='.txt'))

	IF silent LE 2 THEN message, /info, equhide

	; The psfn_sky, fovx_sky and time_sky array are low-res (720 x 360) arrays

	smei_sky, equfile, /psfn   , sky=psfn_sky, /noplot, camera=camera, /usetime, /exists, /degrees, hdr=hdr
	smei_sky, equfile, /fovx   , sky=fovx_sky, /noplot, camera=camera, /usetime, /exists, /degrees

	destroyvar, tmp
	smei_sky, equfile, /orbsecs, sky=time_sky, /noplot, camera=camera, /usetime, /exists, exten_no=tmp
	torigin = TimeSet( fxpar(headfits(equfile,exten=tmp),'TORIGIN') )

	destroyvar, new_star

	; ======================

	FOR jmap=0,2 DO BEGIN			; Equatorial, north pole, south pole

		ikeys = [0,0,0]
		ikeys[jmap] = 1

		ckeys = ['eq','np','sp']

		smei_sky, equfile, $		; Get star-subtracted skymap
			equatraw	= ikeys[0]	, $
			northraw	= ikeys[1]	, $
			southraw    = ikeys[2]	, $
			sky 		= sky		, $
			nbin		= 1 		, $
			/noplot 				, $
			hdr 		= hdr		, $
			camera		= camera	, $
			/usetime				, $
			/exact

		hdr = hdr[0]

		IF jmap EQ 0 THEN tt = TimeSet(smei=hdr.orbit)

		rr = cvsmei(from_map='MAP',from_mode=ckeys[jmap],/to_map,to_mode='lores',/silent)
		rr = round(rr)				; Nearest lowres grid point

		; Direction cosine angle for all bins in hires sky map

		fovx = fovx_sky[reform(rr[0,*]),reform(rr[1,*])]

		; This removes a small section of sky near the edge of the fov in the long dimension.
		; This area tends to be filled with residuals from bad subtractions resulting from
		; stars so close to edge that the PSFs are only partially present, or heavily deformed.

		i = where(finite(sky) AND abs(fovx) GT 29)
		IF i[0] NE -1 THEN sky[i] = BadValue(sky)

		; Skip empty maps

		i = where(finite(sky))
		IF i[0] EQ -1 THEN continue

		; 1. Smooth the sky maps.

		sky -= smooth(sky,large_smooth,/nan)		; Subtract large scale smoothed map
		sky  = smooth(sky,small_smooth,/nan)		; Smooth over 1-deg

		; Find local maxima (as array indices into skymap)

		rr = FindPeaks(sky, mindist=mindist, flat=flat, relative=relative, $
			fraction=fraction, count=nfit, npeak=npeak)

		rr = ArrayLocation(rr,dim=size(sky,/dim))
		rr = cvsmei(from_map=rr,mode=ckeys[jmap],/to_equatorial,/degrees,/silent)

		rr = AngleUnits(from_degrees=rr, /to_almanac, /singlesign)
		rr = float(rr)

		rr[2,*,*] += rr[3,*,*]/1000.0
		rr = rr[0:2,*,*]

		star_list = replicate({smei_star_list},nfit)

		star_list.name	= strupcase(ckeys[jmap])+string(format='(I4.4)', lindgen(nfit)+1)+'      '

		star_list.ra	= reform(rr[*,0,*])
		star_list.dec	= reform(rr[*,1,*])

		; Try to fit local maxima with PSF
		; with /fix_centroid set

		smei_star_fit, equfile		, $
			northpole	= ikeys[1]	, $
			southpole	= ikeys[2]	, $ 
			star_list	= star_list , $
			/degrees				, $
			psf_map 	= psfn_sky	, $
			fovx_map	= fovx_sky	, $
			time_map	= time_sky	, $
			torigin		= torigin	, $
			/use_weights			, $
			/auto_wing				, $
			fix_centroid=fix_centroid, $
			star_fit	= star_fit	, $
			star_info	= i			, $
			/silent

		; Save info for equatorial map

		IF jmap EQ 0 THEN star_info = i

		; Add centroid offset to centroid position
		; Map RA back to range [0,360]

		star_fit.radec += star_fit.dradec
		star_fit.radec[0] = AngleRange(star_fit.radec[0],/degrees)
		star_fit.dradec = 0

		; Retain only stars with correlation better than cvmin

		i = where(star_fit.cvfit GT cvmin,nfit)
		IF nfit EQ 0 THEN continue

		; If use_catalogue is NOT set then no more work needs to be done.

		IF NOT use_catalogue THEN BEGIN
			boost, new_star, star_fit
			continue
		ENDIF

		; use_catalogue IS set.
		; Continue by trying to match each of the objects against the
		; list of catalogues star. If a catalogue star is found within
		; 0.5 a degree or so, we are probably picking up a catalogue
		; object that was subtracted imperfectly.
		; Give the object the name of the catalogue star
		; This step is very time consuming if cvmin is set very low,
		; which is why the use_catalogue keyword was introduced.

		star_fit = star_fit[i]

		; Check for stars already on the star catalogue. Quite a few
		; of the stars just detected are probably stars that were not
		; fitted very well (due to bad quaternions for instance).

		star_cat = smei_star_info(/get_struct,count=ncat)

		cat_pos = smei_star_info(star_cat,/radec,/degrees)
		fit_pos = star_fit.radec

		IF nfit GT 1 THEN BEGIN
			cat_pos = SuperArray(cat_pos,nfit,after=2)
			fit_pos = SuperArray(fit_pos,ncat,after=1)
		ENDIF

		; Elongations of new stars vs catalogue stars

		elo = sphere_distance(cat_pos,fit_pos,/degrees) ; array[ncat,nfit]

		destroyvar, cat_pos, fit_pos

		; Find catalogue star closes to new stars

		tmp = min(elo,p,dim=1)		; p = array[nfit]

		p = ArrayLocation(p,dim=[ncat,nfit])
		p = reform(p[0,*])

		; Closest catalogue stars

		star_cat = star_cat[p]

		; Fit closest stars in original (unsubtracted) skymap
		; with /fix_centroid set to get the best possible location


		;names = star_cat.name


		smei_star_fit, skyfile		, $
			northpole	= ikeys[1]	, $
			southpole	= ikeys[2]	, $ 
			star_list	= star_cat	, $
			/degrees				, $
			psf_map 	= psfn_sky	, $
			fovx_map	= fovx_sky	, $
			time_map	= time_sky	, $
			torigin		= torigin	, $
			/use_weights			, $
			/auto_wing				, $
			fix_centroid=fix_centroid, $
			count		= count_cat	, $
			star_fit	= close_cat , $
			max_stars	= 1			, $
			star_remove	= replicate(0,nfit), $
			star_index	= star_index, $
			skip_dist	= 0.0		, $
			incl_dist	= 0.0		, $
			/silent

		IF count_cat EQ 0 THEN continue

		star_cat = star_cat[star_index]
		star_fit = star_fit[star_index]

		; Add centroid offset to centroid position

		close_cat.radec += close_cat.dradec
		close_cat.radec[0] = AngleRange(close_cat.radec[0],/degrees)
		close_cat.dradec = 0

		elo = sphere_distance(star_fit.radec,close_cat.radec,/degrees)

		p = where(elo LT 0.5,n)

		FOR i=0,n-1 DO BEGIN
			q = p[i]
			print, star_fit[q].name,'=',close_cat[q].name,' @ ',elo[q],' deg'
			star_fit[q].name = close_cat[q].name
		ENDFOR

		boost, new_star, star_fit

	ENDFOR

	; Write the output file with candidate objects
	; Objects are written sorted by correlation coefficient (high first).

	i = n_elements(new_star)
	IF i NE 0 THEN BEGIN
		message, /info, 'writing '+hide_env(bolfile)+' ('+strcompress(i,/rem)+' stars)'
		smei_star_writepnt, bolfile, new_star[reverse(sort(new_star.cvfit))], star_info=star_info
	ENDIF

ENDFOR

RETURN  &  END