;+
; NAME:
;	PlotSphereCut
; PURPOSE:
;	Plot lng-lat map of sun-centered sphere
; CATEGORY:
;	sat/idl/toolbox/graphics
; CALLING SEQUENCE:
	PRO PlotSphereCut, Map, ut=ut	, $
		radius		= radius		, $
		euler_angles= euler_angles	, $
		degrees 	= degrees		, $
		euler_info	= euler_info	, $
		breakval	= breakval		, $
		title		= title 		, $
		upto		= upto			, $
		body		= body			, $
		user_position=user_position , $
		user_align	= user_align	, $
		user_string = user_string	, $

		plotcenter	= plotcenter	, $
		plotsize	= plotsize		, $
		grid		= grid			, $

		_extra=_extra
; INPUTS:
;	zz				array[n,n]; type: float
;						function value on disk
;	ut=ut			array[1]; type: time structure
;						UT time used for ephemeris calculation
; OPTIONAL INPUT PARAMETERS:
;	upto=upto		scalar; type: integer; default: none
;						return value of href=TimeUnit= function.
;						The 'ut' time is plotted in the upper right
;						corner using TimeGet. The 'upto' keyword
;						determines at which time unit the string
;						is terminated.
;	radius=radius	scalar; type: float; default: 1.0
;						radius of sphere (in AU)
;	euler_angles=euler_angles
;	/degrees		if set, angles are in degrees
;						default is radians
;	euler_info=euler_info
;	title=title 	scalar; type: string
;						string plotted in lower-left corner
;	breakval=breakval
;					array; type: integer or float
;						contour levels passed to GetColors
;
;	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.
;
;	_extra=_extra	additional plot keywords
; OUTPUTS:
;	(plot)
; INCLUDE:
	@compile_opt.pro		; On error, return to caller
; CALLS:
;	InitVar, IsType, setup3d, IsType, GetColors, plot3darc, gridgen
;	big_eph, big_orbit, EulerRotate, TimeGet
; PROCEDURE:
;	The appearance of the map is controlled by specifying
;	plot keywords (passed here through the _extra keyword):
;
;	Control character size with keyword charsize
;	Control character thickness with keyword charthick
;	Control size of body symbol with keyword symsize
;	Control thickness of orbit with keyword thick
; MODIFICATION HISTORY:
;	MAY-2010, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;-

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

; Set center of map in device coordinates

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

; Set size of map in device coordinates

InitVar, plotsize, 0.8*plotcenter[0]

; Round center and size to integer

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

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

XLplot = [0,2*pi]

InitVar, XLmap , XLplot
InitVar, YLplot, pi2*[-1,1]
InitVar, YLmap , pi2*[-1,1]
InitVar, YLshow, YLmap

XLshow = [XLmap [0] > XLplot[0], XLmap [1] < XLplot[1]]
YLshow = [YLshow[0] > YLmap [0], YLshow[1] < YLmap [1]]
YLshow = [YLshow[0] > YLplot[0], YLshow[1] < YLplot[1]]

; The tick interval on the horizontal axis is 30 degrees

XLtv = gridgen(13, range=[0,2*pi])
XLtv = XLtv[ where(XLplot[0] LE XLtv AND XLtv LE XLplot[1]) ]

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 'XLplot'.
; 'XLmap' specifies the range of the input map. Converted to pixels this becomes x.

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

x = round(x)
y = round(y)

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

mapcorner = [ [x[0],y[0]], [x[1],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 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, XLplot, YLplot,	$
	/nodata,			$
	/noerase,			$
	;xrange	= XLplot,	$
	position= position,	$

	xstyle	= 9,		$
	xticks	= n_elements(XLtv)-1,$
	xtickv	= XLtv,		$
	xtickn	= strcompress(round(abs(XLtv)*dpm),/rem),$

	ystyle	= 9,		$
	yticks	= n_elements(YLtv)-1,$
	ytickv	= YLtv,		$
	ytickn	= strcompress(round(YLtv*dpm),/rem),$ ; Labels should always be in degrees
	;ytickn	= strcompress(round(abs(YLtv)*dpm),/rem),$ ; 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(XLshow)),$
		gridgen(ndim[1], range=YLshow),			$
		levels	= breakval	, $
		position= position	, $
		/overplot

oplot, XLplot, [0,0]

;rr = fltarr(2,360)


; Label UT time

xyouts, 0.99, 0.95, /norm, align=1, _extra=_extra,	$
	TimeGet(ut, /ymd, upto=upto, /roundt)+' UT'

; Label title

IF IsType(title, /defined) NE 0 THEN	$
	xyouts, 0.05, 0.01, /norm, _extra=_extra, title

; 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, XLplot[1]*[1,1],!y.crange
oplot, XLplot, !y.crange[1]*[1,1]

IF grid THEN BEGIN
	tmp = size(Map, /dim)

	x = gridgen(tmp[0], range=XLMap)
	y = gridgen(tmp[1], range=YLMap)

	FOR i=0,tmp[0]-1 DO			$
		FOR j=0,tmp[1]-1 DO		$
			IF XLplot[0] LE x[i] AND x[i] LE XLplot[1] THEN plots, x[i], y[j]
ENDIF

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

RETURN  &  END