;+
; NAME:
;	PlotEarthSkymap
; PURPOSE:
;	Plot sky map, with locations of point sources overplotted (optional)
; CATEGORY:
;	sat/idl/toolbox/graphics
; CALLING SEQUENCE:
	PRO PlotEarthSkymap, UT, RA, DEC, F, $
		equator 		= equator		, $
		maxelo			= maxelo		, $
		minelo			= minelo		, $
		geolng			= geolng		, $
		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
;					used to decide on map projection
;					> 0: fish-eye map out to elongation 'maxelo'
;					= 0: Mercator projection (similar to synoptic map)
;					< 0: Hammer-Aitoff projection
;					if absent then a Hammer-Aitoff map is drawn
;	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		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' 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: 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.
;	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.
;
;	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, HammerAitoff, PlotCurve, AngleRange, gridgen
;	GeographicInfo, TimeGet, TimeUnit, ANiceNum, IsType
;	EulerRotate, MercatorProj, InitVar, vu_point_source
;	sphere_distance
; 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
; >	Longitude/RA is plotted increasing right-to-left across the map
;	(as it is for a viewer at Earth looking up at the sky). For
;	sun-centered maps this means that east is left and west is right.
;	A mirror image of the map can be made by setting keyword /mirror.
; > The data range for the horizontal axis is -180,+180 degrees
;	relative to 'zero_point' as needed for the Hammer-Aitoff and
;	fisheye projections. I.e. longitude/RA 'zero_point' will
;	appear at the center of the map.
; > For the fisheye maps the center is at location (longitude,latitude
;	or RA/dec) zero_point+dabg[0],-dabg[1]
;
; 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.
;	SEP-2006, Paul Hick (UCSD/CASS)
;		Modified default for zero_point. Was zero, now is ecliptic
;		longitude/RA of Sun. This automatically centers the Sun
;		if only point sources are plotted
;		(as by href=vu_nagoyasourcemap=).
;	JUL-2008, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;		Modified way the Sun is marked for fish-eye plots.
;		If outside the plot for fishe-eye maps then the Sun is
;		not marked at all, or (if dabg exists) the direction to
;		the Sun is marked with a tickmark on the outer edge.
;-

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

sweep = IsType(geolng, /defined)

CASE IsType(maxelo, /defined) OF
0: BEGIN								; Hammer-Aitoff if maxelo undefined
	hammer	 =  1
	fish	 =  0
	mercator =  0
	skyedge  = -1
END
1: BEGIN								; 'maxelo' is defined
	hammer   = maxelo LT 0.0
	mercator = maxelo EQ 0.0
	fish     = maxelo GT 0.0
	skyedge  = maxelo
END
ENDCASE

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
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	= 'Sun'		, $
		center	= 'Earth'	, $
		/to_equatorial		, $
		/to_sphere			, $
		/precess			, $
		degrees = degrees	, $
		/onebody			, $
		/silent				)
	Sun = [Sun[0],equator*Sun[1]]

END

ENDCASE

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.

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

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)



; For fisheye maps 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
IF fish THEN BEGIN
	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
ENDIF

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)

		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

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

	ENDIF

	; Mark location of Sun

	CASE 1 OF
	fish    : Sun = FishEye		(Sun, zero_phase=zero_point, dabg=dabg, degrees=degrees, maxelo=skyedge)
	hammer  : Sun = HammerAitoff(Sun, zero_phase=zero_point, dabg=dabg, degrees=degrees)
	mercator: Sun = MercatorProj(Sun, zero_phase=zero_point, dabg=dabg, degrees=degrees)
	ENDCASE

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

END

fish: 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 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,*] -= tmp			; Subtract decl of ecliptic at RA pos[0,*]
		END
		ENDCASE

	ENDIF

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

	say, 'plotting position of '+body, threshold=0, silent=silent

	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

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

	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,*] -= 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,*])

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

		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.

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

	END
	ENDCASE

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

	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=colors[i]
			n++
			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

	say, 'plot'+strcompress(n)+'/'+strcompress(pnt_n_,/rem)+' sources', threshold=0, silent=silent

	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

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


	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,/silent)
			1: pos = CvSky(UT,from_ecliptic  =pos,/to_equatorial,degrees=degrees,/silent)
			ENDCASE
		END
		1: BEGIN
			CASE equator OF
			0: pos = CvSky(UT,from_galactic=pos,/to_ecliptic  ,degrees=degrees,/silent)
			1: pos = CvSky(UT,from_galactic=pos,/to_equatorial,degrees=degrees,/silent)
			ENDCASE
		END
		ENDCASE

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

		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] 

		IF hammer THEN pos = HammerAitoff(pos, dabg=dabg, 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

	ENDFOR

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