;+
; NAME:
;	PlotCoronagraph
; PURPOSE:
;	Plot sky map, with locations of point sources overplotted (optional)
; CATEGORY:
;	sat/idl/toolbox/graphics
; CALLING SEQUENCE:
	PRO PlotCoronagraph, UT, RA, DEC, F, $
		equator 	= equator		, $
		maxelo		= maxelo		, $
		minelo		= minelo		, $
		drawelo		= drawelo		, $
		degrees 	= degrees		, $
		mirror		= mirror		, $
		noerase 	= noerase		, $
		breakval	= breakval		, $
		basebreak	= basebreak 	, $
		logscale	= logscale		, $
		charsize	= charsize		, $
		title		= title 		, $
		upto		= upto			, $
		scale		= scale 		, $
		zero_point	= zero_point	, $
		zero_phase	= zero_phase	, $
		dabg		= dabg			, $
		format		= format		, $
		fill2edge	= fill2edge 	, $
		_extra		= _extra		, $
		galactic	= galactic		, $
		silent		= silent		, $
		goodcolor	= goodcolor 	, $
		badcolor	= badcolor		, $

		sn_position = sn_position	, $
		we_position = we_position	, $
		ut_position = ut_position	, $
		user_position=user_position , $
		user_align	= user_align	, $
		user_string = user_string	, $

		body		= body			, $

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

		position	= position		, $
		legend		= legend		, $
		compass 	= compass		, $
		max_dec 	= max_dec		, $
		ra_step 	= ra_step		, $
		dec_step	= dec_step		, $
		naked		= naked
; INPUTS:
;	UT			scalar; type: standard time structure
;					time; determines position of Sun in the sky
;	RA			array[n]; type: float
;					array of ecliptic longitudes or right ascensions
;					(if /equator set) in the range [-180,+180]
;					(values are mapped into this range if they are not)
;	DEC			array[m]; type: float
;					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.
;	maxelo		scalar; type: float
;					fish-eye map out to elongation 'maxelo'
;	minelo		scalar; type: float
;					Erases plot inside of minelo
;	drawelo		scalar or array; type: float; default: maxelo
;					elongations of circles to be plotted
;	/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		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.
;	ut_position=ut_position
;				array[2]; type: real
;					x,y coordinates for shifts of UT time 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.
;
;	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'.
;
;	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: centered on Sun
;					(used internally and passed to FishEye,
;					HammerAitoff and MercatorProj). Defines the
;					longitude or RA of the center of the FishEye,
;					HammerAitoff or Mercator map.
;
;	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.
;					zero_phase is not used if F does not exist on
;					input, i.e. if only point sources are plotted
;					(in this case the RA and DEC arrays don't exist
;					either).
;
;	For sky 'snapshots' zero_phase is the same as zero_point.
;
;	Keywords passed to ColorSkybox:
;
;	use_mask = use_mask
;	mask=mask
; OUTPUTS:
; INCLUDE:
	@compile_opt.pro		; On error, return to caller
	@vu_fnc_include.pro 	;
; CALLS:
;	ToRadians, ToDegrees, ColorSkybox, GetColors, big_eph, CvSky
;	FishEye, PlotCurve, AngleRange, gridgen
;	GeographicInfo, TimeGet, TimeUnit, ANiceNum, IsType
;	EulerRotate, InitVar, vu_point_source
;	sphere_distance
; SEE AlSO:
;	ips_point_source__define
; PROCEDURE:
;	Based on PlotEarthSkymap
; MODIFICATION HISTORY:
;	MAR-2011, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;-

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, galactic, /key
InitVar, silent  , 0
InitVar, charsize, 1
InitVar, legend	 , /key
InitVar, naked	 , /key

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 = point_sz[0]
	point_sz  = point_sz[1]

ENDIF

IF pnt_n NE 0 THEN BEGIN

	tmp = vu_los_nodata
	IF IsType(title, /defined) THEN	BEGIN
		tmp = vu_fnc(/los,/all,/name)
		i = min(strlen(tmp))
		tmp = (where(strpos(strmid(tmp,0,i), strmid(title,0,i)) EQ 0))[0]
		IF tmp NE -1 THEN tmp = (vu_fnc(/get_order,/los))[tmp] ELSE tmp = vu_los_nodata
	ENDIF
	IF tmp EQ vu_los_nodata THEN tmp = vu_los_ips_speed

	CASE tmp OF
	vu_los_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
	vu_los_ips_glevel: BEGIN
		pos_fnc = (where( tag_names( point_sources[0] ) EQ 'GG' ))[0]
		pos_sz  = (where( tag_names( point_sources[0] ) EQ 'SZGG' ))[0]
	END
	vu_los_thomson_b: BEGIN
		pos_fnc = (where( tag_names( point_sources[0] ) EQ 'BB' ))[0]
		pos_sz  = (where( tag_names( point_sources[0] ) EQ 'SZBB' ))[0]
	END
	vu_los_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_sz0,/defined) THEN BEGIN
		point_sz0 = replicate(point_sz0, pnt_n)
		point_sz  = replicate(point_sz , 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

InitVar, point_sz0, 1.0/dpm

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

skyedge = maxelo

InitVar, scale		, 1.0
InitVar, logscale	, /key
InitVar, equator	, /key
InitVar, mirror 	, /key

mirror 	= 2*mirror-1


SunEqu = big_eph(UT 	, $ 			; Position of Sun at time UT
	body	= 'Sun'		, $
	center	= 'Earth'	, $
	/to_equatorial		, $
	/to_sphere			, $
	/precess			, $
	degrees = degrees	, $
	/onebody			, $
	/silent				)
SunEqu = SunEqu[0:1]					; Equatorial coordinates
SunEcl = CvSky(UT, from_equatorial=SunEqu, /to_ecliptic, degrees=degrees, /silent)
										; Ecliptic coordinates
; A 'snapshot' skymap is always centered on the Sun
; (indicated by zero_point)

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

InitVar, zero_point, Sun[0]

; Range of Hammer-Aitoff projection is -2*sqrt(2),2*sqrt(2) in horizontal
; and -1*sqrt(2),1*sqrt(2) in vertical direction.

projsz = FishEye( [pi2,0], maxelo=pi2, degrees=degrees )
projsz = projsz[0]

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

xrange = scale*mirror*projsz*[-1,1]
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

plot, xrange, xrange*(float(!d.y_size)/!d.x_size),		$
	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]+1		; Fixes bug with IDL clip window

CASE IsType(breakval, /defined) OF
0: 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

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

nRA = n_elements(RA )
nDE = n_elements(DEC)

CASE 1 OF
skymap: BEGIN		; ColorSkybox may change skyedge if it's too big
	InitVar, zero_phase, zero_point
	ColorSkybox, RA, DEC,	$
		GetColors(F, breakval, legend=legend, logscale=logscale, format=format, /open, charsize=charsize, _extra=_extra),$
		skyedge		= skyedge	, $
		zero_phase	= zero_phase, $
		dabg		= dabg		, $
		degrees		= degrees	, $
		/black					, $
		fill2edge	= fill2edge	, $
		_extra		= _extra
END
pnt_n NE 0:	good = GetColors (point_fnc, breakval, legend=legend, 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)



; minelo can lead to trouble if zero_point
; becomes very big, (for zero_point=180 the polyfill call
; will erase the entire plot!). We avoid this by skipping
; the polyfill call if the sun is outside the fisheye map.

InitVar, minelo, 0.0, set=innerelo

show_sun = 1
origin   = FishEye([0.0d0,0.0d0], /inverse, zero_phase=zero_point, dabg=dabg, degrees=degrees, maxelo=skyedge)
elo_sun  = sphere_distance(sun,origin,degrees=degrees)
show_sun = elo_sun LT maxelo+innerelo

CASE 1 OF

show_sun: BEGIN

	IF innerelo GT 0.0 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 = FishEye(pos, zero_phase=zero_point, dabg=dabg, degrees=degrees, maxelo=skyedge)

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

	ENDIF

	; Mark location of Sun

	Sun = FishEye(Sun, zero_phase=zero_point, dabg=dabg, degrees=degrees, maxelo=skyedge)

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

END

ELSE: BEGIN

	IF IsType(dabg,/defined) THEN BEGIN
		Sun = FishEye(Sun, zero_phase=zero_point, dabg=dabg, degrees=degrees, maxelo=skyedge)
		Sun /= sqrt(total(Sun*Sun)) 	; Unit vector
		Sun *= projsz					; Scale to data units
		oplot, Sun[0]*[0.9,0.96], Sun[1]*[0.9,0.96], thick=2
	ENDIF

END

ENDCASE



dx = 0.01
dy = 0.01

InitVar, upto, TimeUnit(/min)
InitVar, ut_position, [0,0]

xyouts, /normal, $
	0.99+dx*ut_position[0]	, $
	0.95+dy*ut_position[1]	, $
	TimeGet(UT,/ymd,upto=upto,/roundt)+' UT' , $
	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


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

	IF NOT equator THEN 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]

	dG = FishEye( dG, degrees=degrees, maxelo=pi2 )

	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
	pos = FishEye(pos, zero_phase=zero_point, dabg=dabg, degrees=degrees, maxelo=skyedge)

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

IF pnt_n NE 0 THEN BEGIN			; Overplot point sources over skymap

	pos = point_pos					; RA and declination of sources

	IF NOT equator THEN pos = CvSky(UT, from_equatorial=pos, /to_ecliptic, degrees=degrees, /silent)

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

		dG = FishEye( dG, maxelo=pi2, degrees=degrees )

		dG = reform(dG[0,*])

		colors = 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)]])

		dG = FishEye( dG, maxelo=pi2, degrees=degrees )

		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.

		colors = (dF GE 0 OR 1-finite(dF))*!p.color+(dF LT 0)*!p.background

	END
	ENDCASE

	;print, projsz, dg[0], dg[0]/projsz
	pos = FishEye(pos, zero_phase=zero_point, dabg=dabg, degrees=degrees, maxelo=skyedge)

	; 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 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=colors[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=-0.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

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

f = projsz/skyedge

;InitVar, drawelo, maxelo
;IF (where(maxelo EQ drawelo))[0] EQ -1 THEN drawelo = [drawelo,maxelo]

FOR i=0,n_elements(drawelo)-1 DO BEGIN
	a = drawelo[i]
	oplot, f*a*xcircle,f*a*ycircle, color=badcolor, _extra=_extra
	elo_label = strcompress(round(a*dpm),/rem)
	;IF ew THEN BEGIN
		tmp = (convert_coord(f*[ [-a,0],[a,0] ], /data, /to_normal))[0:1,*]
		xyouts, /normal, tmp[0]-0.005, tmp[1]+0.005, elo_label,		$
			charsize=charsize, color=badcolor, align=1, _extra=_extra
		xyouts, /normal, tmp[2]+0.005, tmp[3]+0.005, elo_label,		$
			charsize=charsize, color=badcolor, align=0, _extra=_extra
	;ENDIF ELSE BEGIN
	;	tmp = (convert_coord(f*[ [0,-a],[0,a] ], /data, /to_normal))[0:1,*]
	;	xyouts, /normal, tmp[0]+0.005, tmp[1]+0.005, elo_label,		$
	;		charsize=charsize, color=badcolor, _extra=_extra
	;	xyouts, /normal, tmp[2]+0.005, tmp[3]+0.005, elo_label,		$
	;		charsize=charsize, color=badcolor, _extra=_extra
	;ENDELSE
ENDFOR

tmp = FishEye( [[-pi2,0],[pi2,0],[0,-pi2],[0,pi2]], maxelo=pi2, degrees=degrees )

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 '])[equator]+['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