;+
; NAME:
;	vu_lineofsight
; PURPOSE:
; CATEGORY:
;	user/phick/idl
; CALLING SEQUENCE:
	PRO vu_lineofsight, hdr, ff	, $
		hdr_bb		= hdr_bb			, $
		bb			= bb				, $
		ut0			= ut0				, $
		ihdr		= ihdr				, $

		filter_bb	= filter_bb			, $

		ephemeris	= ephemeris			, $

		RA			= RA				, $
		DEC 		= DEC				, $
		nRR			= nRR				, $
		dRR			= dRR				, $

		equator 	= equator			, $
		degrees 	= degrees			, $
		silent		= silent			, $
		_extra		= _extra
; INPUTS:
;	hdr			array; 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.
;	ra				array[n];  ra/ecl long of line(s) of sight
;	dec 			array[n]; decl/ecl latitude of line(s) of sight
;	ephemeris		array[2,n]; alternative way of specifying ra,dec
;						(if used, this overrides ra,dec) 
; OPTIONAL INPUT PARAMETERS:
;	ut0=ut0			array[1]; type: time structure
;						UT time 
;	nrr=nrr 		scalar; type: integer
;						nr of segments along line(s) of sight
;	drr=drr 		scalar; type: float
;						stepsize along line(s) of sight (AU)
;	/degrees		if set, indicates that all angles are in degrees. Default: radians.
;	/equator		if set, ra and dec are assumed to be equatorial
;						coordinates (default is ecliptic long and latitude)
;
;	/sw_density		solar wind density
;	/sw_speed		solar wind speed
;	/ips_speed		IPS speed
;	/ips_glevel 	IPS g-level
;	/thomson_b		Thomson scattering brightness
;	/thomson_pb 	Thomson scattering polarization brightness
;
; OUTPUTS:
;	Print to standard output
; INCLUDE:
	@compile_opt.pro			; On error, return to caller
	@vu_fnc_include.pro 		;
; EXTERNAL:
;	IPS_velocity, IPS_velocity_w, IPS_hlevel_w, ThomsonBrightness
; CALLS:
;	InitVar, IsType
;	EarthTransit3DLoc, EarthSky3DLoc, IPS_params, SuperArray
;	ThomsonUBVConst, IntegrateLOS, AngleRange
;	ThomsonRadialFilter, TimeGet, TimeUnit
; PROCEDURE:
; >	Several parameters define the grid used to calculate the sky map.
;	These are passed to href=EarthSky3DLoc= and href=EarthTransit3DLoc=:
;
;	nRR		scalar; type: int; default: 20
;				# steps along lines of sight at each point in the sky
;	dRR		scalar; type: float; default: 0.1 AU
;				step size along line of sight (AU)
;
; MODIFICATION HISTORY:
;	FEB-2007, Paul Hick (UCSD/CASS)
;-

InitVar, silent , 0
InitVar, equator, /key

ihdr_set = IsType(ihdr, /defined)

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

NeedBB = 1
HaveBB = IsType(hdr_bb,/defined) AND IsType(bb,/defined)

IF 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, silent=silent)
		HaveBB = count_bb EQ count
	ENDIF

	IF NOT HaveBB THEN BEGIN
		say, 'no magnetic field files found', threshold=2, silent=silent
    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

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

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)

nv_norm = vu_get(hdri, /power_norm)
bb_norm = HaveBB ? vu_get(hdri_bb, /power_norm) : [0,0,0]
IF n_elements(bb_norm) EQ 2 THEN BEGIN
	Bn = SubArray(bbi, /lastdim,  element=2)
	bmin = min(Bn, /nan, max=bmax)
	IF bmin NE 0 OR bmax NE 0 THEN BEGIN
		say, 'non-zero Bn found, but normalization is not specified in header'
		RETURN
	ENDIF
	bb_norm = [bb_norm,0]		; Value doesn't matter since Bn=0
ENDIF

InitVar, nRR, 20
InitVar, dRR, 0.1

; The input ra and dec are either in equatorial coordinates
; (/equator set) or in ecliptic coordinates (/equator NOT set)..

; The output R vectors (locations of all line-of-sight segments
; The return array R has longitude angles in R[0,*] (not Carrington variables)

IF IsType(ephemeris,/defined) THEN BEGIN
	RA  = SubArray(ephemeris, element=0)
	DEC = SubArray(ephemeris, element=1)
ENDIF

say, 'input ephemeris is in '+(equator ? 'equatorial' : 'ecliptic')+' coordinates'

; Get the locations along the line of sight:
; (heliographic longitude and latitude, heliocentric distance)

R = EarthSky3DLoc(utt		, $
	nRA 		= 0 		, $ 	; Forces use of input ra and dec
	nRR 		= nRR		, $ 	; Los # segments
	dRR 		= dRR		, $ 	; Los step size
	RA			= RA		, $ 	; Los ra
	DEC			= DEC		, $		; Los dec
	rr_los		= rr_los	, $
	rr_earth	= rr_earth	, $
	ut_earth	= ut_earth	, $
	pa_earth	= pa_earth	, $
	elo_earth	= elo_earth	, $
	elo_sun 	= elo_sun	, $
	equator 	= equator	, $ 	; Equatorial/ecliptic ra,dec?
	degrees 	= degrees	, $
	/to_heliographic		)

nagoya  =  1
fillbad =  0
minRR	= 0;10

rr_sun = SubArray(R,dim=0,element=2)

; Velocity and density along line of sight; remove normalization

VN = InterpolateHeliosphere(R, ffi, R3D, dR3D, xcrange=xcrange, /cv2carrington, degrees=degrees, fillbad=fillbad)
FOR i=0,1 DO IF nv_norm[i] NE 0 THEN VN = SubArray(VN, /lastdim, element=i, divide=rr_sun^nv_norm[i])

; Velocity perpendicular to line of sight

VVperp = IPS_velocity(rr_earth, R, SubArray(VN,/lastdim,element=0), elo_earth=elo_earth, elo_sun=elo_sun, degrees=degrees)

; Magetic field component along line of sight; remove normalization

BRTN = HaveBB ? InterpolateHeliosphere(R, bbi, R3D, dR3D,xcrange=xcrange, /cv2carrington, degrees=degrees, fillbad=fillbad) : replicate( 0.0, [(size(R,/dim))[1:*],3] )
FOR i=0,2 DO IF bb_norm[i] NE 0 THEN BRTN = SubArray(BRTN, /lastdim,  element=i, divide=rr_sun^bb_norm[i])	; Remove normalization

; Magnetic field component along line of sight 

BLOS = RotationMeasure(rr_earth, R, BRTN, replicate(0,6), ut_earth=ut_earth, pa_earth=pa_earth, elo_earth=elo_earth, elo_sun=elo_sun, degrees=degrees, /get_blos)

; IPS velocity

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]

Vips = IntegrateLOS(R, ffi[*,*,*,0], D3D, $
	R3D			= R3D				, $
	dR3D		= dR3D				, $
	xcrange		= xcrange			, $
	/cv2carrington					, $
	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			)

; G-level

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

Glevel = IntegrateLOS(R, D3D 		, $
	R3D			= R3D				, $
	dR3D		= dR3D				, $
	xcrange		= xcrange			, $
	/cv2carrington					, $
	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			)

Glevel = sqrt(Glevel)

; Thomson scattering brightness

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

Thomson_B = IntegrateLOS(R, ffi[*,*,*,1], $
	R3D			= R3D				, $
	dR3D		= dR3D				, $
	xcrange		= xcrange			, $
	/cv2carrington					, $
	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 			)

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


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

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

Itb = IntegrateLOS(R, Fb			, $
	R3D			= R3D				, $
	dR3D		= dR3D				, $
	xcrange		= xcrange			, $
	/cv2carrington					, $
	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]	, $
	R3D			= R3D				, $
	dR3D		= dR3D				, $
	xcrange		= xcrange			, $
	/cv2carrington					, $
	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,nv_norm[1]]

Irb = IntegrateLOS(R, Fb			, $
	R3D			= R3D				, $
	dR3D		= dR3D				, $
	xcrange		= xcrange			, $
	/cv2carrington					, $
	minRR		= minRR 			, $
	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]	, $
	R3D			= R3D				, $
	dR3D		= dR3D				, $
	xcrange		= xcrange			, $
	/cv2carrington					, $
	minRR		= minRR 			, $
	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			)

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

; Column density

;ColumnDensity = IntegrateLOS(R, ffi[*,*,*,1] 	, $
;	R3D			= R3D				, $
;	dR3D		= dR3D				, $
;	xcrange		= xcrange			, $
;	/cv2carrington					, $
;	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			)

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

D3D = (HaveBB ? bbi : 0)*SuperArray(ffi[*,*,*,1],3,/trail)
f3darg = [dRR,BadValue(dRR),nv_norm[1],bb_norm]

RM = IntegrateLOS(R, D3D 	, $
	R3D			= R3D				, $
	dR3D		= dR3D				, $
	xcrange		= xcrange			, $
	/cv2carrington					, $
	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)

FOR n=0,n_elements(RA)-1 DO BEGIN	; All lines of sight

	print
	print, "Time                 =", ctt
	IF equator THEN BEGIN
		print, "RA,Dec               =", RA[n],',',DEC[n]
	ENDIF ELSE BEGIN
		print, "Ecl lng,lat          =", RA[n],',',DEC[n]
	ENDELSE

	print, "Elongation           =", elo_earth[n,0]
	print, "PA angle             =", pa_earth[n,0]
	point_p = rr_earth[2]*sin(elo_earth[n,0]*ToRadians(degrees=degrees))
	print, "Point-P              =", point_p, ' AU (',point_p/!sun.rau,' Rsun)'
	print, "IPS velocity         =", Vips  [n], ' km/s'
	print, "IPS g-level          =", Glevel[n]
	print, "Thomson scattering B =", Thomson_B[n], ' S10'
	print, "Thomson scattering pB=", Thomson_pB[n]
	;print, "Column density=", ColumnDensity[n]
	print, "Rotation measure     =", RM[n], ' deg/m^2'
	print

	print, "  d      hlng    hlat    R          V       Vperp       n          Br        Bt        Bn       B||"
	FOR i=0,nRR-1 DO BEGIN		; Steps along line of sight
		print, format='(F5.2,3F8.3,2F10.2,F10.2,4F10.2)', rr_los[i],R[*,n,i],VN[n,i,0],VVperp[n,i],VN[n,i,1],BRTN[n,i,*],BLOS[n,i]
	ENDFOR

ENDFOR

RETURN  &  END