;+
; NAME:
;	ColorSkybox
; PURPOSE:
;	Draws a color display for a set of 'skyboxes' as defined by a grid
;	of spherical coordinates. Draws either a rectangular map
;	(x=phase,y=latitude), a fish-eye map or a Hammer-Aitoff map
; CALLING SEQUENCE:
	PRO ColorSkybox, Phase, Latitude, Color,	$
		skyedge 	= skyedge	, $
		zero_phase	= zero_phase, $
		dabg		= dabg		, $
		degrees		= degrees	, $
		black		= black		, $
		void		= void		, $
		fill2edge	= fill2edge	, $
		use_mask	= use_mask	, $
		mask		= mask		, $
		_extra		= _extra
; INPUTS:
;	Phase		array[n [+1] ] or array[n [+1],m [+1] ]; type: float
;					Phase angle; 'longitude'
;					Dimension n and/or m is used to specify the center of
;					sky box n,m (the box edges will be calculated internally
;					Dimension n+1 and/or m+1 is used to specify the edges of
;					sky box n,m
;	Latitude	array[m [+1] ] or array[n [+1],m [+1] ]; type: float
;					Latitude angle
;					Box center or box edge can be specified, as for 'Phase'
;	Color		array[n,m]; type: integer
;					2D array with color indices
;					boxes with negative color indices are never colored in
;					boxes with index zero are colored only if keyword /black is set
;	skyedge		scalar; type: float
;					if set and positive, data are plotted as a 'fish-eye' map.
;					The value is used as cut-off for the range of colatitude (polar)
;					angles plotted (polar angles above 'skyedge' not plotted);
;					if set and negative, a Hammer-Aitoff projection is used;
;					if not set or zero, a rectangular map is drawn
;	mask		array[n,m]; type: float
;					array with same dimensions as Color.
;					Usually this is a "mask" identifying "good" and "bad" areas
;					in the Color map.
;	use_mask	scalar or array; type: float
;					specifies contour levels in the "mask" array
;					These contour leves are overplotted on the color map
;					using the IDL contour function.
;	_extra=_extra	can be used to pass keywords to IDL contour function
;					processing the mask and use_mask keywords
; OPTIONAL INPUTS:
;	zero_phase=zero_phase
;				scalar, or array with same structure as 'Phase'; type: float
;					The input arrays Phase, Latitude, Color are rearranged to put
;					zero_phase in the center of the map.	
;	dabg=dabg	array[3]; type: float; default: [0,0,0]
;					Passed to href=FishEye= and href=HammerAitoff=
;					Determines the direction of the center of the plot
;	/degrees	if set, all angles are in degrees (default: radians)
;	/black		if set, then color index 0 (black) is drawn with a call to polyfill
;					By default color index 0 is not drawn (i.e. remains in the
;					background color).
;	/void		if set, then negative color index is drawn using the foreground
;					color !p.color. By default these are not drawn (i.e. remain
;					in the background color).
;	/fill2edge	if set, then the phase angle boundaries of the outermost bins are set to
;					+/- 180 degrees. Only used if bin centers are specified for
;					the phase angles.
; OUTPUTS:
;	skyedge		scalar
;					updated only if positive and bigger than 170 deg
;					(returned value is 170 deg)
; INCLUDE:
	@compile_opt.pro					; On error, return to caller
; CALLS:
;	InitVar, IsType, ToRadians, SuperArray, SubArray, AngleRange, FishEye, HammerAitoff
; RESTRICTIONS:
; >	Color must be a 2D array with dimension [N,M]. The first dimension
;	(N) is the phase angle dimension; the second (M) is the latitude
;	angle dimension
; >	The Phase and Latitude arrays can be 1-dim or 2-dim.
;	A 1-dim Phase array (n) is interpreted as a 2-dim array (n,m) with
;	all m rows the same. A 1-dim Latitude array (m) is interpreted as
;	a 2-dim array (n,m) with all n columns the same.
; >	The second dimension of Phase should be the same (n=N) or one larger
;	(n=N+1) than for Color. If n=N+1 then Phase contains the phase angles
;	of the edges of sky boxes; if n=N it contains the center phase angles,
;	and the edges are calculated internally.
; >	The second dimension of Latitude should be the same (m=M) or one
;	larger (m=M+1) than for Color. If m=M+1 then Latitude contains the
;	latitude angles of the edges of sky boxes; if m=M it contains
;	the center latitude angles, and the edges are calculated internally.
; >	The value of 'skyedge' is changed to skyedge = (skyedge < 170).
;	The method of plotting the skyboxes (by connecting the corners by
;	straight lines in the x-y plane of the plot) does not work for large
;	elongations in the fish-eye maps. Problems are avoiding by not
;	permitting 'skyedge' to be larger than 170.
; PROCEDURE:
; >	The color array indices contains color indices for NxM 'skyboxes'.
;	The corners of the boxes are stored in the Phase and Latitude arrays
; >	Phase and Latitude are usually equatorial (RA, dec) or ecliptic
;	coordinates
; >	The center of fish-eye and Hammer-Aitoff has phase angle zero_lng
;	and latitude 0.
; >	The arrays Phase and Latitude will usually be obtained by a call
;	to href=EulerRotate=
; >	The array color can be obtained by a call to href=GetColors=
; >	If the 'skyedge' keyword is set, then the angles Phase and Latitude
;	are plotted on the screen as a 'fish-eye' view.
; >	The boxes in the x-y plane of the plot are defined by connecting the
;	corners by straight lines (and using the appropriate color).
; >	Boxes with color index 0 (black) are skipped by default. If the keyword
;	'black' is set color index 0 is explicitly colored with polyfill
; MODIFICATION HISTORY:
;	AUG-1999, Paul Hick (UCSD/CASS)
;		added check for negative color indices, these are now ignored
;		(GetColors now checks for bad values using the 'finite'
;		function and sets corresponding boxes to a negative color
;		index).
;	JAN-2002, Paul Hick (UCSD/CASS)
;		Improved tests to decide which sky boxes to plot in a fish-eye
;		plot at large polar angles.
;	APR-2002, Paul Hick (UCSD/CASS)
;		Added /fill2edge keyword.
;	MAR-2011, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;		Dropped restriction on zero_phase being applied only if
;		Phase is 1-dim array.
;-

N = size(Color)
IF N[0] NE 2 THEN message, 'Color must be 2D array'

InitVar, fill2edge	, /key
InitVar, void		, /key
InitVar, black		, /key
InitVar, skyedge	, 0

has_mask = IsType(use_mask,/defined) AND IsType(mask,/defined)

rpm	= ToRadians(degrees=degrees)
dpm	= ToDegrees(degrees=degrees)
pi	= !pi/rpm
pi2	= pi/2

Xp = Phase
Yp = Latitude

szXp = size(Xp)
szYp = size(Yp)

; If zero_phase is set, then rearrange the Phase array to put zero_phase in the center.
; (but only if box centers are specified for phase angle. Why??)

IF szXp[1] EQ N[1] AND IsType(zero_phase, /defined) THEN BEGIN
	Xp -= zero_phase
	Xp = AngleRange(Xp, /pi, degrees=degrees)

	IF szXp[0] EQ 1 THEN BEGIN			; 1-dim Phase array (same phases for all latitudes)
		dp = sort(Xp)
		Xp = Xp[dp]
		Color = Color[dp,*]				; Sort all latitudes
		IF has_mask THEN mask = mask[dp,*]
		IF szYp[0] EQ 2 THEN IF szYp[1] EQ N[1] THEN Yp = Yp[dp,*]
	ENDIF ELSE BEGIN
		FOR i=0,szXp[2]-1 DO BEGIN		; Adjust all latitudes separately
			dp = sort(Xp[*,i])			; Here we assume that the input Xp is monotonic and covers at most 360 degrees.
			Xp   [*,i] = Xp   [dp,i]
			Color[*,i] = Color[dp,i]
			IF has_mask THEN mask[*,i] = mask[dp,i]
			IF szYp[0] EQ 2 THEN IF szYp[1] EQ N[1] THEN Yp[*,i] = Yp[dp,i]
		ENDFOR
	ENDELSE
ENDIF

;IF szXp[0] EQ 1 AND IsType(zero_phase, /defined) THEN BEGIN
;	IF szXp[1] EQ N[1] THEN BEGIN		; Box centers specified for phase angle
;		Xp -= zero_phase
;		Xp  = AngleRange(Xp, /pi, degrees=degrees)
;
;		; Here we assume that the input Xp is monotonic and covers at most 360 degrees.
;
;		dp = sort(Xp)
;		Xp = Xp[dp]
;		Color = Color[dp,*]
;		IF has_mask THEN mask = mask[dp,*]
;		IF szYp[0] EQ 2 THEN IF szYp[1] EQ N[1] THEN Yp = Yp[dp,*]
;	ENDIF
;ENDIF

; If necessary, make Xp and Yp 2D grids, combining every
; Xp with every Yp

IF szXp[0] EQ 1 THEN Xp = SuperArray(Xp, szYp[szYp[0]], /trail)
IF szYp[0] EQ 1 THEN Yp = SuperArray(Yp, szXp[1]      , /lead )

szXp = size(Xp)
szYp = size(Yp)

IF has_mask THEN BEGIN
	XpC = Xp
	YpC = Yp
	M = szXp[1]						; # phase angles
	IF M EQ N[1]+1 THEN BEGIN		; Box edges specified for phase angles: calculate centers
		XpC = (XpC+shift(XpC,1,0))[1:*,*]/2
		YpC = YpC[1:*,*]
	ENDIF
	M = szYp[2]						; # latitudes
	IF M EQ N[2]+1 THEN BEGIN
		XpC = XpC[*,1:*]			; Box edges specified for latitudes: calculate centers
		YpC = (YpC+shift(YpC,0,1))[*,1:*]/2
	ENDIF
ENDIF

M = szXp[1]							; # phase angles

IF M EQ N[1] THEN BEGIN				; Box centers specified for phase angles: calculate edges

	dp = (Xp[M-1,*]-Xp[0,*])/(M-1)	; Average phase angle separation of centers at at all latitudes

	; Internal edges are set halfway between the centers

	Xp = (Xp+shift(Xp,1,0))[1:*,*]/2

	; Use the average separation to set the outer edges

	XpMin = ( (Xp[  0,*]-dp) > (-pi) )*(1-fill2edge)-pi*fill2edge
	XpMax = ( (Xp[M-2,*]+dp) < (+pi) )*(1-fill2edge)+pi*fill2edge

	;Xp = [ (Xp[0,*]-dp) > (-pi), Xp, (Xp[M-2,*]+dp) < pi
	Xp = [ XpMin, Xp, XpMax ]
	dp = (Yp[M-1,*]-Yp[0,*])/(M-1)
	Yp = (Yp+shift(Yp,1,0))[1:*,*]/2
	Yp = [ (Yp[0,*]-dp) > (-pi2), Yp, (Yp[M-2,*]+dp) < pi2 ]

	szXp = size(Xp)
	szYp = size(Yp)
ENDIF


M = szYp[2]							; # latitudes
IF M EQ N[2] THEN BEGIN				; Box centers specified for latitude angles: calculate edges
	Xp = transpose(Xp)
	Yp = transpose(Yp)

	dp = (Xp[M-1,*]-Xp[0,*])/(M-1)
	Xp = (Xp+shift(Xp,1,0))[1:*,*]/2
	Xp = [ (Xp[0,*]-dp) > (-pi), Xp, (Xp[M-2,*]+dp) < pi ]

	dp = (Yp[M-1,*]-Yp[0,*])/(M-1)
	Yp = (Yp+shift(Yp,1,0))[1:*,*]/2
	Yp = [ (Yp[0,*]-dp) > (-pi2), Yp, (Yp[M-2,*]+dp) < pi2 ]

	Xp = transpose(Xp)
	Yp = transpose(Yp)

	szXp = size(Xp)
	szYp = size(Yp)
ENDIF

; Xp is a 2D array with phase angles and Yp a 2D array with latitude angles
; They refer to the edges of the sky boxes specified in Color
; XpC and YpC refer to the centers of the boxes

BinOK = Color GT 0 OR (black AND Color EQ 0) OR (void AND Color LT 0)

Xp = SuperArray(Xp,2,/lead)
Xp[1,*,*] = Yp

IF has_mask THEN BEGIN
	XpC = SuperArray(XpC,2,/lead)
	XpC[1,*,*] = YpC
ENDIF

IF skyedge GT 0 THEN BEGIN			; Fish eye plot

	skyedge = (skyedge < 170/dpm) 	; Exclude large latitude angles
									; Scales skyedge to 2*sqrt(2)
	Xp  = FishEye(Xp , dabg=dabg, maxelo=skyedge, degrees=degrees, polar=Yp)
	IF has_mask THEN XpC = FishEye(XpC, dabg=dabg, maxelo=skyedge, degrees=degrees)

	; Polar coordinates of 4 corners

	c00 = (      Yp         )[*,0:N[1]-1,0:N[2]-1]
	c10 = (shift(Yp,0,-1, 0))[*,0:N[1]-1,0:N[2]-1]
	c01 = (shift(Yp,0, 0,-1))[*,0:N[1]-1,0:N[2]-1]
	c11 = (shift(Yp,0,-1,-1))[*,0:N[1]-1,0:N[2]-1]

	; Phase differences with first coordinates

	c10[0,*,*] = abs(AngleRange(c10[0,*,*]-c00[0,*,*], /pi, degrees=degrees))
	c01[0,*,*] = abs(AngleRange(c01[0,*,*]-c00[0,*,*], /pi, degrees=degrees))
	c11[0,*,*] = abs(AngleRange(c11[0,*,*]-c00[0,*,*], /pi, degrees=degrees))
	c00[0,*,*] = 0

	; Calculate maximum phase difference and maximum polar angle

	cmax = c00 > c10 > c01 > c11

	; Max. allowed polar angle, and max. allowed phase difference
	; This should reject boxes at the outside limits of the fish-eye plot

	rmax = 0.9*max(!x.crange)/(2*sqrt(2.0))*skyedge
	pmax = 90/dpm
		
	BinOK AND= cmax[0,*,*] LE pmax AND cmax[1,*,*] LE rmax

	;rpm = 1.0							; Needed for polyfill command

ENDIF ELSE IF skyedge LT 0 THEN BEGIN		; Hammer-Aitoff plot

	Xp  = HammerAitoff(Xp , dabg=dabg, degrees=degrees)
	IF has_mask THEN XpC = HammerAitoff(XpC, dabg=dabg, degrees=degrees)

	;rpm = 1.0							; Needed for polyfill command

; ENDIF ELSE BEGIN
; 	Why isn't dabg applied for Mercator projections?????????
; 	Xp  = MercatorProj(Xp , dabg=dabg, degrees=degrees)
; 	IF has_mask THEN XpC = MercatorProj(XpC, dabg=dabg, degrees=degrees)
; ENDELSE

ENDIF

FOR i=0L,N[2]-1 DO BEGIN
	i1 = i+1
	FOR j=0L,N[1]-1 DO BEGIN
		IF BinOK[j,i] THEN BEGIN
			j1 = j+1
			IF Color[j,i] LT 0 THEN col = !p.color ELSE col = Color[j,i]

			; The division by rpm is really needed !!!

			polyfill,	[Xp[0,j,i],Xp[0,j1,i],Xp[0,j1,i1],Xp[0,j,i1]],	$
				[Xp[1,j,i],Xp[1,j1,i],Xp[1,j1,i1],Xp[1,j,i1]], col=col

			;polyfill,	[Xp[0,j,i],Xp[0,j1,i],Xp[0,j1,i1],Xp[0,j,i1]]/rpm,	$
			;	[Xp[1,j,i],Xp[1,j1,i],Xp[1,j1,i1],Xp[1,j,i1]]/rpm, col=col
		ENDIF
	ENDFOR
ENDFOR

IF has_mask THEN BEGIN
	CASE use_mask[0] LT 0 OF
	0: contour, /overplot, mask     , XpC[0,*,*], XpC[1,*,*], levels=    use_mask , _extra=_extra
	1: contour, /overplot, mask NE 0, XpC[0,*,*], XpC[1,*,*], levels=abs(use_mask), _extra=_extra
	ENDCASE
ENDIF

RETURN  &  END