;+
; NAME:
;	PlotPolarSkymap
; PURPOSE:
;	Plot sky map in polar projection, with locations of point sources overplotted (optional)
; CATEGORY:
;	gen/idl/graphics
; CALLING SEQUENCE:
	PRO PlotPolarSkymap, UT, cosRA, sinRA, F,	$
		equator 	= equator	, $
		geolng		= geolng	, $
		degrees 	= degrees	, $
		breakval	= BreakVal	, $
		logscale	= logscale	, $
		charsize	= charsize	, $
		title		= title 	, $
		upto		= upto		, $
		scale		= scale 	, $
		minelo		= minelo	, $
		zero_point	= zero_point, $
		zero_phase	= zero_phase, $
		dabg		= dabg		, $
		format		= format	, $
		fill2edge	= fill2edge , $
		_extra		= _extra	, $
		galactic	= galactic	, $
		silent		= silent	, $

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

		point_sources=point_sources,$
		point_names = point_names  ,$
		point_size	= point_size   ,$
		point_onesize=point_onesize,$

		mirror		= mirror	, $
		noerase 	= noerase	, $
		xrange		= xrange	, $
		yrange		= yrange	, $
		position	= position	, $
		legend		= legend	, $
		dec_step	= dec_step	, $
		ra_step 	= ra_step	, $

		ra_position = ra_position, $
		naked		= naked
; INPUTS:
;	UT			scalar; type: standard time structure
;					time; determines position of Sun in the sky
;	cosRA		array[n]; type: float
;				+/- (90-DEC)*cos(RA)
;					array of ecliptic longitudes or right ascensions
;					(if /equator set) in the range [-180,+180]??NOT REALLY??
;	sinRA		array[m]; type: float
;				+/- (90-DEC)*sin(RA)
;					array of ecliptic latitudes or declinations
;					(if /equator is set) in the range [-90.90]
;	F			array[n,m]; type: float
;					array of function values in sky map; each function value
;					refers to 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
;	/logscale	if set, changes to logarithmic scale
;	/degrees	if set, indicates that all angles are in degrees
;					Default: radians.
;	minelo		scalar; type: float
;					Erases plot inside of minelo
;	/equator	if set, the 'horizontal' in the sky map is parallel to the
;					equatorial plane. By default the horizontal is along the ecliptic
;	upto=upto	scalar; type: integer; default: TimeUnit(/hour)
;					controls the length of the UT string plotted
;	scale=scale	scalar; type: float; default: 1.0
;					controls the overall size of the skymap relative to the plot window
;	/galactic	if set, add a line for the galactic plane
;
;	/naked		if set, then plot skymap without labeling and axes.
;
;	point_sources=point_sources
;				array; type: sky_point_source structure
;					contains information about point sources to be overplotted
;					on the skymap. See href=vu_point_source=.
;	/point_names
;				if set, then the names of the point sources are plotted
;
;	point_size=point_size
;				scalar or array[2]; controls the size of the circles used to plot the
;				point source position. The size is specified in radians or degrees
;				(depending on setting of /degrees).
;
;				point_size[0]: minimim size of circle.
;				point_size[1]: increment in circle size per unit of the function
;					value in skymap F.
;
;				point_size[1] is used to increase the circle size proportional to
;				the difference between the point source fnc value and the value in
;				the skymap F.
;
;				If point_size[1] is not set bad (and /point_onesize is NOT set)
;				then only sources with valid function values will be plotted.
;				(I.e. to plot only good sources with the same cirle size set
;				point_size[1] to zero.)
;
;				If point_size[1] is bad or omitted then /point_onesize is assumed set.
;
;	/point_onesize
;				if set then all sources (including those with bad fnc values)
;				are plotted with the same size circle point_size[0].
;				This is useful to override the information stored
;				in 'point_sources' or in keyword element point_size[1].

;				If only IPS sources are plotted (i.e. if no skymap F is specified) then
;				setting /point_onesize will plot all source position with the same size.
;
;	The following three keywords are clumsy to use. The idea is to control the orientation
;	of the map (i.e. direction of origin, and tilt of horizontal on sky).
;	href=vu_earthskymap= uses these keywords to set up sky 'snapshots' and sky 'sweeps'.
;
;	dabg		array[3]; type: float; default: [0,0,0]
;					(used by the projection functions FishEye, HammerAitoff and
;					MercatorProj, either directly in this procedure or indirectly in
;					ColorSkybox)
;
;	zero_point	scalar; type: float; default: 0.0
;					(used internally and passed to FishEye, HammerAitoff and MercatorProj)
;					Defines the longitude or RA of the origin.
;
;	zero_phase	scalar, or array with same structure as 'RA'; type: float; default: zero_point
;					(passed to ColorSkybox)
;					ColorSkybox uses this keyword to rearrange its first three arguments
;					(RA,DEC and fnc-values) to put zero_phase in the center of the map.
;					This modifies input arrays RA and DEC, but not F.
;
;	For sky 'snapshots' zero_phase is the same as zero_point.
;	For sky 'sweeps' zero_phase is a monotonic array with a value for each sweep, and
;	the value of 'zero_point' (almost) exactly in the center of 'zero_phase'. Note that
;	inconsistent settings of zero_point and zero_phase make an incorrect map.
;
;	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
; CALLS:
;	ToRadians, ToDegrees, ColorSkybox, GetColors, CvSky, FishEye, HammerAitoff
;	PlotCurve, AngleRange, gridgen, GeographicInfo, TimeGet, TimeUnit, ANiceNum
;	IsType, EulerRotate, MercatorProj, InitVar, vu_point_source, BadValue
;	big_eph, jpl_body
; SEE AlSO:
;	ips_point_source__define
; PROCEDURE:
; >	Time UT is needed to calculate the location of equator and ecliptic
; >	The ecliptic is drawn on an equatorial skymap and v.v
; MODIFICATION HISTORY:
;	SEP-1999, Paul Hick (UCSD/CASS)
;	FEB-2002, Paul Hick (UCSD/CASS)
;		Added option to plot a sky map in Mercator projection by setting maxelo=0.
;	APR-2002, Paul Hick (UCSD/CASS)
;		Added /fill2edge keyword (passed to ColorSkybox)
;	SEP-2003, Paul Hick (UCSD/CASS)
;		Minor tweaking of labels to get Mercator projection to look better
;		Added /galactic keyword.
;	AUG-2004, Paul Hick (UCSD/CASS)
;		Started adding some code to deal with Thomson scattering brightness
;		(mainly for overplotting point sources).
;		The time plotted at the top is now rounded instead of truncating
;		to timeunit 'upto'.
;	SEP-2004, Paul Hick (UCSD/CASS)
;		For fisheye plots point sources with outside maxelo are not plotted anymore
;		Reworked the determination of the circle sizes for the point sources, and
;		added some keywords to control them.
;	OCT-2006, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;		Replaced NewcombSun by big_eph calls
;-

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
pnt_n	= n_elements(point_sources)		; Check whether any point sources are specified

InitVar, point_names  , /key
InitVar, point_onesize, /key
InitVar, format 	  , '(F4.1)'
InitVar, legend		  , /key
InitVar, naked		  , /key

InitVar, galactic	, /key
InitVar, zero_point , 0.0
InitVar, silent 	, 0

ncosRA = n_elements(cosRA)
nsinRA = n_elements(sinRA)


IF pnt_n NE 0 THEN BEGIN

	CASE 1 OF
	IsType(title, /undefined): BEGIN
		pos_fnc = (where( tag_names( point_sources[0] ) EQ 'VV'   ))[0]
		pos_sz  = (where( tag_names( point_sources[0] ) EQ 'SZVV' ))[0]
	END
	strmid(title,0,9) EQ 'IPS speed': BEGIN
		pos_fnc = (where( tag_names( point_sources[0] ) EQ 'VV' ))[0]
		pos_sz  = (where( tag_names( point_sources[0] ) EQ 'SZVV' ))[0]
	END
	strmid(title,0,11) EQ 'IPS g-level': BEGIN
		pos_fnc = (where( tag_names( point_sources[0] ) EQ 'GG' ))[0]
		pos_sz  = (where( tag_names( point_sources[0] ) EQ 'SZGG' ))[0]
	END
	strmid(title,0,11) EQ 'Thomson sca': BEGIN
		pos_fnc = (where( tag_names( point_sources[0] ) EQ 'BB' ))[0]
		pos_sz  = (where( tag_names( point_sources[0] ) EQ 'SZBB' ))[0]
	END
	strmid(title,0,10) EQ 'Thomson pB': BEGIN
		pos_fnc = (where( tag_names( point_sources[0] ) EQ 'PB' ))[0]
		pos_sz  = (where( tag_names( point_sources[0] ) EQ 'SZPB' ))[0]
	END
	ELSE: BEGIN					; Safety belt; probably redundant
		pos_fnc = (where( tag_names( point_sources[0] ) EQ 'VV' ))[0]
		pos_sz  = (where( tag_names( point_sources[0] ) EQ 'SZVV' ))[0]
	END
	ENDCASE

	IF IsType(point_size,/defined) THEN BEGIN

		; If point_size[1] is not present or set bad then assume point_onesize=1
		; i.e. plot all sources (including bad ones) with the same circle size.

		point_sz = point_size
		IF n_elements(point_sz) EQ 1 THEN point_sz = [point_sz, BadValue(point_sz)]
		point_onesize = point_onesize OR point_sz[1] EQ BadValue(point_sz)

		point_sz0 = replicate(point_sz[0], pnt_n)
		point_sz  = replicate(point_sz[1], pnt_n)

	ENDIF

	pnt_n_ = pnt_n
	CASE point_onesize OF
	0: good_points = where(finite(vu_point_source(point_sources,field=pos_fnc)), pnt_n)
	1: good_points = indgen(pnt_n)				; Plot all sources
	ENDCASE

	IF pnt_n NE 0 THEN BEGIN

		tmp = point_sources[good_points]

		point_nam = strcompress(vu_point_source(tmp,/name), /remove_all)
		point_pos = vu_point_source(tmp, /location, degrees=degrees)

		point_fnc = vu_point_source(tmp, field=pos_fnc)
		point_tim = vu_point_source(tmp, /time)

		IF IsType(point_sz, /undefined) THEN BEGIN
			point_sz0 = vu_point_source(tmp, /size_base, degrees=degrees)
			IF total(point_sz0) LE 0.0 THEN point_sz0 = replicate(2.5/dpm,pnt_n)
			point_sz  = vu_point_source(tmp, field=pos_sz, degrees=degrees, /size_ratio)
		ENDIF

	ENDIF

ENDIF


IF NOT skymap AND pnt_n EQ 0 THEN message, 'no sky map or point source info specified'

sweep = IsType(geolng, /defined)

InitVar, scale		, 1
InitVar, logscale	, /key
InitVar, equator	, /key
InitVar, mirror 	, 0


; Calculate position of Sun in angular coordinates at time UT

SunEqu = big_eph(UT,	$
	body	= jpl_body(/sun  ,/string)	, $
	center	= jpl_body(/earth,/string)	, $
	/to_equatorial		, $
	/to_sphere			, $
	/precess			, $
	degrees = degrees	, $
	/onebody			, $
	/silent				)
SunEqu = SunEqu[0:1]
SunEcl = CvSky(UT, from_equatorial=SunEqu, /to_ecliptic, degrees=degrees, /silent)

CASE sweep OF

0: BEGIN

	; A 'snapshot' skymap is always centered on the Sun (indicated by zero_point)

	Sun = SunEqu*equator+SunEcl*(1-equator)	; Snapshot

END

1: BEGIN									; Sweep of sky

	; Sweep of the sky centered on the local meridian at time UT, or on the Sun
	; utnoon is nearest noon to UT

	GeographicInfo, UT, location=geolng, utnoon=utnoon, degrees=degrees

	; Why utnoon??

	Sun = big_eph(utnoon,	$
		body	= jpl_body(/sun  ,/string)	, $
		center	= jpl_body(/earth,/string)	, $
		/to_equatorial		, $
		/to_sphere			, $
		/precess			, $
		degrees = degrees	, $
		/onebody			, $
		/silent				)
	Sun = [Sun[0],equator*Sun[1]]

END

ENDCASE

IF IsType(BreakVal, /undefined) THEN BEGIN

	CASE 1 OF
	skymap		: fmin = min(F, /nan, max=fmax)
	pnt_n ne 0	: fmin = min(point_fnc, /nan, max=fmax)
	ENDCASE

	IF logscale 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 logscale THEN BreakVal = 10^BreakVal

ENDIF


CASE mirror OF
0: BEGIN
	xdis =  cosRA
	ydis =  sinRA
	IF IsType(xrange, /defined) THEN xrange_ =  xrange
	IF IsType(yrange, /defined) THEN yrange_ =  yrange
END
1: BEGIN
	xdis = -sinRA
	ydis =  cosRA
	IF IsType(yrange, /defined) THEN xrange_ = -yrange
	IF IsType(xrange, /defined) THEN yrange_ =  xrange
END
2: BEGIN
	xdis = -cosRA
	ydis = -sinRA
	IF IsType(xrange, /defined) THEN xrange_ = -xrange
	IF IsType(yrange, /defined) THEN yrange_ = -yrange
END
3: BEGIN
	xdis =  sinRA
	ydis = -cosRA
	IF IsType(yrange, /defined) THEN xrange_ =  yrange
	IF IsType(xrange, /defined) THEN yrange_ = -xrange
END
4: BEGIN
	xdis =  sinRA
	ydis =  cosRA
	IF IsType(yrange, /defined) THEN xrange_ =  yrange
	IF IsType(xrange, /defined) THEN yrange_ =  xrange
END
5: BEGIN
	xdis = -cosRA
	ydis =  sinRA
	IF IsType(xrange, /defined) THEN xrange_ = -xrange
	IF IsType(yrange, /defined) THEN yrange_ =  yrange
END
6: BEGIN
	xdis = -sinRA
	ydis = -cosRA
	IF IsType(yrange, /defined) THEN xrange_ = -yrange
	IF IsType(xrange, /defined) THEN yrange_ = -xrange
END
7: BEGIN
	xdis =  cosRA
	ydis = -sinRA
	IF IsType(xrange, /defined) THEN xrange_ =  xrange
	IF IsType(yrange, /defined) THEN yrange_ = -yrange
END
ENDCASE

F = rotate(F, mirror)

nxdis = n_elements(xdis)
nydis = n_elements(ydis)

InitVar, xrange_, xdis[[0,nxdis-1]]
InitVar, yrange_, ydis[[0,nydis-1]]

pdis_max = max([xrange_,yrange_])

InitVar, ra_step , 0.5*pi
InitVar, dec_step, pdis_max

pdis = dec_step*(1+indgen(fix(pdis_max/dec_step)))
lras = ra_step*indgen(fix(2*pi/ra_step))

CASE (where(mirror EQ [0,2,5,7]))[0] EQ -1 OF
0: BEGIN
	xras = pdis_max*cos(lras*rpm)
	yras = pdis_max*sin(lras*rpm)
END
1: BEGIN
	xras = pdis_max*sin(lras*rpm)
	yras = pdis_max*cos(lras*rpm)
END
ENDCASE

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

InitVar, noerase, /keyword
IF NOT noerase THEN erase

plot, scale*xrange_, scale*yrange_*(float(!d.y_size)/!d.x_size),	$
	xrange	= scale*xrange_ ,	$
	yrange	= scale*yrange_*(float(!d.y_size)/!d.x_size),	$
	position= position		, $
	xstyle	= 5				, $
	ystyle	= 5 			, $
	/nodata 				, $
	/noerase, clip=0

!p.clip[2] = !p.clip[2]+1				; Fixes bug with IDL clip window

CASE 1 OF
skymap: BEGIN							; ColorSkybox may change skyedge if it's too big
	InitVar, zero_phase, zero_point
	ColorSkyBox, xdis, ydis,	$
		GetColors(F, BreakVal, legend=legend AND NOT naked, logscale=logscale, format=format,	$
		/open, charsize=charsize, _extra=_extra),							$
		skyedge=0, zero_phase=0, dabg=dabg, degrees=degrees,				$
		/black, fill2edge=fill2edge, _extra=_extra
END
pnt_n NE 0:	good = GetColors (point_fnc, BreakVal, legend=legend AND NOT naked, logscale=logscale, $
		format='(F6.2)',/open, charsize=charsize)
ENDCASE

nangle	= 61
angle   = gridgen(nangle, range=[0,2*pi])
xcircle = cos(angle*rpm)
ycircle = sin(angle*rpm)

;IF IsType(minelo, /defined) THEN BEGIN

	; Erase innermost minelo elongations from plot.
	; Set up a cone with opening angle minelo around the Sun in a coordinate system
	; with the z-axis pointing to the Sun, and the x-axis in the ecliptic/equatorial plane.
	; Then convert to ecliptic/equatorial coordinates.

	;pos = transpose([[angle], [replicate(pi2-minelo, nangle)]])
	;pos = EulerRotate([-pi2,-(pi2-Sun[1]),-Sun[0]], pos, degrees=degrees)

	;pos = MercatorProj(pos, zero_phase=zero_point, dabg=dabg, degrees=degrees)

	;polyfill, pos[0,*], pos[1,*], color=!p.background

;ENDIF




IF naked THEN RETURN



;Sun = MercatorProj(Sun, zero_phase=zero_point, dabg=dabg, degrees=degrees)

;rr = 0.015*projsz						; Mark location Sun
;polyfill, Sun[0]+rr*xcircle, Sun[1]+rr*ycircle

; Overplot point sources over skymap

IF pnt_n NE 0 THEN BEGIN
	pos = point_pos						; RA and declination of sources

	IF NOT equator THEN BEGIN

		CASE sweep OF
		0: pos = CvSky(UT, from_equatorial=pos, /to_ecliptic, degrees=degrees, /silent)
		1: BEGIN						; This is the same trick as used in EarthTransit3DLoc
			B   = 23.439291/!radeg		; Angle equator/ecliptic
			AB  = reform(pos[0,*])*rpm
			tmp = sin( atan( sin(AB)/cos(B), cos(AB) ) )
			tmp = atan( tmp*sin(B), tmp*cos(B)/sin(AB) )/rpm
			pos[1,*] = pos[1,*]-tmp		; Subtract decl of ecliptic at RA pos[0,*]
		END
		ENDCASE

	ENDIF

	CASE skymap OF
	0: BEGIN

		; If /point_onesize is set then point_fnc[i] may be bad
		; If NOT set then all point_fnc values will be valid.

		CASE point_onesize OF
		0: dG = point_sz0+point_sz*point_fnc
		1: dG = point_sz0
		ENDCASE

		dG = transpose([[dG],[replicate(0,pnt_n)]])

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

		dG = reform(dG[0,*])

		color = replicate(!p.color, pnt_n)

	END
	1: 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.

		dF = replicate(BadValue(0.0),pnt_n)
		IF a[0] NE -1 THEN dF[a] = F[(i[a]+izero+1) mod nRA,j[a]]
		dF = point_fnc-dF		; Difference with nearest grid point

		; Set up size of circle in degrees. Then convert to data units.
		; Don't merge this to dG = point_sz0+(1-point_onesize)*point_sz*abs(dF)
		; If point_onesize is set we don't want bad values in dG.

		CASE point_onesize OF
		0: dG = point_sz0+point_sz*abs(dF)
		1: dG = point_sz0
		ENDCASE

		dG = transpose([[dG],[replicate(0,pnt_n)]])

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

		dG = reform(dG[0,*])

		; Select colors for point sources from the same color table as the skymap

		good = GetColors(point_fnc, BreakVal, logscale=logscale, /open)

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

		color = (dF GE 0 OR 1-finite(dF))*!p.color+(dF LT 0)*!p.background
	END
	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).

	point_elo = vu_point_source( vu_point_source(point_sources[good_points],UT,/time), /elong, degrees=degrees)

	n = 0
	FOR i=0,pnt_n-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]
			oplot   , pos[0,i]+dG[i]*xcircle, pos[1,i]+dG[i]*ycircle,color=color[i]
			n = n+1
			IF point_names THEN BEGIN
				sgn = pos[0,i] LT 0
				xyouts, pos[0,i]+(1-2*sgn)*dG[i], pos[1,i], point_nam[i], align=-.1*sgn + 1.1*(1-sgn)
			ENDIF
		ENDIF
	ENDFOR

	IF silent LE 0 THEN message, /info, 'plotting'+strcompress(n)+'/'+strcompress(pnt_n_,/rem)+' sources'

	pos = convert_coord(pos, /data, /to_device)
	point_sources[good_points] = vu_point_source(point_sources[good_points], round(pos[0:1,*]), /xy_pixel )
ENDIF

goodcolor = !d.n_colors-1
badcolor  = !p.color

;tmp = max([cosRA,sinRA])
;oplot, tmp*xcircle, tmp*ycircle

; Plot circles around origin at intervals dec_step

FOR k=0,n_elements(pdis)-1 DO BEGIN

	xx = pdis[k]*xcircle
	yy = pdis[k]*ycircle

	IF skymap THEN BEGIN
		i = round( interpol(indgen(nxdis), xdis, xx) )
		j = round( interpol(indgen(nydis), ydis, yy) )

		a = where(0 LE i AND i LE nxdis-1 AND 0 LE j AND j LE nydis-1)

		good = bytarr(nangle)
		IF a[0] NE -1 THEN good[a] = finite(F[i[a],j[a]])
	ENDIF

	CASE goodcolor EQ badcolor OF
	0: PlotCurve, /oplot, xx, yy, good, color=[goodcolor,badcolor], /changecolor, /silent, /noclip
	1: PlotCurve, /oplot, xx, yy, color=goodcolor, /silent, /noclip
	ENDCASE

ENDFOR

; Plot radial vectors from origin to xrange[1] at steps ra_step

dx = 0.01
dy = 0.01

InitVar, ra_position, replicate(0.0,2*n_elements(xras))

one = gridgen(nangle,/one)

FOR k=0,n_elements(xras)-1 DO BEGIN

	xx = xras[k]*one
	yy = yras[k]*one

	IF skymap THEN BEGIN
		i = round( interpol(indgen(nxdis), xdis, xx) )
		j = round( interpol(indgen(nydis), ydis, yy) )

		a = where(0 LE i AND i LE nxdis-1 AND 0 LE j AND j LE nydis-1)

		good = bytarr(nangle)
		IF a[0] NE -1 THEN good[a] = finite(F[i[a],j[a]])
	ENDIF

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

	tmp = [xx[nangle-1],yy[nangle-1]]
	tmp = (convert_coord(tmp,/data,/to_normal))[0:1]

	IF finite(ra_position[2*k+0]) THEN	$
		xyouts, /normal, tmp[0]+dx*ra_position[2*k+0], tmp[1]+dy*ra_position[2*k+1],	$
			strcompress(round(dpm*lras[k]),/rem), charsize=charsize, align=0.5

ENDFOR


F = rotate(F,([0,3,2,1,4,5,6,7])[mirror])

RETURN











FOR plane=0,galactic DO BEGIN

	n = 120								; Plot ecliptic/equator
	pos = transpose( [[gridgen(n, /open, range=pi*[-1,1])],[fltarr(n)]] )

	CASE plane OF
	0: BEGIN
		CASE Equator OF
		0: pos = CvSky(UT,from_equatorial=pos,/to_ecliptic  ,degrees=degrees)
		1: pos = CvSky(UT,from_ecliptic  =pos,/to_equatorial,degrees=degrees)
		ENDCASE
	END
	1: BEGIN
		CASE Equator OF
		0: pos = CvSky(UT,from_galactic=pos,/to_ecliptic  ,degrees=degrees)
		1: pos = CvSky(UT,from_galactic=pos,/to_equatorial,degrees=degrees)
		ENDCASE
	END
	ENDCASE

	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)
		good = bytarr(n)
		IF a[0] NE -1 THEN good[a] = finite(F[(i[a]+izero+1) mod nRA,j[a]])
	ENDIF

	pos[0,*] = pos[0,*]-zero_point
	pos[0,*] = AngleRange(pos[0,*], /pi, degrees=degrees)

	;a	= max(pos[0,*],i)				; Make sure lng decreases
	;pos = shift(pos,0,-i-1)
	;IF skymap THEN good = shift(good,-i-1) 

	a   = sort(pos[0,*])				; This does the same
	pos = pos[*,a]
	IF skymap THEN good = good[a] 

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

ENDFOR



b = 1

pos = projsz*[[-b,0],[b,0]]
oplot, pos[0,*], pos[1,*], color=badcolor	; Draw horizontal axis

a = 1

pos = projsz*[[0,-a],[0,a]]					; Draw vertical axis
oplot, pos[0,*], pos[1,*], color=badcolor

str = (['Ecl ','Eq '])[Equator]				; Label plot with directions

pos = projsz*[[b,0],[-b,0],[0,a],[0,-a]]



xyouts, pos[0,*]+0.0566*projsz*b*[1,-1,0,0], pos[1,*]+a*0.0177*projsz*[1,1,2,-5],	$
	['E','W',str+'N',str+'S'], alignment=0.5, charsize=charsize



InitVar, upto, TimeUnit(/min)

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


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

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

RETURN  &  END