;+
; NAME:
;	vu_lineofsight
; PURPOSE:
; CATEGORY:
;	user/phick/idl
; CALLING SEQUENCE:
	FUNCTION vu_lineofsight, ff	, $
		xcgrid		= xcgrid	, $
		xlgrid		= xlgrid	, $
		rrgrid		= rrgrid	, $
		utt 		= utt		, $
		ra			= ra		, $
		dec 		= dec		, $
		equator 	= equator	, $
		nrr 		= nrr		, $
		drr 		= drr		, $
		degrees 	= degrees	, $
		ips_speed	= ips_speed , $
		ips_glevel	= ips_glevel, $
		thomson_b	= thomson_b , $
		thomson_pb	= thomson_pb, $
		column_density=column_density, $
		faraday_rm	= faraday_rm, $
		betan		= betan 	, $
		betar		= betar 	, $
		silent		= silent	, $
		_extra		= _extra
; INPUTS:
;	ff				array[n,m,l,k]; type: float
;						k=2,3,4 or 5
;						3D volume data on grid specified in
;						a spherical coordinate system
;					 ff[*,*,*,0]: solar wind speed (km/s)
;					 ff[*,*,*,1]: normalized density (cm^-3); n x (r/r0)^2
;					 ff[*,*,*,2]: ? Bx (not sure what coordinate frame)
;					 ff[*,*,*,3]: ? By
;					 ff[*,*,*,4]: ? Bz
;	xcgrid=xcgrid	array[n]; longitudinal grid
;	xlgrid=xlgrid	array[m]; latitudinal grid
;	rrgrid=rrgrid	array[l]; radial grid
;
;	ra				array[*]; ra/ecl long of line(s) of sight
;	dec 			array[*]; decl/ecl latituded of line(s) of sight
; OPTIONAL INPUT PARAMETERS:
;	utt=utt			array[1]; type: time structure
;						UT time assigned to volume data ff
;	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)
;
;	/ips_speed		IPS speed
;	/ips_glevel 	IPS g-level
;	/thomson_b		Thomson scattering brightness
;	/thomson_pb 	Thomson scattering polarization brightness
;
; OUTPUTS:
;	(none)
; 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

InitVar, betan	, 0.3
InitVar, betar	, 0.5
InitVar, n0		, 5.0

InitVar, nRR	, 20				; Nr. segments along los
InitVar, dRR	, 0.1				; Step size along los

; The input ra and dec are either in equatorial coordinates or in
; ecliptic coordinates.

; The output R vectors (locations of all line-of-sight segments
; is in the coordinate frame selected through one of the to_*
; keyword to CvSky (the default is ecliptic coordinates).
; The to_* keywords also affect rr_rearth

R = EarthSky3DLoc(utt		, $
	nRA 		= 0 		, $ 	; Forces use of input ra and dec
	ra			= ra		, $ 	; Los ra
	dec 		= dec		, $ 	; Los decl
	nrr 		= nrr		, $ 	; Los # segments
	drr 		= drr		, $ 	; Los step size
	equator 	= equator	, $ 	; Equatorial/ecliptic ra,dec?
	degrees 	= degrees	, $

	rr_earth	= rr_earth	, $
	ut_earth	= ut_earth	, $
	pa_earth	= pa_earth	, $
	elo_earth	= elo_earth	, $
	elo_sun 	= elo_sun	, $

	_extra		= _extra	)		; Used to pass to_* to CvSky

nagoya  =  1
fillbad =  0
minRR	= 10


CASE 1 OF
IsType(ips_speed	 , /defined): id_type = vu_los_ips_speed
IsType(ips_glevel	 , /defined): id_type = vu_los_ips_glevel
IsType(thomson_b 	 , /defined): id_type = vu_los_thomson_b
IsType(thomson_pb	 , /defined): id_type = vu_los_thomson_pb
IsType(column_density, /defined): id_type = vu_los_column_density
IsType(faraday_rm	 , /defined): id_type = vu_los_faraday_rm
ELSE							: id_type = vu_los_nodata
ENDCASE

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

CASE id_type OF

vu_los_nodata: IF silent LE 1 THEN message, /info, 'no los integration selected'

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

	D3D = ff[*,*,*,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, ff[*,*,*,0], D3D, $
		xcgrid		= xcgrid			, $
		xlgrid		= xlgrid			, $
		rrgrid		= rrgrid			, $
		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 = ff[*,*,*,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 			, $
		xcgrid		= xcgrid			, $
		xlgrid		= xlgrid			, $
		rrgrid		= rrgrid			, $
		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, ff[*,*,*,1] 	, $
		xcgrid		= xcgrid			, $
		xlgrid		= xlgrid			, $
		rrgrid		= rrgrid			, $
		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 = F/ThomsonRadialFilter(elo_earth[*,*,0], degrees=degrees)

	;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			, $
		xcgrid		= xcgrid			, $
		xlgrid		= xlgrid			, $
		rrgrid		= rrgrid			, $
		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, ff[*,*,*,1]	, $
		xcgrid		= xcgrid			, $
		xlgrid		= xlgrid			, $
		rrgrid		= rrgrid			, $
		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			, $
		xcgrid		= xcgrid			, $
		xlgrid		= xlgrid			, $
		rrgrid		= rrgrid			, $
		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, ff[*,*,*,1]	, $
		xcgrid		= xcgrid			, $
		xlgrid		= xlgrid			, $
		rrgrid		= rrgrid			, $
		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			)

	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, ff[*,*,*,1] 	, $
		xcgrid		= xcgrid			, $
		xlgrid		= xlgrid			, $
		rrgrid		= rrgrid			, $
		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)

	;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 = ff[*,*,*,2:4]*SuperArray(ff[*,*,*,1],3,/trail)

	F = IntegrateLOS(R, D3D 			, $
		xcgrid		= xcgrid			, $
		xlgrid		= xlgrid			, $
		rrgrid		= rrgrid			, $
		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)

END

ENDCASE

RETURN, F  &  END