;+
; NAME:
;	ColorPolarBox
; 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=polar), a fish-eye map or a Hammer-Aitoff map
; CALLING SEQUENCE:
	PRO ColorPolarBox, cosP, sinP, Color,	$
		skyedge		= skyedge	, $
		zero_phase	= zero_phase, $
		dabg		= dabg		, $
		degrees		= degrees	, $
		black		= black		, $
		void		= void		, $
		fill2edge	= fill2edge
; 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
;	Polar		array[m [+1] ] or array[n [+1],m [+1] ]; type: float
;					Polar angle; 'colatitude'
;					Box center of 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 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
; OPTIONAL INPUTS:
;	zero_phase=zero_phase
;				scalar, or array with same structure as 'Phase'; type: float
;					Only used if Xp is a 1-dim array of box centers.
;					The input arrays Xp, Yp, 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 polar angle dimension
; >	The Phase and Polar 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 Polar 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 Polar should be the same (m=M) or one larger (m=M+1) than for Color.
;	If m=M+1 then Polar contains the polar angles of the edges of sky boxes; if m=M it contains
;	the center polar 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. 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 Polar arrays
; >	Phase and Polar are spherical coordinates in a coordinate system where
;	the reference direction (e.g. the direction to the Sun) is the Z-axis.
; >	Synonyms:
;		phase angle = longitude, position angle
;		polar angle = colatitude, elongation
; >	The arrays Phase and Polar 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 Polar are plotted
;	on the screen as a 'fish-eye' view, interpreting the Polar angle
;	as 'radius' and the Phase angle as an phase angle in the x-y plane of
;	the plot, i.e. x = Polar*cos(Phase) and y = Polar*sin(Phase).
; >	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; pphick@ucsd.edu)
;		Added /fill2edge keyword.
;-

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

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

Xp = cosP
Yp = sinP

szX = size(Xp)
szY = size(Yp)

print, zero_phase
;IF szX[0] EQ 1 AND IsType(zero_phase, /defined) THEN BEGIN
;	IF szX[1] EQ N[1] THEN BEGIN		; Box centers specified for phase angle
;		Xp = 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 szY[0] EQ 2 THEN IF szY[1] EQ N[1] THEN Yp = Yp[dp,*]
;	ENDIF
;ENDIF

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

szX = size(Xp)
szY = size(Yp)

M = szX[1]							; # polar 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 separation of centers

	; 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 ]

	szX = size(Xp)
	szY = size(Yp)
ENDIF


M = szY[2]							; # declinations
IF M EQ N[2] THEN BEGIN				; Box centers specified for polar 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)

	szX = size(Xp)
	szY = size(Yp)
ENDIF

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

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

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

; 	Xp = MercatorProj(Xp, dabg=dabg, degrees=degrees)


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

RETURN  &  END