;+
; NAME:
;	PlotSynopticMap
; PURPOSE:
;	Plot synoptic map
; CATEGORY:
;	WWW: graphics
; CALLING SEQUENCE:
	PRO PlotSynopticMap, XCplot, Map, $	
		XCmap	= XCmap		, $
		YLmap	= YLmap		, $
		YLshow	= YLshow	, $
		YLplot	= YLplot	, $
		degrees	= degrees	, $
		breakval= breakval	, $
		plotcenter=plotcenter,$
		plotsize= plotsize	, $
		fill	= fill		, $
		grid	= grid		, $
		draw_contour=draw_contour, $
		timeaxis= timeaxis	, $
		user_position=user_position , $
		user_align	= user_align	, $
		user_string = user_string	, $
		_extra	= _extra
; INPUTS:
;	XCplot			scalar or array[2]; type: float
;						start and end Carrington variable to be displayed
;						if a scalar is specified this is used as start value, and the end
;						value is set to the start value, plus one (i.e. one rotation
;						is assumed)
;	Map				2D array; any type
;
; OPTIONAL INPUT PARAMETERS:
;	xcmap=XCmap		array[2]; type: float; default: XCplot
;						range of Carrington variable of input array Map
;	ylmap=YLmap		array[2]; type: float; default: [-90,90] degrees
;						latitude range of input array Map
;	ylplot=YLplot	array[2]; type: float: default: ylmap
;						latitude range covered by vertical axis
;	ylshow=YLshow	array[2]; type: float: default: ylplot
;						latitude range of input map actually shown on plot
;	/degrees		scalar; type: integer; default:0 (radians)
;						if set YLmap should be in degrees
;	plotcenter=plotcenter
;					array[2]; type: integer; default: [!d.x_size,!d.y_size]/2
;						pixel coordinates for center of map
;	plotsize=plotsize
;					array; type: integer; default: 0.8*plotcenter[0]
;						size of map in pixels
;
;	breakval=breakval
;					array[*]; type: float
;						contour levels used to plot map
;	/draw_contour	if set then the IDL 'contour' function is used to plot
;						the map instead of the 'tv' function
;
;	/fill			if set then 'map_grid' is called to fill in bad values
;						before plotting the map
;	/grid			if set then the grid of the input map is marked by
;						putting crosses at the corners of grid squares.
;	timeaxis=timeaxis	if set, adds a time axis at the top of the map, specifying the time at
;				    which the corresponding heliographic longitude crossed the
;				    central meridian on the solar disk
;
;	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:
;	(none)
; OPTIONAL OUTPUT PARAMETERS:
; INCLUDE:
	@compile_opt.pro			; On error, return to caller
; CALLS:
;	InitVar, ToRadians, ToDegrees, Carrington, TimeXAxis, GetColors
;	gridfill, gridgen, IsType, TimeUnit
; PROCEDURE:
; MODIFICATION HISTORY:
;	APR-1999, Paul Hick (UCSD/CASS)
;	FEB-2002, Paul Hick (UCSD/CASS)
;		Added xsize and ysize keywords to 'tv' call. This is needed for
;		devices with scalable pixels, such as PS.
;	JUL-2002, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;		Added 'timeaxis' keyword to add a time axis
;		at top of map. This was always done before.
;-

InitVar, draw_contour, /key
InitVar, fill		 , /key
InitVar, grid		 , /key
InitVar, timeaxis	 , /key

CASE n_elements(plotcenter) OF
0   : plotcenter = [!d.x_size,!d.y_size]/2
1   : plotcenter = plotcenter*[1,1]
ELSE: plotcenter = plotcenter[0:1]
ENDCASE

InitVar, plotsize, 0.8*plotcenter[0]

plotcenter = round(plotcenter)
plotsize   = round(plotsize  )

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

InitVar, XCplot, XCplot[0]+[0,1], count=1
InitVar, XCmap , XCplot
InitVar, YLplot, pi2*[-1,1]
InitVar, YLmap , pi2*[-1,1]
InitVar, YLshow, YLmap

XCshow = [XCmap [0] > XCplot[0], XCmap [1] < XCplot[1]]
YLshow = [YLshow[0] > YLmap [0], YLshow[1] < YLmap [1]]
YLshow = [YLshow[0] > YLplot[0], YLshow[1] < YLplot[1]]

; Get times for start, end and center of map

T  = Carrington([XCplot,mean(XCplot),XCmap,mean(XCmap)])

; The tick interval on the horizontal axis is 30 degrees for 1 rotation
; and increases proportional to the number of rotations plotted.
; The interval is always a multiple of 15.

dXtick = 30.0
IF !d.x_size LE 400 THEN dXtick = 60.0

dXtick = round((XCplot[1]-XCplot[0])*dXtick/15)*15

dX = dXtick/360.		; Tickmarks every 30 degree

ICtv = floor(XCplot)
XCtv = XCplot-ICtv
XCtv = ICtv+[ceil(XCtv[0]/dX),floor(XCtv[1]/dX)]*dX
XCtv = XCtv[0]+dX*(findgen(1+round((XCtv[1]-XCtv[0])/dX)))
XCtn = (360-round((XCtv-floor(XCtv))*360)) mod 360
XCtn = strcompress(XCtn,/rem)

YLtv = gridgen(7, range=pi2*[-1,1])
YLtv = YLtv[ where(YLplot[0] LE YLtv AND YLtv LE YLplot[1]) ]

; The displayed map has is a rectangle of 'tmp' pixels centered on 'plotcenter'

tmp = plotsize*[2,1]/2
corner = [ [plotcenter-tmp], [plotcenter+tmp] ]
corner = corner/2*2+1				; Corners of the plot area (in pixels)

; The plot area covered by 'corner' (in pixels) corresponds to 'XCplot'.
; 'XCmap' specifies the Carrington range of the input map. Converted to pixels this becomes x.

x = corner[0,1]+(XCmap-XCplot[0])/(XCplot[1]-XCplot[0])*(corner[0,0]-corner[0,1])
y = corner[1,0]+(YLmap-YLplot[0])/(YLplot[1]-YLplot[0])*(corner[1,1]-corner[1,0])

tmp = x-round(x)
if abs(tmp[0]) LE abs(tmp[1]) THEN tmp = tmp[0] ELSE tmp = tmp[1]
x = round(x-tmp)

y = round(y)

; The lower-left corner of the input map is put at pixel 'mapcorner'

mapcorner = [ [x[1],y[0]], [x[0],y[1]] ]
mapsize   = [ mapcorner[0,1]-mapcorner[0,0] , mapcorner[1,1]-mapcorner[1,0] ]+1

showcorner = mapcorner

; Calculate range of pixels in vertical direction corresponding to YLshow
; and adjust showcorner accordingly

showcorner[1,*] = showcorner[1,0]+round( (YLshow-YLmap[0])/(YLmap[1]-YLmap[0])*(mapsize[1]-1) )

; Determine corners of map inside the plot area

showcorner = [ [ showcorner[*,0] > corner[*,0] ], [ showcorner[*,1] < corner[*,1] ] ]

erase

; convert_coord needs an open window. This is created at the latest with the erase statement

position = reform((convert_coord(corner,/device,/to_normal))[0:1,*],4)

; Check whether part of the synoptic map is still inside the plot area

plot_visible = showcorner[0,0] LE showcorner[0,1] AND showcorner[1,0] LE showcorner[1,1]

IF plot_visible THEN BEGIN

	Bad = where( 1B-finite(Map) )

	CASE fill OF
	0: MapWin = Map
	1: MapWin = GridFill(Map)
	ENDCASE

	MapWin = congrid(MapWin,mapsize[0],mapsize[1],/interp,/minus_one)

	; Extract part of MapWin inside plot area

	tmp = showcorner-[mapcorner[*,0],mapcorner[*,0]]
	MapWin = MapWin[tmp[0,0]:tmp[0,1],tmp[1,0]:tmp[1,1]]
	ndim = tmp[*,1]-tmp[*,0]+1

ENDIF

; If image is plotted using the IDL tv function we put the image on the screen
; prior to plotting all the other annotation (i.e. prior to establishing a
; user scale).

IF NOT draw_contour THEN BEGIN

	IF plot_visible THEN BEGIN

		CASE IsType(breakval, /defined) OF
		0: MapWin = bytscl(MapWin, _extra=_extra, /nan)
		1: MapWin = GetColors(MapWin, breakval, /open, /badbackground, /legend, _extra=_extra)
		ENDCASE

		; The xsize and ysize keywords are needed for PostScript output.

		tv, MapWin, showcorner[0,0],showcorner[1,0], xsize=ndim[0], ysize=ndim[1]

	ENDIF ELSE IF IsType(breakval, /defined) THEN	$

		tmp = GetColors(0, breakval, /open, /badbackground, /legend, _extra=_extra)

ENDIF

plot, XCplot, YLplot,	$
	/nodata,			$
	/noerase,			$
	xrange	= reverse(XCplot),	 $
	xstyle	= 9,		$
	ystyle	= 9,		$
	position= position,	$
	xticks	= n_elements(XCtv)-1,$
	xtickv	= XCtv,		$
	xtickn	= XCtn,		$
	yticks	= n_elements(YLtv)-1,$
	ytickv	= YLtv,		$
	ytickn	= round(abs(YLtv)*dpm),$ ; Labels should always be in degrees
	_extra	= _extra


; If line contours are drawn this is done after the axis have been plotted
; (i.e. after the user scale has been established).

IF draw_contour AND plot_visible THEN	$

	contour, MapWin,							$
		gridgen(ndim[0], range=reverse(XCshow)),$
		gridgen(ndim[1], range=YLshow),			$
		levels	= breakval	, $
		position= position	, $
		/overplot

	; Draw vertical line at start of new rotation (0 deg longitude)

FOR i=ceil(XCplot[0]),floor(XCplot[1]) DO oplot, [i,i], !y.crange
oplot, XCplot, [0,0]

; Draw vertical line at right edge of plot

; IDL BUG:
; The oplot command sometimes does not draw anything, because
; IDL sometimes sets !p.clip just inside the right edge of the plot.
; As a fix !p.clip is explicitly increased by one pixel on the right

!p.clip[2:3] = !p.clip[2:3]+1

oplot, XCplot[0]*[1,1],!y.crange

; Plot Carrington rotation below x-axis

FOR i=floor(XCplot[0]),floor(XCplot[1]) DO BEGIN
	tmp = mean([i > XCplot[0], i+1 < XCplot[1]])
	IF  XCplot[0]+0.1 LT tmp AND tmp LT XCplot[1]-0.1 THEN	$
		xyouts, tmp, total(([0,1]+1.2*[1,-1])*!y.crange), 'CR '+strcompress(i,/rem), align=0.5, _extra=_extra
ENDFOR

CASE timeaxis OF
0: oplot, XCplot, !y.crange[1]*[1,1]
1: TimeXAxis, /axis_only, /exact, T[[1,0]], /noxtitle, upto=TimeUnit(/day), _extra=_extra
ENDCASE

IF grid THEN BEGIN
	tmp = size(Map, /dim)
	x = gridgen(tmp[0], range=XCMap)
	y = gridgen(tmp[1], range=YLMap)
	FOR i=0,tmp[0]-1 DO			$
		FOR j=0,tmp[1]-1 DO		$
			IF XCplot[0] LE x[i] AND x[i] LE XCplot[1] THEN plots, x[i], y[j], psym=3
ENDIF

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

RETURN  &  END