;+
; NAME:
;	RotationMeasure
; PURPOSE:
;	Calculates component of vector field n*B along line of sight
; CATEGORY:
;	sat/idl/toolbox/faraday
; CALLING SEQUENCE:
	FUNCTION RotationMeasure, REarth, R, F, w3darg, $
		ut_earth	= ut_earth	, $
		pa_earth	= pa_earth	, $
		elo_earth	= elo_earth , $
		elo_sun 	= elo_sun	, $
		degrees 	= degrees	, $
		get_blos	= get_blos 
; INPUTS:
;	REarth						not used
;	R			array[3,N,L,M]	heliographic coordinates of points in a 2D grid
;								of NxL lines of sight with M segments along each
;								R[0,*,*,*] : heliographic longitude
;								R[1,*,*,*] : heliographic latitude
;								R[2,*,*,*] : heliocentric distance (AU)
;	F			array[N,L,M,3]	magnetic field vector, multiplied by the
;								density, in RTN coordinates
;	w3darg						w3darg[0]: los step size (AU)
;								w3darg[1]: wavelength (m)
;								w3darg[2]: density normalization (default: 2)
;								w3darg[3]: Br normalization (default: 2)
;								w3darg[4]: Bt normalization (default: 1)
;								w3darg[5]: Bn (default: 1)
; OPTIONAL INPUT PARAMETERS:
;	/degrees					if set all input angles are in degrees
;	pa_earth	array[N,L,M]	position angle from ecliptic north
;	elo_earth	array[N,L,M]	angle Sun-Earth-LOS segment (i.e. conventional
;								elongation angle)
;	elo_sun		array[N,L,M]	angle Earth-Sun-LOS segment
;	/get_blos	if set, then the F-component along the line of sight is
;				returned. Note that this can be used to get the magnetic
;				field component along the line of sight by setting F to the
;				magnetic field (omitting multiplication with density), and
;				setting w3dard[2]=0.
; OUTPUTS:
;	R			array[N,L,M]	rotation measure in degrees/m^2; see PROCEDURE	
; INCLUDE:
	@compile_opt.pro		; On error, return to caller
; CALLS:
;	ToRadians, SubArray, CvSky, EulerRotate, scalarproduct
; PROCEDURE:
;	RM = e^3/(2*pi*m^2*c^4) Integral( n_e * H.dl )	radians cm^-2
;	In Gaussian units:
;		e = 4.8 x 10^-10 esu
;		m = 9.1 x 10^-28 gr
;		c = 3.0 x 10^-10 cm/s
;		n_e in cm^-3; H in Gauss; dl in cm
;
;	RM = 2.6311927 x 10^-17 Integral( n_e x H.dl ) radians cm^-2
;
;	The tomography magnetic files contain magnetic field in nT = 10^-5 Gauss
;	The distance scale for the integration is AU = 1.5 x 10^13 cm
;	The wavelengths are in the meter range, so we change
;	to radians m^-2. Multiplying with the resulting
;
;	10^-5 x 1.5 x 10^13 x 10^4 = 1.5 x 10^12:
;
;	RM = 3.9362082 x 10^-5 Integral( n_e[cm^-3] x H[nT].dl[AU] ) radians m^-2
;
;	The integrand is returned here in degrees m^-2:
;
;	RM = 3.9362082 x 10^-5 n_e[cm^-3] x H[nT].dl[AU]            radians m^-2
;
;	RM = 3.9362082 x 10^-5 n_e[cm^-3] x H[nT].dl[AU] x (180/pi) degrees m^-2
;	   = 2.2552811 x 10^-3 n_e[cm^-3] x H[nT].dl[AU]            degrees m^-2
;
;	The constant 3.9362082 x 10^-05 is calculated from fields in the
;	!sun and !physics structures:
;		3.9362082 x 10^-05 = !sun.au*0.01d0*!physics.chargeofelectron^3/	$
;			(2*!dpi*!physics.massofelectron^2*!physics.speedoflight^4)
;
;	The integration step size dl[AU] is stored in w3darg[0]
;	If w3darg[1] exists it should contain the wavelength in meters.
;	In this case RM*w3darg[1]^2 is returned, i.e. the Faraday rotation
;	in degrees.
;
;	A typical integration will extend over 2 AU from Earth. So in order of
;	magnitude the integrated RM will be:
;
;	RM = 0.045 n_e[cm^-3] x H[nT] degrees m^-2
;
;	Typical background values are n_e ~ 10-20 cm^-3 and H ~ 1-2 nT, so typical
;	background RM values will be in the range 0.5 to 2 degrees m^-2
;
; MODIFICATION HISTORY:
;	OCT-2006, Paul Hick (UCSD/CASS)
;	JUN-2007, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;		Bug identified by Liz Jensen. The integration step size needs to
;		be absorbed in the return value.
;		Changed sign on return value (effectively this changes the direction
;		of the integration, which now runs in the direction of the wave vector
;		of the incoming radio waves (toward Earth).
;-

InitVar, get_blos, /key
rpm = ToRadians(degrees=degrees)

pa  =  pa_earth*rpm 			; Position angle from ecl north for all los
elo = elo_earth*rpm 			; Elongations of all los

sinpa  = sin(pa )
cospa  = cos(pa )
sinelo = sin(elo)
coselo = cos(elo)

; Set up unit vectors along lines of sight in sun-centered
; ecliptic coordinates.
; Array P has a leading dimension of 3 (for three components
; in spherical coordinates). All components are initialized to
; one. Only the third is kept to one (length of a unit
; vector). The first two are overwritten below.

; The vectors P are unit vectors (in spherical coordinates) in
; a geocentric coordinate system with the x-axis towards the Sun
; and the z-axis to ecliptic north.

P = replicate(1+0*pa[0], [3,size(pa,/dim)])

P = SubArray(P,element=0,replace=atan(sinpa*sinelo,coselo)/rpm)	; Ecliptic longitude relative to Sun
P = SubArray(P,element=1,replace=asin(cospa*sinelo)/rpm)		; Ecliptic latitude

; Rotate line-of-sight unit vectors from sun-centered ecliptic to
; heliographic coordinates:
;
; Think of the unit vectors P shifted from the Earth to the Sun;
; the coordinates in P then become spherical coordinates in a
; coordinate system with the Sun in the origin, the x-axis
; pointing away from Earth, and the z-axis to ecliptic north.
; These heliocentric vectors are then converted to heliographic
; coordinates in the next statement.

P = CvSky(ut_earth, from_centered_ecliptic=P, /to_heliographic, degrees=degrees, /silent)

; Get Euler angles for rotation from heliographic to RTN coordinates
; at each line of sight segment.
; Note that the RTN frame depends on the heliospheric position vector
; R of each line-of-sight segment.

tmp = CvSky_RTN(R, degrees=degrees)

; Rotate line-of-sight unit vectors to RTN coordinates

P = EulerRotate(tmp,P,degrees=degrees)

; Convert from spherical to rectangular RTN coordinates

tmp = size(P,/dim)

P = reform(P, 3, n_elements(P)/3, /overwrite)
P = cv_coord(from_sphere=P, /to_rect, degrees=degrees)
P = reform(P, tmp, /overwrite)

; F is the product of two normalized quantities: density and magnetic field
; - density (r^2)
; - radial magnetic field (r^2)
; - tangential magnetic field (r^1)
; - normal magnetic field (r^0) ??? Or maybe not???

rr = SubArray(R, element=2) 		; Sun-Electron distance in AU

; Move the component dimension (3 magnetic field components)
; from last to first dimension. Then undo the normalization

n = size(F, /n_dim)-1
tmp = transpose(F,[n,indgen(n)])

nn_norm = n_elements(w3darg) GT 2 ? w3darg[2] : 2
br_norm = n_elements(w3darg) GT 3 ? w3darg[3] : 2
bt_norm = n_elements(w3darg) GT 4 ? w3darg[4] : 1
bn_norm = n_elements(w3darg) GT 5 ? w3darg[5] : 1

IF nn_norm+br_norm NE 0 THEN tmp = SubArray(tmp, element=0, divide=rr^(nn_norm+br_norm))
IF nn_norm+bt_norm NE 0 THEN tmp = SubArray(tmp, element=1, divide=rr^(nn_norm+bt_norm))
IF nn_norm+bn_norm NE 0 THEN tmp = SubArray(tmp, element=2, divide=rr^(nn_norm+bn_norm))

; Calculate dot-product between los unit vectors and nB
; This is the magnetic field magnitude along the los, multiplied by density

tmp = scalarproduct(P,tmp)

IF NOT get_blos THEN BEGIN
	const = !sun.au*0.01d0*!physics.chargeofelectron^3/(2*!dpi*!physics.massofelectron^2*!physics.speedoflight^4)

	tmp = -w3darg[0]*const*tmp*(180.0d0/!dpi)	; Rotation measure in degrees/m^2.
							; Multiply by wavelength^2 if present
	IF n_elements(w3darg) GE 1 THEN IF finite(w3darg[1]) THEN tmp *= w3darg[1]*w3darg[1]
ENDIF

RETURN, tmp  &  END