;+
; NAME:
;	smei_sky_cleanedge_fix
; PURPOSE:
;	Applies a correction to a hires skymap (equatorial
;	north pole or south pole) for straylight contribution
;	from sky just outside the fov.
; CATEGORY:
;	camera/idl/toolbox
; CALLING SEQUENCE:
	PRO smei_sky_cleanedge_fix		, $
		lo_fovx 	, lo_time		, $ 	; Lores direction cosine angle and time maps
		hi_sky		, 				  $ 	; Hires sky map
		lo_edge 	, hi_edge		, $ 	; Inner and outer edge of band in need of fixing
		fraction	, time_sample	, $ 	; Max fraction of sample time series subtracted
		left_sample , right_sample	, $
		mode	= mode
; INPUTS:
;	lo_fovx 			array[720,360]; type: float
;	lo_time  			array[720,360]; type: float
;							2D lores skymaps: time and fovx angle
;	hi_sky				array[3600,1200,n] or array[800,800]
;							2D hires skymap to be corrected
;							If n=1 the 3rd dimension may be absent
;
;	lo_edge 			scalar; inner edge of fov area that needs
;							to be corrected (lo_edge=20 degrees)
;	hi_edge 			scalar: outer edge of fov area that needs
;							to be corrected (hi_edge=30 degrees)
;	time_sample 		array[n]; type: float
;							times (as seconds since start of orbit)
;							covering the entire orbit
;	left_sample 		array[n]; type: float
;							sky brightness at "left edge" of fov
;							(negative direction cosine angles)
;							set by href=smei_sky_cleanedge_median=)
;	right_sample		array[n]; type: float
;							sky brightness at "right edge" of fov
;							(positive direction cosine angles)
;							set by href=smei_sky_cleanedge_median=)
; OUTPUTS:
;	hi_sky				array[3600,1200] or array[800,800]
;							corrected 2D hires skymap
; INCLUDE:
	@compile_opt.pro		; On error, return to caller
; CALLS:
;	IsType
; RESTRICTIONS:
;	The units of the angles (degrees or radians) must be the
;	same for lo_fovx, lo_edge and hi_edge
; PROCEDURE:
;	All bins in the hires sky map with direction cosine angle
;	greater than lo_edge in absolute value are corrected by
;	subtracted a fraction of the left_sample or right_sample
;	array.
;
;	The correction is largest at the edge of the fov
;	and drops of linearly with direction cosine angle to zero
;	at direction cosine angle lo_edge.
;
;	At direction cosine angle hi_edge the correction is
;	'fraction' times left_sample or right_sample.
;
;	For each skybin the time is used to determine the
;	appropriate left_sample or right_sample value by
;	linear interpolation on the two sample time series.
; MODIFICATION HISTORY:
;	APR-2007, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;-

ndim = size(hi_sky,/n_dim)			; Number of dimensions
sdim = size(hi_sky,/dim)

pos = cvsmei(from_map='MAP',from_mode=mode,/to_map,to_mode='lores',/silent)
pos = round(pos) 					; Nearest lowres grid point

xpos = reform(pos[0,*]) 			; Hor  index of nearest lowres grid point
ypos = reform(pos[1,*]) 			; Vert index of nearest lowres grid point

fovx = lo_fovx[xpos,ypos]			; Direction cosine angles for all hires sky bins

fpd = fraction/(hi_edge-lo_edge)	; Fractional correction per degree
									; Left edge (negative direction cosine angles)
IF IsType(left_sample,/defined) THEN BEGIN

	edge = where(fovx LT -lo_edge)	; Hires sky bins along left edge

	IF edge[0] NE -1 THEN BEGIN
									; Find sample sky brightness by linear interpolation
		edge_sky = fpd*interpol(left_sample,time_sample,lo_time[xpos[edge],ypos[edge]])

		; Set missing data to zero: no correction is made where
		; sample brightness is not available

		tmp = where(1-finite(edge_sky))
		IF tmp[0] NE -1 THEN edge_sky[tmp] = 0.0

		tmp = (-fovx[edge]-lo_edge)*edge_sky

		CASE ndim OF
		2: hi_sky[edge] -= tmp
		3: BEGIN
			hi_sky = reform(hi_sky,sdim[0]*sdim[1],sdim[2],/overwrite)
			FOR n=0,sdim[2]-1 DO hi_sky[edge,n] -= tmp
			hi_sky = reform(hi_sky,sdim,/overwrite)
		END
		ENDCASE

	ENDIF

ENDIF
									; Right edge (positive direction cosine angles)
IF IsType(right_sample,/defined) THEN BEGIN

	edge = where(fovx GT  lo_edge)	; Hires sky bins along right edge

	IF edge[0] NE -1 THEN BEGIN
									; Find sample sky brightness by linear interpolation
		edge_sky  = fpd*interpol(right_sample,time_sample,lo_time[xpos[edge],ypos[edge]])

		; Set missing data to zero: no correction is made where
		; sample brightness is not available

		tmp = where(1-finite(edge_sky))
		IF tmp[0] NE -1 THEN edge_sky[tmp] = 0.0

		tmp = (fovx[edge]-lo_edge)*edge_sky

		CASE ndim OF
		2: hi_sky[edge] -= tmp
		3: BEGIN
			hi_sky = reform(hi_sky,sdim[0]*sdim[1],sdim[2],/overwrite)
			FOR n=0,sdim[2]-1 DO hi_sky[edge,n] -= tmp
			hi_sky = reform(hi_sky,sdim,/overwrite)
		END
		ENDCASE

	ENDIF

ENDIF

RETURN  &  END