;+
; NAME:
;	vu_coronagraph
; PURPOSE:
;	Wrapper for href=PlotCoronagraph=. Displays a view of the sky from Earth
;	using a 3D heliospheric matrix from one of the tomography programs.
;	This is a special version of the fisheye-maps plotted by
;	vu_earthskymap, optimized for small coronagraph-like patches of sky.
; CATEGORY:
;	sat/idl/util/vupack
; CALLING SEQUENCE:
	PRO vu_coronagraph, hdr, ff,	$
		hdr_bb		= hdr_bb		, $
		bb			= bb			, $
		ut0 		= ut0			, $
		ihdr		= ihdr			, $
		type		= type			, $
		nodata		= nodata		, $
		maxelo		= maxelo		, $
		degrees 	= degrees		, $
		equator 	= equator		, $
		_extra		= _extra		, $
		solar		= solar 		, $
		is_sequence = is_sequence	, $
		point_sources=point_sources , $
		printer 	= printer		, $
		count		= count 		, $
		nEE 		= nEE			, $
		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; default: 7 degrees
;						Determines the size of the fisheye map
;	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
;
;	Three parameters define the grid used to calculate the sky map (see PROCEDURE).
;	These are passed to href=EarthSky3DLoc=:
;
;	nEE = nEE		scalar; type: int; default: 100
;						# points in elongation
;	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=
;					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
;
;	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:
;       Output is displayed on the screen or written to file.
;       See href=vu_get_page= for information on the keywords involved.
; 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
;	EarthSky3DLoc, IPS_params, SuperArray
;	ThomsonUBVConst, IntegrateLOS, AngleRange, PlotCoronagraph, vu_get_page
;	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-2011, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;		Based on vu_earthskymap
;-

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

ihdr_set = IsType(ihdr, /defined)

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

rpm = ToRadians(degrees=degrees)
dpm = ToDegrees(degrees=degrees)
InitVar, maxelo, 7.0/dpm
InitVar, nEE, 100

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

	label = ''

	type = vu_type_skymap(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(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)

	; Need to set up the lines of sight here
	; Set up in pa-angle, elo; then convert to ecliptic long/lat?

	; Set up rectangular grid covering a square of maxelo by maxelo centered on the Sun

	rr = gridgen( 2*nEE*[1,1], range=maxelo*[-1,1]#[1,1], /open )

	; Convert to spherical angles in coordinate frame with z-axis towards sun,
	; and x-axis to ecliptic north
	; (the phase angle is measured clockwise from the north, opposite to the
	; convention for the position angle)

	halfpi = !dpi/rpm/2

	rr = [ atan(rr[0,*],rr[1,*])/rpm, halfpi-reform(sqrt(total(rr*rr,1)),size(rr[0,*],/dim)) ]

	; Convert to sun-centered ecliptic coordinates (Euler angles (90,0,180) degrees

	rr = EulerRotate( halfpi*[0,1,2], rr, degrees=degrees)

	; Convert to ecliptic/equatorial coordinates
	; This is what we need as input for EarthSky3DLoc

	rr = CvSky(utt, from_centered_ecliptic=rr, to_ecliptic=1-equator, to_equatorial=equator, degrees=degrees)

	InitVar, nRR, 20
	InitVar, dRR, 0.1

	; 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[nEE*nEE); type: float]
	;	DEC=DEC	array[nEE*nEE]; type: float
	; 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)

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

	RA  = reform(rr[0,*])
	DEC = reform(rr[1,*])

	R = EarthSky3DLoc(utt		, $
		nRA 		= 0			, $
		nDE 		= 0			, $
		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	)

	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 /= 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


old_device = !D.name

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

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

n = 2*nEE
RA  = reform(RA ,n,n,/overwrite)
DEC = reform(DEC,n,n,/overwrite)
F   = reform(F  ,n,n,/overwrite)

PlotCoronagraph, utt, RA, DEC, F, $
	log 		= log			, $
	equator 	= equator		, $
	breakval	= BreakVal		, $
	charsize	= type.charsize , $
	maxelo		= maxelo		, $
	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]+'_'+'coronagraph'+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