;+
; 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
; 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]: (optional) wavelength (m)
; 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
; OUTPUTS:
;	R		array[N,L,M]		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^-05 Integral( n_e[cm^-3] x H[nT].dl[AU] ) radians m^-2
;
;	The quantity returned by this function is the integrand
;
;	RM = 3.9362082 x 10^-05 n_e[cm^-3] x H[nT].dl[AU]            radians m^-2
;	   = 3.9362082 x 10^-05 n_e[cm^-3] x H[nT].dl[AU] x (180/pi) 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.
; 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).
;-
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. This sets the 3rd component 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 coordinate dimension from last to first dimension
; Then undo the normalization

tmp = transpose(F,[3,0,1,2])

tmp = SubArray(tmp, element=0, multiply=1/rr^4)
tmp = SubArray(tmp, element=1, multiply=1/rr^3)
tmp = SubArray(tmp, element=2, multiply=1/rr^2)

const = !sun.au*0.01d0*!physics.chargeofelectron^3/		$
	(2*!dpi*!physics.massofelectron^2*!physics.speedoflight^4)

; Calculate dot-product between los unit vectors and magnetic field
; This is the rotation measure.

tmp = -w3darg[0]*const*scalarproduct(P,tmp)*(180.0d0/!dpi)

; Multiply by wavelength^2 if present

IF n_elements(w3darg) GE 2 THEN tmp *= w3darg[1]*w3darg[1]

RETURN, tmp  &  END