;+
; NAME:
;	PlotEloTimeMap
; PURPOSE:
;	Plot elongation - time map (J-map)
; CATEGORY:
;	sat/idl/toolbox
; CALLING SEQUENCE:
	PRO PlotEloTimeMap, pa_angle, time_grid, elo_grid, F, $
		degrees 	= degrees		, $
		noerase 	= noerase		, $
		breakval	= breakval		, $
		basebreak	= basebreak 	, $
		log 		= log			, $
		charsize	= charsize		, $
		title		= title 		, $
		scale		= scale 		, $
		format		= format		, $
		_extra		= _extra		, $
		silent		= silent		, $
		goodcolor	= goodcolor 	, $
		badcolor	= badcolor		, $

		pa_position = pa_position	, $
		user_position=user_position , $
		user_align	= user_align	, $
		user_string = user_string	, $

		body		= body			, $

		position	= position		, $
		legend		= legend		, $
		compass 	= compass		, $
		max_dec 	= max_dec		, $
		ra_step 	= ra_step		, $
		dec_step	= dec_step		, $
		naked		= naked
; INPUTS:
;	pa_angle	scalar; type: float
;					position angle
;	time_grid	array[n]; type: time structure
;					array of times
;	elo_grid	array[m]; type: float
;					array of solar elongations
;	F			array[n,m]; type: float
;					array of function values in "time vs elo" map;
;					each function value refers to a bin on the sky.
;					The edges of the skybin are specified in RA and DEC
; OPTIONAL INPUT PARAMETERS:
;	breakval	array[*]; type: integer or float
;					levels between colors (passed to ColorSkybox)
;					if not set, a set of break values is calculated
;					equally spaced between minimum and maximum
;	/log		if set, changes to logarithmic scale
;	/degrees	if set, indicates that all angles are in degrees
;					Default: radians.
;	scale=scale	scalar; type: float; default: 1.0
;					controls the overall size of the skymap relative
;					to the plot window
;
;	/naked		plots the skymap without any labeling or axes.
;
;	goodcolor = goodcolor
;				scalar; type: integer; default: !d.n_colors-1
;	badcolor  = badcolor
;				scalar; type: integer; default: !p.color
;
;	/compass
;	compass=compass
;				/compass will add label 'E' to left, and 'W' to
;					right side of plot. To customize labeling specify
;					2-element string (i.e. /compass is the same as
;					compass=['E','W'])
;	sn_position=sn_position
;				array[4]; type: real
;					Pairs of x,y coordinates for shifts of South/North
;					labels. Units are percentages of the window size.
;	we_position=we_position
;				array[4]; type: real
;					Pairs of x,y coordinates for shifts of East/West
;					labels.
;	pa_position=pa_position
;				array[2]; type: real
;					x,y coordinates for shifts of position angle label
;
;	user_string=user_string
;				scalar or array; type: string
;					User specified string(s) to be plotted
;	user_position=user_position
;				array[2,n]; type: real; default: [[0.05,0.95],[0.80,0.05]]
;					Start position(s) of user string(s) in normal
;					coordinates. The default allows plotting of two
;					user-defined strings in upper-left and lower-right
;					corners.
;
; OUTPUTS:
; INCLUDE:
	@compile_opt.pro		; On error, return to caller
	@vu_fnc_include.pro 	;
; CALLS:
;	ToRadians, ToDegrees, ColorSkybox, GetColors, big_eph, CvSky
;	FishEye, HammerAitoff, PlotCurve, AngleRange, gridgen
;	GeographicInfo, TimeGet, TimeUnit, ANiceNum, IsType
;	EulerRotate, MercatorProj, InitVar
;	sphere_distance
; PROCEDURE:
; MODIFICATION HISTORY:
;	SEP-1999, Paul Hick (UCSD/CASS)
;-

rpm	= ToRadians(degrees=Degrees)
dpm	= ToDegrees(degrees=Degrees)
pi	= !dpi/rpm
pi2	= pi/2
pi4	= pi/4

skymap = IsType(F, /defined)			; Check whether skymap is specified
InitVar, format, '(F4.1)'

InitVar, silent  , 0
InitVar, charsize, 1
InitVar, legend	 , /key
InitVar, naked	 , /key

IF NOT skymap THEN BEGIN
	message, 'no sky map info specified'
	;message, /info, 'no sky map or point source info specified'
	;RETURN
ENDIF

hammer	 =  1
fish	 =  0
mercator =  0
skyedge  = -1

InitVar, scale	, 1.0
InitVar, log	, /key

; 'xrange' is the horizontal range in data units.

InitVar, position, [0.55,0.50]#[1,1]+0.4*[1,1]#[-1,1]

; If the window size is not set already the 'erase' command does it here

InitVar, noerase, /key
IF NOT noerase THEN erase

njd_grid = TimeGet(time_grid,/njd)
xrange = [min(njd_grid),max(njd_grid)]
yrange = [min(elo_grid),max(elo_grid)]

plot, xrange, yrange,	$
	;xrange=xrange	,	$			; Needed to flip the x-axis if mirror is set
	/nodata 		,	$
	/noerase		,	$
	position=position,	$
	xstyle=5		,	$
	ystyle=5

!p.clip[2:3] = !p.clip[2:3]+10					; Fixes bug with IDL clip window

CASE IsType(breakval, /defined) OF

0: BEGIN

	fmin = min(F, /nan, max=fmax)

	IF log THEN BEGIN
		fmin = alog10(fmin)
		fmax = alog10(fmax)
	ENDIF
	fmin = ANiceNum(fmin)
	fmax = ANiceNum(fmax, /upper)

	N = 40
	breakval = gridgen(N, /open, range=[fmin,fmax])

	IF log THEN breakval = 10^breakval

END

1: BEGIN

	IF IsType(basebreak,/defined) THEN BEGIN
		CASE basebreak OF
		'mean'  : baseval = mean  (F,/nan)
		'median': baseval = median(F)
		ENDCASE
		message, /info, basebreak+' = '+strcompress(baseval)
		breakval += baseval
	ENDIF 

END

ENDCASE

ColorEloTimeBox, njd_grid, elo_grid, $
	GetColors(F, breakval, legend=legend,	$
		format		= format	, $
		/open					, $
		charsize	= charsize	, $
		_extra		= _extra)	, $
	degrees = degrees	, $
	/black				, $
	_extra	= _extra

plot, xrange, yrange,	$
	;xrange=xrange	,	$			; Needed to flip the x-axis if mirror is set
	/nodata 		,	$
	/noerase		,	$
	position=position,	$
	xstyle=5		,	$
	ystyle=1		,	$			; Make y-axis visible (elongation)
	ytitle="Elongation",$
	charsize=charsize


TimeXAxis, TimeSet(njd=!x.crange), /exact, /axis_only, xaxis=0, charsize=charsize, _extra=_extra
whatis, !x.crange,!y.crange[1]*[1,1]

oplot, !x.crange,!y.crange[1]*[1,1]

dx = 0.01
dy = 0.01

InitVar, pa_position, [0,0]

xyouts, /normal, $
	0.99+dx*pa_position[0]	, $
	0.95+dy*pa_position[1]	, $
	'PA='+string(pa_angle,format='(F6.2)')+([' radians',' deg'])[degrees] , $
	align		= 1 		, $
	charsize	= charsize	, $
	_extra		= _extra

IF IsType(title, /defined) THEN			 $
	xyouts, /normal, 0.05, 0.01 , title ,$
		charsize	= charsize	, $
		_extra		= _extra

PlotUserstring, user_string, user_position, $
	align	= user_align	, $
	charsize= charsize		, $
	_extra=_extra


RETURN
IF naked THEN RETURN



; Overplot bodies supported by big_eph

IF IsType(body,/defined) THEN BEGIN

	pos = big_eph(UT	, $
		body   = body	, $
		center = 'earth', $
		/to_equatorial	, $
		/to_sphere		, $
		names  = body_names, $
		degrees= degrees, $
		/silent 		)

	pos = pos[0:1,*]				; RA and declination of bodies
	point_elo = sphere_distance(pos,SunEqu,degrees=degrees)
									; Solar elongation of bodies
	npos = n_elements(pos)/2		; Number of bodies

	pos = CvSky(UT, from_equatorial=pos, /to_ecliptic, degrees=degrees, /silent)

	; Set up size of circle in degrees. Then convert to data units.

	dG = SuperArray([1,0],npos,/trail)*point_sz0[0]

	CASE 1 OF
	fish	: dG = FishEye	   ( dG, degrees=degrees, maxelo=pi2 )
	hammer	: dG = HammerAitoff( dG, degrees=degrees )
	mercator: dG = MercatorProj( dG, degrees=degrees )
	ENDCASE

	dG = reform(dG[0,*])

	IF skymap THEN BEGIN

		lng = reform(pos[0,*])
		lat = reform(pos[1,*])

		lng = AngleRange(lng, degrees=degrees)

		RAs = AngleRange(RA, degrees=degrees)
		tmp = max(RAs, izero)
		RAs = shift(RAs,-izero-1)

		i = round( interpol(indgen(nRA), RAs, lng) )
		j = round( interpol(indgen(nDE), DEC, lat) )

		a = where(0 LE i AND i LE nRA-1 AND 0 LE j AND j LE nDE-1)

		; Background values in skymap at grid points nearest to the
		; point source location. Note that some of these may be bad.

		good = replicate(BadValue(0.0),npos)
		IF a[0] NE -1 THEN good[a] = F[(i[a]+izero+1) mod nRA,j[a]]

		; The circle edge is drawn using the foreground color for positive,
		; and background color for negative differences.

		good = finite(good)
		good = (1-good)*!p.color+good*!p.background

	ENDIF

	;print, projsz, dg[0], dg[0]/projsz
	CASE 1 OF
	fish    : pos = FishEye		(pos, zero_phase=zero_point, dabg=dabg, degrees=degrees, maxelo=skyedge)
	hammer  : pos = HammerAitoff(pos, zero_phase=zero_point, dabg=dabg, degrees=degrees)
	mercator: pos = MercatorProj(pos, zero_phase=zero_point, dabg=dabg, degrees=degrees)
	ENDCASE

	; For fish-eye maps plot only point sources inside skyedge.
	; Note that the elongation is calculated for the time of the skymap (NOT for
	; the time of observation of the point source).

	IF silent LE 0 THEN message, /info, 'plotting position of '+body

	FOR i=0,npos-1 DO BEGIN
		IF finite(dG[i]) AND (NOT fish OR point_elo[i] LT skyedge) THEN BEGIN
			polyfill, pos[0,i]+dG[i]*xcircle, pos[1,i]+dG[i]*ycircle,color=good[i]
			sgn = pos[0,i] LT 0
			xyouts, pos[0,i]+(1-2*sgn)*dG[i], pos[1,i], body_names[i] , $
				align	= -0.1*sgn+1.1*(1-sgn)	, $
				color	= good[i]
		ENDIF
	ENDFOR

ENDIF

InitVar, goodcolor, !d.n_colors-1
InitVar, badcolor , !p.color

CASE 1 OF
hammer OR mercator: BEGIN

	IF hammer THEN max_dec = 0.5*pi ELSE InitVar, max_dec, 0.5*pi
	max_dec = max_dec < 0.5*pi

	n   = 37
	lng = replicate(1,n)
	lat = gridgen(n, range=[-1,1])

	lng = [-lng, lng,-lng[0]]			; Longitudes relative to zero_phase
	lat = [ lat,-lat, lat[0]]

	n = n_elements(lng)

	; Plot outer boundary: +/- 180 deg meridian and horizontal lines at poles
	; for the mercator projection.

	pos = transpose( [[pi*lng],[max_dec*lat]] )
	IF hammer THEN pos = HammerAitoff(pos, degrees=degrees)
	oplot, pos[0,*], pos[1,*], _extra=_extra

	; Meridians at steps of ra_step

	InitVar, ra_step , 0.5d0*pi
	InitVar, dec_step, max_dec	;0.5d0*pi

	this = fix(pi/ra_step)
	this = -this*ra_step

	n   = 37
	lng = replicate(1,n)
	lat = max_dec*gridgen(n, range=[-1,1])

	WHILE this LE pi DO BEGIN

		pos_this = this
		IF mercator THEN pos_this = AngleRange(this-zero_point,degrees=degrees,/pi)

		IF pos_this NE 0 THEN BEGIN

			pos = transpose( [[pos_this*lng],[lat]] )

			IF skymap THEN BEGIN
				RAs = AngleRange(RA, degrees=degrees)
				tmp = max(RAs, izero)
				RAs = shift(RAs,-izero-1)

				i = round( interpol(indgen(nRA), RAs, AngleRange(reform(pos[0,*])+zero_point, degrees=degrees)) )
				j = round( interpol(indgen(nDE), DEC, reform(pos[1,*])) )
		
				a = where(0 LE i AND i LE nRA-1 AND 0 LE j AND j LE nDE-1)

				good = bytarr(n)
				IF a[0] NE -1 THEN good[a] = finite(F[(i[a]+izero+1) mod nRA,j[a]])
			ENDIF

			IF hammer THEN pos = HammerAitoff(pos, degrees=degrees)

			IF mercator OR abs(this) NE pi THEN BEGIN
				CASE goodcolor EQ badcolor OF
				0: PlotCurve, /oplot, pos[0,*], pos[1,*], good, color=[goodcolor,badcolor], /changecolor, /silent, _extra=_extra
				1: PlotCurve, /oplot, pos[0,*], pos[1,*], color=badcolor, /silent, _extra=_extra
				ENDCASE
			ENDIF

			IF mercator THEN BEGIN
				xyouts, pos_this, -max_dec-4*0.0177*projsz*charsize, charsize=charsize, _extra=_extra, $
					strcompress(round(AngleRange(this*dpm,/degrees)),/rem), alignment=0.5
				IF pos_this EQ pi THEN BEGIN
					xyouts, -pi, -max_dec-4*0.0177*projsz*charsize, charsize=charsize, _extra=_extra, $
						strcompress(round(AngleRange(this*dpm,/degrees)),/rem), alignment=0.5
				ENDIF
			ENDIF

		ENDIF

		this = this+ra_step

	ENDWHILE

	n   = 37
	lng = pi*gridgen(n, range=[-1,1])
	lat = replicate(1,n)

	this = -fix(max_dec/dec_step)*dec_step

	WHILE this LE max_dec DO BEGIN

		IF this NE 0 THEN BEGIN

			pos = transpose( [[lng],[this*lat]] )

			IF skymap THEN BEGIN
				RAs = AngleRange(RA, degrees=degrees)
				tmp = max(RAs, izero)
				RAs = shift(RAs,-izero-1)

				i = round( interpol(indgen(nRA), RAs, AngleRange(reform(pos[0,*])+zero_point, degrees=degrees)) )
				j = round( interpol(indgen(nDE), DEC, reform(pos[1,*])) )
		
				a = where(0 LE i AND i LE nRA-1 AND 0 LE j AND j LE nDE-1)

				good = bytarr(n)
				IF a[0] NE -1 THEN good[a] = finite(F[(i[a]+izero+1) mod nRA,j[a]])
			ENDIF

			IF hammer THEN pos = HammerAitoff(pos, degrees=degrees)

			CASE goodcolor EQ badcolor OF
			0: PlotCurve, /oplot, pos[0,*], pos[1,*], good, color=[goodcolor,badcolor], /changecolor, /silent, _extra=_extra
			1: PlotCurve, /oplot, pos[0,*], pos[1,*], color=badcolor, /silent, _extra=_extra
			ENDCASE

			IF mercator THEN BEGIN
				xyouts, !x.crange+0.0566*projsz*mirror*[-1,1]*charsize, [this,this]-0.0177*projsz*charsize,	$
					strcompress(round(this*dpm),/rem), alignment=0.5, charsize=charsize, _extra=_extra
			ENDIF

		ENDIF

		this = this+dec_step

	ENDWHILE

END
fish: BEGIN

	;origin = [440,399.5]
	;xy = gridgen([!d.x_size,!d.y_size], origin=origin)
	;pp = cv_coord(from_rect=xy, /to_polar)
	;pp = where( 310 le pp[1,*] and pp[1,*] le 330 )
	;for i=0L,n_elements(pp)-1 do tv, [!p.background], origin[0]+xy[0,pp[i]], origin[1]+xy[1,pp[i]]

	f = projsz/skyedge
	FOR a=pi4,pi2,pi4 DO oplot, f*a*xcircle,f*a*ycircle, color=badcolor, _extra=_extra	; draw 45 and 90deg circle

	tmp = (convert_coord(f*[ [-pi2,0],[pi2,0],[0,-pi4],[0,pi4] ], /data, /to_normal))[0:1,*]

	IF skyedge GE pi2 THEN BEGIN
		xyouts, /normal, tmp[0]-0.005, tmp[1]+0.005, '90',		$
			charsize=charsize, color=badcolor, align=1, _extra=_extra
		xyouts, /normal, tmp[2]+0.005, tmp[3]+0.005, '90',		$
			charsize=charsize, color=badcolor, align=0, _extra=_extra
	ENDIF

	IF skyedge GE pi4 THEN BEGIN
		xyouts, /normal, tmp[4]+0.005, tmp[5]+0.005, '45',		$
			charsize=charsize, color=badcolor, _extra=_extra
		xyouts, /normal, tmp[6]+0.005, tmp[7]+0.005, '45',		$
			charsize=charsize, color=badcolor, _extra=_extra
	ENDIF

	;da = 5/dpm
	;IF skyedge GE pi2 THEN FOR a=-1, 1,2 DO xyouts, a*(pi2+da)*f, -da*f, '90', alignment=0.5, charsize=charsize, color=badcolor, _extra=_extra
	;IF skyedge GE pi4 THEN FOR a=-1,-1,2 DO xyouts, -da*f, a*(pi4+da)*f, '45', alignment=0.5, charsize=charsize, color=badcolor, _extra=_extra

END
ENDCASE


CASE 1 OF
fish	: tmp = FishEye	    ( [[-pi2,0],[pi2,0],[0,-pi2],[0,pi2]], maxelo=pi2, degrees=degrees )
hammer	: tmp = HammerAitoff( [[-pi ,0],[pi ,0],[0,-max_dec],[0,max_dec]], degrees=degrees )
mercator: tmp = [[-pi ,0],[pi ,0],[0,-max_dec],[0,max_dec]]
ENDCASE

oplot, [tmp[0,0],tmp[0,1]],[tmp[1,0],tmp[1,1]], color=badcolor, _extra=_extra	; Draw horizontal axis
oplot, [tmp[0,2],tmp[0,3]],[tmp[1,2],tmp[1,3]], color=badcolor, _extra=_extra	; Draw vertical axis

tmp = (convert_coord(tmp,/data,/to_normal))[0:1,*]

; Label with directions (N,S always; E,W optional)

dx = 0.01
dy = 0.01

InitVar, sn_position, [0,0,0,0]

xyouts, /normal, tmp[0,2:3]+dx*sn_position[[0,2]], tmp[1,2:3]+dy*([-4,2]+sn_position[[1,3]]),	$
	(['Ecl ','Eq '])[0]+['S','N'], alignment=0.5, charsize=charsize, _extra=_extra

CASE 1 OF
IsType(compass,/string ): we_string = compass
IsType(compass,/defined): IF compass THEN we_string = ['W','E']
ELSE:
ENDCASE

IF IsType(we_string,/defined) THEN BEGIN
	InitVar, we_position, [0,0,0,0]
	xyouts, /normal, alignment=0.5, charsize=charsize, _extra=_extra,		$
		tmp[0,0:1]+dx*([ 3,-3]+we_position[[0,2]]),	$
		tmp[1,0:1]+dy*([-1,-1]+we_position[[1,3]]),	$
		we_string
ENDIF

RETURN  &  END