;+
; NAME:
;	vu_elotime
; PURPOSE:
;	Produces a "J-map" (time vs elongation for a line-of-sight
;	integrated quantity from a time-dependent tomography solution
; CATEGORY:
;	sat/idl/util/vupack
; CALLING SEQUENCE:
	PRO vu_elotime, hdr, ff, ut0=ut0, ihdr=ihdr,	$

		ts_all 		= ts_all	, $		; Defines times in time series
		ts_range	= ts_range	, $
	 	ts_step		= ts_step	, $

		pa_angle	= pa_angle	, $
		elo_grid	= elo_grid	, $

		nRR 		= nRR		, $
		dRR 		= dRR		, $

		degrees 	= degrees	, $

		silent		= silent	, $
		_extra		= _extra	, $

		delta_time	= delta_time, $
		readall 	= readall	, $
		printer		= printer
; INPUTS:
;	hdr				array[k]; type: structure
;						information about the tomographic reconstruction
;						(single file or time-dependent sequence)
;						if 'hdr' does not exist, then user is prompted for
;						a data file using the 'pickfile' dialog.
;	ff				array[n,m,l,2,k]; type: float
;						3D volume data associated with 'hdr'
;						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 'hdr[ihdr]' is used.
;						See href=Carrington= for more details.
;	ihdr=ihdr		scalar; type: integer; default: 0
;						index of header to be used for display
;
;	Keywords to define the time series (i.e. the times at which data are evaluated).
;	Either ts_all is used, or the pair ts_range and ts_step. If neither of these is set then
;	the times in 'hdr' are used).
;
;	ts_all=ts_all	array; type: float or time structures
;						times at which data are to be evaluated, as Carrington times
;						or UT times.
;	ts_step=ts_step	scalar; type: float or difference time structure: default: 0.25
;						time step in time series in days.
;	ts_range=ts_range
;					array[2]; type: float; default: [-15,15]
;						range of times centered on ut0 in days. The time series covers
;						[ut0-ts_range[0],ut0+ts_range[1]] in steps of ts_step.
;
; INCLUDE:
	@compile_opt.pro
	@vu_fnc_include.pro 		;
; MODIFICATION HISTORY:
;	DEC-2008, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;-

InitVar, silent 	, 0
InitVar, degrees	, /key
InitVar, readall	, /key
InitVar, printer	, /key

uday = TimeUnit(/day )
uhr  = TimeUnit(/hour)
umin = TimeUnit(/min )

running_diff = IsType(delta_time,/defined)
IF running_diff THEN				$
	IF NOT IsTime(delta_time) THEN	$
		delta_time = TimeSet(/diff,delta_time,uhr)

dpm = ToDegrees(degrees=degrees)
rpm = ToRadians(degrees=degrees)

InitVar, pa_angle, 80.0d0/dpm
InitVar, elo_grid, gridgen(111,range=[10,120])/dpm

; If /readall is set then read all the files and store the arrays in ff
; If /readall NOT set then only read all the headers.

CASE readall OF
0: BEGIN
	CASE IsType(hdr,/defined) OF
	0: hdr = vu_select(/pick, /check, silent=silent, _extra=_extra, count=count)
	1: count = n_elements(hdr)
	ENDCASE
	IF count NE 0 THEN utt = vu_get(hdr[0],/uttime)
END
1: BEGIN
	utt = vu_getdata(hdr, ff, ut0=ut0, ihdr=ihdr, silent=silent, _extra=_extra, count=count)
	IF count NE 0 THEN utt = Carrington(utt,/get_time)
END
ENDCASE

IF count EQ 0 THEN RETURN

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

InitVar, nosequence, /key
nosequence OR= n_elements(hdr) EQ 1
IF NOT nosequence THEN nosequence = 1-vu_is_sequence(hdr,silent=silent)

IF IsType(ts_range, /undefined) AND IsType(ts_step, /undefined) THEN BEGIN

	; It neither ts_range nor ts_step are set, use the times from the headers
	; (unless ts_all is input)

	IF nosequence EQ 0 THEN InitVar, ts_all, TimeGet(vu_get(hdr,/uttime),uhr,/round)

ENDIF


nts_all = n_elements(ts_all)

IF nts_all EQ 0 THEN BEGIN

	InitVar, ts_range, [-15.0,15.0]		; Range +/- 15 days
	InitVar, ts_step , 0.25

	IF IsTime(ts_range) THEN ts_range = TimeSet(day=ts_range)
	IF IsTime(ts_step ) THEN ts_step  = TimeSet(day=ts_step )

	nts_all = round( (ts_range[1]-ts_range[0])/float(ts_step)+1 )

	; The next roundoff is not really necessary, but makes the timeseries in
	; movies look less jumpy, because identical times are used.

	;ts_all = TimeGet(utt, /roundt, day=ts_step)
	ts_all = utt
	ts_all = TimeOp(/add, ts_all, TimeSet(day=gridgen(nts_all, range=ts_range)) )
										; utt +/- 15 days with time every 6 hours
ENDIF

nts_all = n_elements(ts_all)

InitVar, power, vu_get(hdr,/power_norm); Needed by vu_timeseries (default is zero)

; rr defines a radial direction away from the Sun at position
; angle pa_angle at elongations elo_grid.

rr = SuperArray(elo_grid,2,/lead)

; Make rr[0,*],rr[1,*] "longitude" and "latitude" in
; coordinate system with z-axis in the anti-solar direction,
; and x-axis to ecliptic north. This makes the "longitude"
; angle equal to the position angle (increasing
; counter-clockwise from ecliptic north, and 90+latitude
; equal to the solar elongation angle
; (i.e. latitude = elongation-90).

rr[0,*]  = pa_angle				; "phase angle"
rr[1,*] -= 90.0d0/dpm			; "latitude angle"

; Change from helio-ecliptic (sun-centered ecliptic)
; longitude/latitude
; From pa,elo-90: z=Anti-solar, x=ecliptic north
; To   lng,lat  : x=Direction Sun, z=ecliptic north

rr = EulerRotate([0.0d0,90.0d0,0.0d0]/dpm,rr,degrees=degrees)

InitVar, nRR, 20
InitVar, dRR, 0.1

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

F = dblarr(nts_all,n_elements(elo_grid))


FOR i=0,nts_all-1 DO BEGIN

	IF running_diff THEN BEGIN

		tdiff = TimeOp(/subtract,ts_all[i],delta_time)

		dd = TimeOp(/subtract,ts_all,tdiff,uday)
		tmp = min(abs(dd),j)

		print, 'Running diff ',strjoin(TimeGet(ts_all[[i,j]],/ymd,upto=umin), ' - ')

	ENDIF


	; Get ecliptic longitude and latitude.
	; This is what EarthSky3DLoc converts to the heliographic
	; coordinate system used in the tomography.

	R = CvSky(ts_all[i],from_centered_ecliptic=rr,degrees=degrees,/to_ecliptic,/silent)

	; Convert to heliographic coordinates

	R = EarthSky3DLoc(ts_all[i]		, $
		nRA 		= 0 			, $ ; Forces use of input ra and dec
		nRR 		= nRR			, $
		dRR 		= dRR			, $
		ra			= reform(R[0,*]), $
		dec 		= reform(R[1,*]), $
		rr_earth	= rr_earth		, $
		ut_earth	= ut_earth		, $
		pa_earth	= pa_earth		, $
		elo_earth	= elo_earth 	, $
		elo_sun 	= elo_sun		, $
		degrees 	= degrees		, $
		/to_heliographic			, $
		zero_point	= zero_point	, $
		zero_phase	= zero_phase	, $
		zero_shift	= zero_shift	, $
		center		= center		)

	ctt = TimeGet(ts_all[i],/ymd,upto=umin)

	CASE readall OF
	0: BEGIN
		destroyvar, ffi
		hdri = hdr[i]
		utt = vu_getdata(hdri,ffi,silent=silent,_extra=_extra)
		IF running_diff THEN BEGIN
			hdrj = hdr[j]
			utt = vu_getdata(hdrj,ffj,silent=silent,_extra=_extra)
			ffi -= ffj
		ENDIF
	END
	1: BEGIN
		hdri = vu_gettime(hdr,ff,ut=ts_all[i],ff_ut=ffi)
		IF running_diff THEN BEGIN
			hdrj = vu_gettime(hdr,ff,ut=ts_all[j],ff_ut=ffj)
			ffi -= ffj
		ENDIF
	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)

	; 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[i,*] = 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[i,*] = 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[i,*] = sqrt(F[i,*])
	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]

		;CASE running_diff OF
		;0: tmp_ff = ffi[*,*,*,1]
		;1: tmp_ff = ffi[*,*,*,1]-ffi[*,*,*,1]
		;ENDCASE

		F[i,*] = 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[i,*] /= 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[i,*] = (It-Ir)/(Itb-Irb)
		;F[i,*] = ((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[i,*] = 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[i,*] /= 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[i,*] = 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

ENDFOR


old_device = !D.name

set_page, $
	zbuffer		= 1-printer		, $
	printer		= printer		, $
	xysize		= type.xysize	, $
	_extra		= _extra		, $
	old_device	= old_device	, $
	ctable		= type.ctable	, $
	/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

; F represents a "J-map" (time on the x-axis, vs.
; elongation on the y-axis).

PlotEloTimeMap, pa_angle, ts_all, elo_grid, f, $
	log 		= log			, $
	breakval	= BreakVal		, $
	charsize	= type.charsize , $
	degrees 	= degrees		, $
	title		= title 		, $
	format		= format		, $
	silent		= silent		, $
	_extra		= _extra

module = type.display+'_elotime'
vu_get_page, hdri, module	, $
	title	= ''			, $
	ut		= utt			, $
	device	= old_device	, $
	silent	= silent		, $
	_extra	= _extra


;whatis, f
;view, in=magnifyarray(f,2),/stretch,/new


RETURN  &  END