;+
; NAME:
;	vu_earthskymap
; PURPOSE:
;	Wrapper for href=PlotEarthSkymap=. Displays a view of the sky from Earth
;	using a 3D heliospheric matrix from one of the tomography programs.
; CATEGORY:
;	WWW: graphics
; CALLING SEQUENCE:
	PRO vu_earthskymap, hdr, ff,	$
		hdr_bb		= hdr_bb		, $
		bb			= bb			, $
		ut0 		= ut0			, $
		ihdr		= ihdr			, $
		type		= type			, $
		nodata		= nodata		, $
		maxelo		= maxelo		, $
		geolng		= geolng		, $
		degrees 	= degrees		, $
		equator 	= equator		, $
		_extra		= _extra		, $
		solar		= solar 		, $
		band		= band			, $
		F_full		= FF_full		, $
		RA_full 	= RA_full		, $
		is_sequence = is_sequence	, $
		point_sources=point_sources , $
		hardcopy	= hardcopy		, $
		count		= count 		, $
		nRA 		= nRA			, $
		nDE 		= nDE			, $
		nRR 		= nRR			, $
		dRR 		= dRR			, $
		silent		= silent		, $
		module		= module		, $
		zero_shift	= zero_shift	, $
		wavelength	= wavelength	, $
		filter_bb	= filter_bb 	, $
		center		= center	, $
		use_mask	= use_mask
; INPUTS:
;	hdr				array[1]; type: structure
;						information about the tomographic reconstruction.
;						if 'hdr' does not exist, then user is prompted for
;						a data file using the 'pickfile' dialog.
;	ff				array[n,m,l,2]; type: float
;						3D volume data; if not present then the files
;						specified in the 'hdr' argument are read.
; OPTIONAL INPUT PARAMETERS:
;	ut0=ut0			scalar or array[1]; type: time structure, float or integer
;						Specifies the time on which the time series is centered.
;						If not specified then the time from the 'hdr' argument is used.
;						See href=Carrington= for more details.
;	maxelo			scalar; type: float
;						if set and non-zero then a fish-eye map is plotted out to
;						elongation 'maxelo'. By default a Hammer-Aitoff map is drawn
;	destination=destination
;					if set to valid directory then graphics output is saved as GIF files
;	/degrees		if set, indicates that all angles are in degrees. Default: radians.
;	/equator		if set, the 'horizontal' in the sky map is parallel to the
;						equatorial plane. By default the horizontal is along the ecliptic
;	geolng			scalar; type: float
;						geographic longitude of location on Earth (in units implied
;						by /degrees keyword). If 'geolng' is specified then a 'transit map'
;						is calculated for a one-day period centered on 'ut'
;
;	Four parameters define the grid used to calculate the sky map (see PROCEDURE).
;	These are passed to href=EarthSky3DLoc= and href=EarthTransit3DLoc=:
;
;	nRA = nRA		scalar; type: int; default: 72
;						# points in skymap in longitudinal direction
;						(either right ascension or ecliptic longitude)
;	nDE	= nDE		scalar; type: int; default: 36
;						# points in skymap in latitudinal direction
;						(either declination or ecliptic latitude)
;	nRR	= nRR		scalar; type: int; default: 20
;						# steps along lines of sight at each point in the sky
;	dRR	= dRR		scalar; type: float; default: 0.1 AU
;						step size along line of sight (AU)
;
;	zero_shift		passed to href=EarthSky3DLoc=
;					(NOT used for sky sweeps, i.e. NOT used if 'geolng' is set).
;					By default the skymap is centered on the Sun. Setting 'zero_shift' moves
;					the center to increasing RA (or ecliptic longitude) by 'zero_shift'.
;					E.g. to produces an anti-solar skymap set zero_shift=!pi
;					(or zero_shift=180, /degrees).
;
;	point_sources=point_sources
;					array[n]; type: structure
;						information about IPS observation (location in sky, IPS speed
;						and g-level; see href=getnagoyasources=.
;	/names			if set, then the names of the point sources are plotted
;
;	/solar
;
;	(only if keyword 'band' is used)
;
;	band=band		scalar; type float
;						used to replace a limited strip in RA in an
;						already existing skymap.
;						band: width of strip in hours
;						The strip is inserted into arrays specified
;						in keywords F_full and RA_full.
;	RA_full=RA_full	array(nRA); type: float
;	 F_full=FF_full	array(nRA,nDE); type: float
;						skysweep arrays (RA and function value) where
;						only a strip of width band needs to be updated.
;						Usually these are obtained from a previous call
;						to this procedure. If they don't exist yet then
;						keyword 'band' is ignored and a full skysweep
;						map is intialized.
;
;	Several keywords are passed to href=vu_type_skymap=:
;	The most important ones set the data type to be plotted:
;		/ips_speed			IPS speed
;		/ips_glevel 		IPS g-level
;		/brightness 		Thomson scattering brightness
;		/pb 				Thomson scattering polarization brightness
;		/rotationmeasure	Faraday rotation measure
;
;	type=type		array[1]; type: structure
;
;	wavelength=wavelength
;					scalar; type: any; default: none
;						radio wavelength in meters (used only for Faraday
;						rotation maps to convert the rotation measure
;						(in degrees/m^2) to rotation (degrees).
; OUTPUTS:
;	(none)
; OPTIONAL OUTPUTS:
;	(only if keyword 'band' is used)
;
;	RA_full=RA_full	array(nRA); type: float: default: none
;	 F_full=FF_full	array(nRA,nDE); type: float; default: none
;						updated skysweep arrays (RA and function value)
;						with new strip of width band inserted
; INCLUDE:
	@compile_opt.pro			; On error, return to caller
	@vu_fnc_include.pro 		;
; EXTERNAL:
;	IPS_velocity, IPS_velocity_w, IPS_hlevel_w, ThomsonBrightness
; CALLS:
;	Carrington, IsType, vu_getdata, vu_gettime, vu_type_skymap, set_page
;	EarthTransit3DLoc, EarthSky3DLoc, IPS_params, SuperArray
;	ThomsonUBVConst, IntegrateLOS, AngleRange, PlotEarthSkymap, vu_get_page
;	reset_colors, InitVar, vu_get, vu_point_source, ThomsonRadialFilter
;	TimeGet, TimeUnit
; PROCEDURE:
;	The sky map is build from nRA*nDE lines of sight integrations
;	through the F3D array.
;	The locations in the sky of the lines of sight are
;		RA/Lng    = (RA/Lng Sun) -180 + (360./nRA)*[0.5, 1.5, .. , nRA-0.5 ]
;		Decl./Lat = - 90 + (180./nDE)*[0.5, 1.5, .. , nDE-0.5 ]
;		Dist along los = dRR*[0.5,1.5,...,nRR-1]
;	Note that each line of sight is centered on a 'box' of 360/nRA by
;	180/nDE degrees, and that all lines of sight together cover the
;	entire sky. The edges of the boxes for all lines of sight are
;	returned in the RA and DE keywords
; MODIFICATION HISTORY:
;	MAR-2000, Paul Hick (UCSD/CASS)
;	AUG-2000, Paul Hick (UCSD/CASS)
;		Added option to accept a vu_read structure as a 1st argument
;	SEP-2001, Paul Hick (UCSD/CASS)
;		Fixed bug in vu_type_skymap calls
;	JAN-2002, Paul Hick (UCSD/CASS)
;		Changed the initialization procedure for a skysweep with the 'band'
;		keyword set. The initialization now ignores the 'band' value and simply
;		calculates a full skysweep map (the earlier initialization put a new
;		strip in an otherwise empty sky).
;	JUL-2002, Paul Hick (UCSD/CASS)
;		Added interpolation between neigbouring times of ut0, instead of
;		truncating to nearest grid time.
;	AUG-2004, Paul Hick (UCSD/CASS)
;		Modifed IntegrateLOS calls. For Thomson scattering brightness
;		maps with point sources overplotted a radial filter was applied
;		to the maps but not to the point sources. Fixed.
;	OCT-2004, Paul Hick (UCSD/CASS)
;		Added keyword zero_shift.
;	JUN-2007, Paul Hick (UCSD/CASS)
;		Added keyword wavelength
;	SEP-2007, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;		Bugfix. Added filter_bb as keyword.
;-

InitVar, silent  , 0
InitVar, nodata  , /key
InitVar, equator , /key
InitVar, hardcopy, /key

ihdr_set = IsType(ihdr, /defined)

; Set maxelo to non-zero scalar value to get fisheye map (e.g. !pi/2*1.25)

hammer = IsType(maxelo, /undefined)			; If maxelo undefined then plot Hammer-Aitoff map
IF NOT hammer THEN hammer = maxelo LE 0		; If maxelo exists plot fisheye map only if maxelo>0

sweep = IsType(geolng, /defined)

CASE nodata OF
1: BEGIN
	utt = Carrington(ut0,/get_time)

	merge_band = 0
	label = ''

	type = vu_type_skymap(hammer=hammer, type=type,_extra=_extra)
	id_type = vu_los_nodata

END
0: BEGIN

	; Select and read the tomography files

	utt = vu_getdata(hdr, ff, ut0=ut0, ihdr=ihdr, silent=silent, count=count, _extra=_extra)
	IF count EQ 0 THEN RETURN				; Return if no data files were read

	utt = Carrington(utt,/get_time)

	type = vu_type_skymap(hammer=hammer, type=type, _extra=_extra)
	id_type = (where(type.data))[0]

	; Get magnetic field data if necessary

	NeedBB = id_type EQ vu_los_faraday_rm
	HaveBB = NeedBB AND IsType(hdr_bb,/defined) AND IsType(bb,/defined)

	IF NeedBB AND NOT HaveBB THEN BEGIN

		SetFileSpec, vu_get(hdr, /file)
		InitVar, filter_bb, vu_prefix(/magnetic,/silent)
		nearest = 0

		CASE nearest OF
		0: file = GetFileSpec(upto='directory')+filter_bb+strmid(GetFileSpec(from='name'),4)
		1: BEGIN
			file = (GetFileSpec(upto='directory'))[0]+filter_bb+'*.*'
			file = file_search(file, count=count_bb)
			NeedBB = count_bb GT 0
		END
		ENDCASE

		IF NeedBB THEN BEGIN
			hdr_bb = vu_select(file, /check, /read, count=count_bb, ff=bb)
			HaveBB = count_bb EQ count
		ENDIF

		IF NOT HaveBB THEN BEGIN
			message, /info, 'no magnetic field files found'
			RETURN
		ENDIF

	ENDIF

	; Pick out the matrix at time utt or hdr[ihdr].time by interpolation

	CASE ihdr_set OF
	0: BEGIN
		hdri = vu_gettime(hdr, ff, ut=utt, ff_ut=ffi)
		IF HaveBB THEN hdri_bb = vu_gettime(hdr_bb, bb, ut=utt, ff_ut=bbi)
	END
	1: BEGIN
		hdri = vu_gettime(hdr, ff, ut=vu_get(hdr[ihdr], /uttime), ff_ut=ffi, /fixgrid); ihdr set explicitly
		IF HaveBB THEN hdri_bb = vu_gettime(hdr_bb, bb, ut=vu_get(hdr[ihdr], /uttime), ff_ut=bbi, /fixgrid)
	END
	ENDCASE


	xcrange = vu_get(hdri, /array_range)

	nLng	= vu_get(hdri, /nlng)
	nLat	= vu_get(hdri, /nlat)
	nRad	= vu_get(hdri, /nrad)

	R3D		= vu_get(hdri, /start_distance)
	dR3D	= vu_get(hdri, /distance_step)

	betan	= (vu_get(hdri, /power_density ))[0]
	betar	= (vu_get(hdri, /power_distance))[0]

	n0		= vu_get(hdri, /density_at_1AU)

	InitVar, nRA, 72
	InitVar, nDE, 36
	InitVar, nRR, 20
	InitVar, dRR, 0.1

	; For a skysweep the keyword 'band' indicates that only a strip of sky needs to be inserted
	; into an existing skysweep map.
	; If RA_full (and hence FF_full) doesn't exist yet, ignore the 'band' setting and initialize
	; a full skysweep map (the band0 variable passed to EarthTransit3DLoc won't exist).

	merge_band = sweep AND IsType(RA_full, /defined) AND IsType(band, /defined)
	IF merge_band THEN band0 = band

	; Calculate locations (in heliographic coordinates for all points along all
	; lines of sight. Note that R[0,*,*,*] contains heliographic longitudes in [0,360]
	;	RA =RA	array[nRA]; type: float
	;		RA  = -pi+2*pi/nRA*[0.5,1.5,...,nRA-0.5]
	;	DEC=DEC	array[nDE]; type: float
	;		DEC = -pi/2+pi/nDE*[0.5,1.5,...,nDE-0.5]
	; RA and DEC define a regular grid in sky coordinates for the centers of
	; the boxes for the lines of sight used to build the skymap
	; (they are output from EarthSky3DLoc or EarthTransit3DLoc)

	; The return array R has longitude angles in R[0,*] (not Carrington variables)

	CASE sweep OF

	0: R = EarthSky3DLoc(utt	, $
		nRA 		= nRA		, $
		nDE 		= nDE		, $
		nRR 		= nRR		, $
		dRR 		= dRR		, $
		RA			= RA		, $
		DEC 		= DEC		, $
		rr_earth	= rr_earth	, $
		ut_earth	= ut_earth	, $
		pa_earth	= pa_earth	, $
		elo_earth	= elo_earth , $
		elo_sun 	= elo_sun	, $
		equator 	= equator	, $
		degrees 	= degrees	, $
		/to_heliographic		, $
		zero_point	= zero_point, $
		zero_phase	= zero_phase, $
		zero_shift	= zero_shift, $
		center		= center	)

	1: R = EarthTransit3DLoc(utt, $
		nRA 		= nRA		, $
		nDE 		= nDE		, $
		nRR 		= nRR		, $
		dRR 		= dRR		, $
		RA			= RA		, $
		DEC 		= DEC		, $
		rr_earth	= rr_earth	, $
		ut_earth	= ut_earth	, $
		pa_earth	= pa_earth	, $
		elo_earth	= elo_earth , $
		elo_sun 	= elo_sun	, $
		equator 	= equator	, $
		degrees 	= degrees	, $
		/to_heliographic		, $
		zero_point	= zero_point, $
		zero_phase	= zero_phase, $
		geolng		= geolng	, $
		solar		= solar 	, $
		band		= band0 	)

	ENDCASE

	label = type.label

END
ENDCASE

; This gets the SMEI mask for the specified day.
; Only makes sense for processsing SMEI tomography files

IF IsType(use_mask,/defined) THEN BEGIN
	mask = smei_sky_getmask(utt,ra,dec,equator=equator,degrees=degrees)
	IF IsType(mask,/scalar) THEN destroyvar, use_mask, mask
ENDIF

nagoya  = 1
fillbad = 0
minRR	= 10

BreakVal = type.breakval

i = where(finite(BreakVal), nbreak)
IF nbreak GT 0 THEN BEGIN
	BreakVal = type.breakval[i]
	CASE type.format of
	''  : BreakVal = round(BreakVal)
	ELSE: format = type.format
	ENDCASE
ENDIF

ctt = TimeGet(utt,/ymd,upto=TimeUnit(/minute))

; Since R[0,*] are longitude angles, and xcrange is given in
; Carrington variables, IntegrateLOS needs to convert the
; longitude angles to Carrington variables.
; Hence, /cv2carrington is set in all IntegrateLOS calls.

CASE id_type OF

vu_los_nodata: IF silent LE 1 THEN message, /info, 'no sky data'

vu_los_ips_speed: BEGIN					; IPS speed
	IF silent LE 1 THEN message, /info, 'IPS velocity at '+ctt

	D3D = ffi[*,*,*,1]					; Normalized density

	id = where( finite(D3D) )
	IF id[0] NE -1 THEN D3D[id] = (D3D[id]/n0)^(2*betan)

	ips_info = [IPS_params(nagoya=nagoya, cambridg=cambridge),betar]

	F = IntegrateLOS(R, ffi[*,*,*,0], D3D,$
		/cv2carrington				 	, $
		R3D 		= R3D				, $
		dR3D		= dR3D				, $
		xcrange 	= xcrange			, $
		minRR		= minRR 			, $
		fillbad 	= fillbad			, $
		degrees 	= degrees			, $
		f3dfnc		= 'IPS_velocity'	, $ 
		f3darg		= ips_info			, $
		w3dfnc		= 'IPS_velocity_w'	, $ 
		w3darg		= ips_info			, $
		rr_earth	= rr_earth			, $
		ut_earth	= ut_earth			, $
		pa_earth	= pa_earth			, $
		elo_earth	= elo_earth 		, $
		elo_sun 	= elo_sun			, $
		silent		= silent			)
END

vu_los_ips_glevel: BEGIN				; IPS g-level
	IF silent LE 1 THEN message, /info, 'IPS g-level at '+ctt

	D3D = ffi[*,*,*,1]					; Normalized density

	id = where( finite(D3D) )
	IF id[0] NE -1 THEN D3D[id] = (D3D[id]/n0)^(2*betan)

	ips_info = [IPS_params(nagoya=nagoya, cambridg=cambridge),betar]

	F = IntegrateLOS(R, D3D 			, $
		/cv2carrington					, $
		R3D 		= R3D				, $
		dR3D		= dR3D				, $
		xcrange 	= xcrange			, $
		minRR		= minRR 			, $
		fillbad 	= fillbad			, $
		degrees 	= degrees			, $
		w3dfnc		= 'IPS_hlevel_w'	, $
		w3darg		= ips_info			, $
		rr_earth	= rr_earth			, $
		ut_earth	= ut_earth			, $
		pa_earth	= pa_earth			, $
		elo_earth	= elo_earth 		, $
		elo_sun 	= elo_sun			, $
		silent		= silent			)

	F = sqrt(F)
END

vu_los_thomson_b: BEGIN					; Thomson scattering brightness
	IF silent LE 1 then message, /info, 'Thomson scattering brightness at '+ctt

	f3darg = [dRR,ThomsonUBVConst(3),0]

	F = IntegrateLOS(R, ffi[*,*,*,1]	, $
		/cv2carrington				 	, $
		R3D 		= R3D				, $
		dR3D		= dR3D				, $
		xcrange 	= xcrange			, $
		minRR		= minRR 			, $
		fillbad 	= fillbad			, $
		degrees 	= degrees			, $
		f3dfnc		= 'ThomsonBrightness',$
		f3darg		= f3darg			, $
		rr_earth	= rr_earth			, $
		ut_earth	= ut_earth			, $
		pa_earth	= pa_earth			, $
		elo_earth	= elo_earth 		, $
		elo_sun 	= elo_sun			, $
		silent		= silent			)

	F /= ThomsonRadialFilter(elo_earth[*,*,0], degrees=degrees)

	IF IsType(point_sources,/defined) THEN BEGIN

		; Apply same radial filter for point sources and skymap.
		; Note that we modify the brightness in the 'point_sources' structures.
		; The input values are restored before returning.

		point_brightness = vu_point_source(point_sources, /brightness)
		tmp = ThomsonRadialFilter(vu_point_source(point_sources, /elo, degrees=degrees), degrees=degrees)
		point_sources = vu_point_source(point_sources, point_brightness/tmp, /brightness)
	ENDIF

	;log = 1
END

vu_los_thomson_pb: BEGIN				; Polarization brightness
	IF silent LE 1 THEN message, /info, 'Thomson scattering pB at '+ctt

	Fb = replicate(n0, nLng,nLat,nRad)	; Normalized density in 1/r^2 density

	f3darg = [dRR,ThomsonUBVConst(3),1]

	Itb = IntegrateLOS(R, Fb			, $
		/cv2carrington				 	, $
		R3D 		= R3D				, $
		dR3D		= dR3D				, $
		xcrange 	= xcrange			, $
		minRR		= minRR 			, $
		fillbad 	= fillbad			, $
		degrees 	= degrees			, $
		f3dfnc		= 'ThomsonBrightness',$
		f3darg		= f3darg			, $
		rr_earth	= rr_earth			, $
		ut_earth	= ut_earth			, $
		pa_earth	= pa_earth			, $
		elo_earth	= elo_earth 		, $
		elo_sun 	= elo_sun			, $
		silent		= silent			)

	It  = IntegrateLOS(R, ffi[*,*,*,1]	, $
		/cv2carrington				 	, $
		R3D 		= R3D				, $
		dR3D		= dR3D				, $
		xcrange 	= xcrange			, $
		minRR		= minRR 			, $
		fillbad 	= fillbad			, $
		degrees 	= degrees			, $
		f3dfnc		= 'ThomsonBrightness',$
		f3darg		= f3darg			, $
		rr_earth	= rr_earth			, $
		ut_earth	= ut_earth			, $
		pa_earth	= pa_earth			, $
		elo_earth	= elo_earth 		, $
		elo_sun 	= elo_sun			, $
		silent		= silent			)

	f3darg = [dRR,ThomsonUBVConst(3),2]

	Irb = IntegrateLOS(R, Fb			, $
		/cv2carrington				 	, $
		R3D 		= R3D				, $
		dR3D		= dR3D				, $
		xcrange 	= xcrange			, $
		minRR		= minRR 			, $
		fillbad 	= fillbad			, $
		degrees 	= degrees			, $
		f3dfnc		= 'ThomsonBrightness',$
		f3darg		= f3darg			, $
		rr_earth	= rr_earth			, $
		ut_earth	= ut_earth			, $
		pa_earth	= pa_earth			, $
		elo_earth	= elo_earth 		, $
		elo_sun 	= elo_sun			, $
		silent		= silent			)

	Ir  = IntegrateLOS(R, ffi[*,*,*,1]	, $
		/cv2carrington				 	, $
		R3D 		= R3D				, $
		dR3D		= dR3D				, $
		xcrange 	= xcrange			, $
		minRR		= minRR 			, $
		fillbad 	= fillbad			, $
		degrees 	= degrees			, $
		f3dfnc		= 'ThomsonBrightness',$
		f3darg		= f3darg			, $
		rr_earth	= rr_earth			, $
		ut_earth	= ut_earth			, $
		pa_earth	= pa_earth			, $
		elo_earth	= elo_earth 		, $
		elo_sun 	= elo_sun			, $
		silent		= silent			)

	F = (It-Ir)/(Itb-Irb)
	;F = ((It-Ir)/(It+Ir))/((Itb-Irb)/(Itb+Irb))

	;log = 1
END

vu_los_column_density: BEGIN			; Column density
	IF silent LE 1 then message, /info, 'Column density at '+ctt

	f3darg = [dRR,ThomsonUBVConst(3),0]

	F = IntegrateLOS(R, ffi[*,*,*,1]	, $
		/cv2carrington				 	, $
		R3D 		= R3D				, $
		dR3D		= dR3D				, $
		xcrange 	= xcrange			, $
		minRR		= minRR 			, $
		fillbad 	= fillbad			, $
		degrees 	= degrees			, $
		f3dfnc		= 'ThomsonColumnDensity',$
		f3darg		= f3darg			, $
		rr_earth	= rr_earth			, $
		ut_earth	= ut_earth			, $
		pa_earth	= pa_earth			, $
		elo_earth	= elo_earth 		, $
		elo_sun 	= elo_sun			, $
		silent		= silent			)

	F = F/ThomsonRadialFilter(elo_earth[*,*,0], degrees=degrees)

	IF IsType(point_sources,/defined) THEN BEGIN

		; Apply same radial filter for point sources and skymap.
		; Note that we modify the brightness in the 'point_sources' structures.
		; The input values are restored before returning.

		point_brightness = vu_point_source(point_sources, /brightness)
		tmp = ThomsonRadialFilter(vu_point_source(point_sources, /elo, degrees=degrees), degrees=degrees)
		point_sources = vu_point_source(point_sources, point_brightness/tmp, /brightness)
	ENDIF

	;log = 1
END

vu_los_faraday_rm: BEGIN

	IF silent LE 1 then message, /info, 'Faraday rotation measure at '+ctt

	; Multiply vector magnetic field with density
	; The product n*B is passed to IntegrateLos for integration along line of sight

	D3D = bb*SuperArray(ffi[*,*,*,1],3,/trail)

	f3darg = [dRR]
	IF IsType(wavelength,/defined) THEN f3darg = [f3darg,wavelength]

	F = IntegrateLOS(R, D3D 			, $
		/cv2carrington				 	, $
		R3D 		= R3D				, $
		dR3D		= dR3D				, $
		xcrange 	= xcrange			, $
		minRR		= minRR 			, $
		fillbad 	= fillbad			, $
		degrees 	= degrees			, $
		f3dfnc		= 'RotationMeasure' , $
		f3darg		= f3darg			, $
		rr_earth	= rr_earth			, $
		ut_earth	= ut_earth			, $
		pa_earth	= pa_earth			, $
		elo_earth	= elo_earth 		, $
		elo_sun 	= elo_sun			, $
		silent		= silent			)

	;dpm = ToDegrees(degrees=degrees)
	;F *= (elo_earth[*,*,0]*dpm/45.0)^3
END

ENDCASE


IF merge_band THEN BEGIN			; Insert strip into skysweep map

	; Keyword 'band' is set. Only a strip of band hours in RA is updated.
	; RA (RA in band) is set in EarthTransit3DLoc using utt and band.

	nRA = n_elements(RA)

	; Insert the strip of band hours. For each of the nRA meridian scans find the location
	; in RA nearest to it. Insert at that position.

	np = intarr(nRA, /nozero)
	FOR i=0,nRA-1 DO BEGIN
		tmp = min(abs(RA[i]-RA_full), n)
		np[i] = n
	ENDFOR
	FF_full[np,*] = F
	RA_full[np  ] = RA

	F  = FF_full
	RA = RA_full

ENDIF


old_device = !D.name

set_page, $
	zbuffer 	= 1-hardcopy	, $
	printer 	= hardcopy		, $
	xysize		= type.xysize	, $
	fullsize	= hardcopy		, $
	_extra		= _extra		, $
	old_device	= old_device	, $
	ctable		= type.ctable	, $
	/silent

reset_colors, /reverse
IF NOT hardcopy THEN stretch, !d.n_colors-1, 0

IF id_type NE -1 THEN BEGIN

	title = vu_fnc(id_type,/name,/los)
	IF vu_fnc(id_type,/unit,/los) NE '' THEN title += ' ('+vu_fnc(id_type,/unit,/los)+')'

	; If the wavelength is specified for the rotation measure
	; then change the label here from "Rotation Measure (degrees/m^2)"
	; to "Faraday rotation (degrees)"
	; (the multiplication by lambda^2 is done in RotationMeasure
	; as called in IntegrateLos).

	IF id_type EQ vu_los_faraday_rm THEN	$
		IF IsType(wavelength,/defined) THEN $
			title = 'Faraday rotation '+strmid(title,strpos(title,'('),strpos(title,'/')-strpos(title,('(')))+')'

ENDIF

PlotEarthSkymap, utt, RA, DEC, F, $
	log 		= log			, $
	equator 	= equator		, $
	breakval	= BreakVal		, $
	charsize	= type.charsize , $
	maxelo		= maxelo		, $
	geolng		= geolng		, $
	degrees 	= degrees		, $
	title		= title 		, $
	zero_point	= zero_point	, $
	zero_phase	= zero_phase	, $
	format		= format		, $
	point_sources=point_sources , $
	silent		= silent		, $
	use_mask	= use_mask		, $
	mask		= mask			, $
	_extra		= _extra

CASE nodata OF
0: hdri = hdr[ihdr]
1: hdri = filepath(root=getenv('TUB'),'source_map.tmp')
ENDCASE

module = type.display+'_'+(['ecl','equ'])[equator]+'_'+(['fisheye','hammer'])[hammer]+label
vu_get_page, hdri, title='', module, ut=utt, device=old_device, silent=silent, _extra=_extra

IF IsType(point_brightness,/defined) THEN	$
	point_sources = vu_point_source(point_sources, point_brightness, /brightness)

RETURN  &  END