;+
; NAME:
;	sphere_smooth
; PURPOSE:
;	Smooths data on a sphere
; CATEGORY:
;	sat/idl/toolbox/math
; CALLING SEQUENCE:
	FUNCTION smooth_sphere, FF	, $

		irange		= irange 	, $
		igrid		= igrid		, $ 	; Input grid

		odim		= odim_		, $ 	; Output grid
		orange		= orange	, $
		ogrid		= ogrid 	, $

		width		= width		, $
		gauss		= gauss		, $
		degrees		= degrees	, $

		weightfnc	= weightfnc , $

		fillall 	= fillall	, $
		fillbad 	= fillbad	, $
		fillgood	= fillgood	, $

		min_weight	= min_weight
; INPUTS:
;	FF				array[N,M,L]; type: any numerical type
;						N: longitudinal dimension
;						M: latitudinal dimension
;						L: number of arrays to be smoothed
; OPTIONAL INPUT PARAMETERS:
;	irange=irange	array[2,2]; type: float; default: [[0,360],[-90,90]]
;					Specifies the range of the input grid on which
;					F is defined
;						irange[*,0]: longitude range
;						irange[*,1]: latitude range
;					The grid specifies a skymap in
;					equirectangular projection
;	igrid=igrid		array[2,N,M]; type: float
;					Instead of giving irange, the grid can be specified
;					explicitly:
;						igrid[0,*,*]: longitudes
;						igrid[1,*,*]; latitudes
;					Note that the grid must still be equirectangular
;
;	odim=odim		array[2]; type: integer; default: [N,M]
;						dimensions of smoothed output skymap
;						The dimensions must be smaller then/equal to
;						the input dimensions and must be a factor of
;						the input dimensions.
;	orange=orange	array[2,2]; type: float; default: irange
;					Specifies the range of the output grid on which
;					R is defined.
;						orange[*,0]; longitude range
;						orange[*,1]: latitude range
;					Again, this is a equirectangular grid.
;	ogrid=ogrid 	array[2,n,m]; type: float
;					Instead of giving odim and orange, the output
;					grid can be specified explicitly:
;						grid[0,*,*]: longitudes
;						grid[1,*,*]; latitudes
;	/degrees		if set, all angles are in degrees; default: radians
;
;	min_weight		threshold for the weights. All contributing skybins
;					with weight less then min_weight will be ignored.
; OUTPUTS:
;	R				array[n,m]; type: same as FF
;						smoothed output skymap
; OPTIONAL OUTPUT PARAMETERS:
;	igrid=igrid		if not defined on input, they are returned as calculated	
;	ogrid=ogrid		.. from the *dim and *range keywords.
; INCLUDE:
	@compile_opt.pro		; On error, return to caller
; CALLS:
;	InitVar, ToRadians, ToDegrees, sphere_distance, AngleRange
;	ArrayLocation, gridgen
; PROCEDURE:
; >	If neither /gauss, nor weight_fnc is set than all points within
;	'width' (angular separation) of a grid point will be averaged with
;	the same weight.
; >	If /gauss is set a Gaussina weight is used in the calculation
;	of the average (with weight exp(-(elo/width)^2). Contributing points
;	will extend to about 3*width away in angular separation.
; >	A user-defined function weight_fnc can be used to calculate the weight.
;	The function should have the form
;		W = weight_fnc(elo, width=width, degrees=degrees)
;	with 'elo' (in units indicated by keyword 'degrees') the elongation,
;	and 'width' the input argument 'wdith' to this function.
; MODIFICATION HISTORY:
;	JAN-2008, Paul Hick (UCSD/CASS)
;	APR-2008, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;		Bugfix. Weighting was done wrong.
;-

InitVar, fillbad , /key
InitVar, fillgood, /key
InitVar, fillall , /key

fillall OR= 1-fillbad AND 1-fillgood

InitVar, degrees, /key
rpm = float(ToRadians(degrees=degrees))
dpm = float(ToDegrees(degrees=degrees))

InitVar, min_weight, 0.0

; Pick up input dimensions.
; Output dimensions are same as input, if not specified

idim = size(FF,/dim)
InitVar, odim_, idim, set=odim

more_than_one = n_elements(idim) EQ 3

CASE more_than_one OF
0: n_layer = 1
1: BEGIN
	n_layer = idim[2]
	IF n_elements(odim) EQ 2 THEN odim = [odim,n_layer]
END
ENDCASE

; If no function to calculate weights was defined explicitly,
; then do either a "box" smoothing (weight 1 inside, weight 0
; outside the box), or us a Gaussian smoothing

have_fnc = IsType(weightfnc,/defined)
IF NOT have_fnc THEN BEGIn
	InitVar, width, 1.0/dpm
	InitVar, gauss, /key
ENDIF

; The input grid

CASE IsType(igrid,/defined) OF

0: BEGIN
	InitVar, irange, [ [0.0,360.0], [-90.0,90.0] ]/dpm
	igrid = gridgen(idim[0:1], range=irange, /open, /multid)
END

1: irange = [ [igrid[0,[0,idim[0]-1],0]],[igrid[1,0,[0,idim[1]-1]]] ]

ENDCASE

lng_igrid = reform(igrid[0,*,*])
lat_igrid = reform(igrid[1,*,*])

ilng = reform(igrid[0,*,0])
ilat = reform(igrid[1,0,*])

; Grid resolutions (deg/bin) for input grid

dilng = (igrid[0,idim[0]-1,0]-igrid[0,0,0])/(idim[0]-1)
dilat = (igrid[1,0,idim[1]-1]-igrid[1,0,0])/(idim[1]-1)

; Index for latitude 0 in input array (not an integer in general)

ilat_zero = -0.5-irange[0,1]/(irange[1,1]-irange[0,1])*idim[1]

; Number of skybins covering 360 degrees in
; longitude in input grid

ilng_360 = round(360.0/dpm/dilng)

ilng_full360 = idim[0] EQ ilng_360
IF ilng_full360 THEN message, /info, 'input map covers 360 deg in longitude'

CASE IsType(ogrid,/defined) OF

0: BEGIN
	InitVar, orange, irange
	ogrid = gridgen(odim[0:1], range=orange, /open, /multid)
END

1: orange = [ [ogrid[0,[0,odim[0]-1],0]],[ogrid[1,0,[0,odim[1]-1]]] ]

ENDCASE

olng = reform(ogrid[0,*,0]) 		; Output longitudes
olat = reform(ogrid[1,0,*]) 		; Output latitudes

; Grid resolutions (deg/bin) for output grid

dolng = (ogrid[0,odim[0]-1,0]-ogrid[0,0,0])/(odim[0]-1)
dolat = (ogrid[1,0,odim[1]-1]-ogrid[1,0,0])/(odim[1]-1)

; Only allow rebinning to lower resolution
; (smaller array size)

IF dolng LT dilng OR dolat LT dilat THEN BEGIN
	message, /info, 'resolution of results must be coarser than input resolution'
	RETURN, -1
ENDIF

tmp = [dolng,dolat]/[dilng,dilat]
nbin = round(tmp)

; Only allow rebinning to dimensions that are a factor
; of the input dimension.

IF abs(tmp[0]-nbin[0]) GT 0.0001 OR abs(tmp[1]-nbin[1]) GT 0.001 THEN BEGIN
	message, /info, 'dimension of result must be factor of the input dimension'
	RETURN, -1
ENDIF

; For Gaussian weights or if an explicit weight function
; is to be used, make sure to pick a big enough box.

box_width = width*(1+2*max([gauss,have_fnc]))

; Set a box size on the equirectangular map that include
; all the bins that are potentially to be included.
; In the longitudinal direction the width scales as cos(lat)^-1
; with a maximum of 180 degrees.

lat_max = acos(box_width*dpm/180.0)/rpm
lng_width = box_width/cos(((abs(olat)+0.5*box_width) < lat_max)*rpm)
lat_width = box_width

; Find range in output grid that overlaps with the
; input grid and determine which part of the output
; skymap needs to be filled.

; Nearest points in input grid for points in output grid

ilng_near_olng = round((olng-ilng[0])/dilng)
ilat_near_olat = round((olat-ilat[0])/dilat)

CASE ilng_full360 OF

0: BEGIN

	tmp = where( 0 LE ilng_near_olng AND ilng_near_olng LT idim[0])
	IF tmp[0] EQ -1 THEN BEGIN
		message, /info, 'output longitudes do not overlap input longitudes'
		RETURN, -1
	ENDIF

END

1: tmp = indgen(odim[0])

ENDCASE

min_olng_index = min(tmp,max=max_olng_index) 
ref_olng_index = (min_olng_index+max_olng_index)/2
n_olng = max_olng_index-min_olng_index+1

unit = replicate(1,n_olng)
step = gridgen(n_olng, range=[min_olng_index,max_olng_index]-ref_olng_index)

tmp = where( 0 LE ilat_near_olat AND ilat_near_olat LT idim[1])
IF tmp[0] EQ -1 THEN BEGIN
	message, /info, 'output latitudes do not overlap input latitudes'
	RETURN, -1
ENDIF

min_olat_index = min(tmp,max=max_olat_index) 
n_olat = max_olat_index-min_olat_index+1

bad = BadValue(FF)

R = make_array(dim=odim, value=bad)

igrid = reform(igrid,2,idim[0]*idim[1],/overwrite)

; In principle, the do-loop can be reduced further by using North/South symmetry.

olat_done = bytarr(n_olat)
epsilon = 1.0e-6

FOR j=min_olat_index,max_olat_index DO BEGIN 	; Loop over output latitudes

	; Check whether latitude has been done already

	IF olat_done[j] THEN continue

	; Check whether opposite latitude is present.
	; If it is this gets done here too.
	; (This potentially reduces the loop over output latitudes by
	; a factor of two. There is some speed gain, but not much,
	; and this doubles the memory requirements, so I'm not so sure
	; this is a good idea; disable for now).
 
	j_other = -1
	;CASE abs(olat[j]) GT epsilon OF
	;0: j_other = -1
	;1: j_other = (where(abs(olat+olat[j]) LT epsilon))[0]
	;ENDCASE

	olat_done[j] = 1
	IF j_other NE -1 THEN olat_done[j_other] = 1

	; Check the reference longitude at each output latitude:
	; Reference point is P = [olng[ref_olng_index],olat[j]]
	; Find all points in the input map inside the box centered on P.

	ipos = AngleRange(lng_igrid-olng[ref_olng_index],/pi,degrees=degrees)

	ipos = where(abs(ipos) LE lng_width[j] AND abs(lat_igrid-olat[j]) LE lat_width, n_ipos)
												; Array [n_ipos]
	IF n_ipos EQ 0 THEN continue 				; Shouldn't happen

	; Elongations for all n points in box.

	W = sphere_distance( igrid[*,ipos], [olng[ref_olng_index],olat[j]], degrees=degrees)
	W = float(W)

	; Weights for all n points in box

	CASE 1 OF

	have_fnc:		$
		W = call_function(weight_fnc, W, width=width, degrees=degrees)

	gauss: BEGIN								; Gaussian weights
		W /= width
		W = exp(-W*W)
	END

	ELSE:			$ 							; "Box" weighting
		W LE= width

	ENDCASE

	; This is primarily intended to reduce memory requirements
	; for subsequent calculations: ditch points with very low weights

	tmp = where(W GT min_weight)
	IF tmp[0] NE -1 THEN BEGIN
		n_ipos = n_elements(tmp)
		ipos = ipos[tmp]
		W    = W   [tmp]
	ENDIF

	; Array indices into input array FF for n points used
	; for P = [olng[ref_olng_index],olat[j]]

	ipos = ArrayLocation(ipos,dim=idim[0:1]) 	; Array [2,n_ipos]
	opos = [ref_olng_index,j]			; Array [2]

	; If the opposite latitude is to be included, set up an 2nd group.
	; The olng indices remain the same; the olat index differences
	; with olat[j_other] have the opposite sign as for olat[j].

	n_opos = 1+(j_other NE -1)

	IF n_opos NE 1 THEN BEGIN
		ipos = SuperArray(ipos,n_opos,/trail)	; Array [2,n_ipos,n_opos]
		ipos[1,*,1] = round(2*ilat_zero-ipos[1,*,1]); Adjust lat index for opposite latitude
		W = SuperArray(W,n_opos,/trail)			; Array [n_ipos,n_opos]
		opos = SuperArray(opos,n_opos,/trail)	; Add ouput indices for opposite latitude
		opos[1,1] = j_other 					; Array [2,n_opos]
	ENDIF

	; Replicate the weights as n_olng identical rows.
	; Each row represents one of the longitudes between
	; min_olng_index and max_olng_index.

	W = SuperArray(W,n_olng,/trail) 			; Array [n_ipos,n_opos,n_olng]

	; There are n_olng longitudes in the output array at latitute olat[j].
	; We just calculated the points to be included in calculating the
	; smoothed value at lng[ref_olng_index].
	; The other longitudes will include the same grouping of points, but
	; shifted in longitude by a multiple of nbin[0] input skybins.

	; Set up arrays for n_olng groups of n_ipos skybins with longitudinal
	; and latitudinal array indices into the input array FF.

	; Longitudes: n_olng rows, increasing by nbin[0] from row to row.
	; Latitudes: n_olng identical rows.

	ipos = SuperArray(ipos,n_olng,/trail)		; Array [2,n_ipos,n_opos,n_olng]
	ipos = SubArray  (ipos,element=0,add=replicate(1,n_ipos*n_opos)#(nbin[0]*step))

	opos = SuperArray(opos,n_olng,/trail)		; Array [2,n_opos,n_olng]
	opos = SubArray  (opos,element=0,add=replicate(1,n_opos)#step)

	; The longitudinal indices may be outside the valid
	; range [0,idim[0]-1].
	; If irange[1]-irange[0] = 360.0 this can be fixed by
	; adding/subtracting the equivalent of 360 to the
	; array index.

	IF ilng_full360 THEN BEGIN
		tmp = SubArray(ipos,element=0)
		tmp mod= ilng_360
		tmp   += ilng_360
		tmp mod= ilng_360
		ipos = SubArray(ipos,element=0,replace=tmp)
	ENDIF

	x = SubArray(ipos,element=0)
	y = SubArray(ipos,element=1)
	tmp = where(0 LE x AND x LT idim[0] AND 0 LE y AND y LT idim[1])

	IF tmp[0] EQ -1 THEN continue				; Shouldn't happen

	; Convert to linear index into input array

	ipos = reform(ipos,[2,n_ipos*n_opos*n_olng],/overwrite)
	ipos = ArrayLocation(ipos[*,tmp], dim=idim, /onedim)

	; Pick up function values from input array

	F = make_array(dim=[n_ipos,n_opos,n_olng,n_layer],value=bad)

	n_off = n_ipos*n_opos*n_olng
	m_off = idim[0]*idim[1]
	FOR i=0,n_layer-1 DO F[tmp+n_off*i] = FF[ipos+m_off*i]

	W = SuperArray(W,n_layer,/trail)

	; Bad values are set to zero, and are given a zero
	; weight. This effectively excludes them from
	; contributing to the weighted mean.

	tmp = where(1-finite(F))

	IF tmp[0] NE -1 THEN BEGIN
		F[tmp] = 0
		W[tmp] = 0
	ENDIF

	; Calculate weighted mean for all points at olat[j]
	; F,W = Array[n_ipos,n_opos,n_olng,n_layer]

	; Sum over n_ipos dimension

	F = total(W*F,1)/total(W,1)	 				; Array [n_opos, n_olng, n_layer]

	x = SubArray(opos,element=0) 				; Array [n_opos, n_olng]
	y = SubArray(opos,element=1) 				; Array [n_opos, n_olng]
	z = lonarr(n_opos,n_olng)					; Array [n_opos, n_olng]

	x_near = ilng_near_olng[x]
	y_near = ilat_near_olat[y]

	F = reform(F,n_opos*n_olng,n_layer,/overwrite)

	FOR i=0,n_layer-1 DO BEGIN

		CASE 1 OF

		fillall : R[x,y,i+z] = F[*,i]

		fillgood: BEGIN
			tmp = where(  finite(FF[x_near,y_near,i]) )
			IF tmp[0] NE -1 THEN R[x[tmp],y[tmp],i+z[tmp]] = F[tmp,i]
		END

		fillbad : BEGIN
			tmp = where(1-finite(FF[x_near,y_near,i]) )
			IF tmp[0] NE -1 THEN R[x[tmp],y[tmp],i+z[tmp]] = F[tmp,i]
		END

		ENDCASE

	ENDFOR

ENDFOR

RETURN, R  &  END