;+
; NAME:
;	haystack_2006
; CALLING SEQUENCE:
	PRO haystack_2006
; INCLUDE:
	@compile_opt.pro
;-
restore, 'cr2006data.idl'

; UT Time

t = TimeSet(yr=2003, mon='aug', day=16, hr=16)


; Pick up density matrix

rho = data.rhoi

; Pick up x,y,z components of magnetic fields
; (presumable these are components in heliographic coordinate
; system with x-axis along heliographic longitude zero, and
; z-axis along solar rotation axis.

bb = SuperArray(rho, 3, /lead)

bb = SubArray(bb,element=0,replace=data.bxi)
bb = SubArray(bb,element=1,replace=data.byi)
bb = SubArray(bb,element=2,replace=data.bzi)

; Put dimension in order longitude, polar angle, radius
; (needed for IntegrateLOS which sums over the 3rd dimension)

rho = transpose(rho,[2,1,0])
bb  = transpose(bb ,[0,3,2,1])

; Pick up heliographic grid

xcgrid = reform(data.pi[0,0,*])	; Heliographic longitude (radians)
xlgrid = reform(data.ti[0,*,0])	; Polar angle (radians)
rrgrid = reform(data.ri[*,0,0])	; Heliocentric distance (solar radii)

xcgrid = AngleUnits(fromradian=xcgrid, /todegree)		; Heliographic longitude in degrees
xlgrid = 90.0d0-AngleUnits(fromradian=xlgrid, /todegree); Heliographic latitude in degrees
rrgrid = rrgrid*!sun.rau								; Heliocentric distance in AU

xlgrid = (xlgrid < 90.0d0) > (-90.0d0)

; Reverse polar angle grid, so that the grid runs from
; south to north pole (not necessary, but this is what I'm used too).

rho = reverse(rho,2)
bb  = reverse(bb ,3)
xlgrid = reverse(xlgrid)

; Line of sight to source.

los = AngleUnits( fromalmanac=[ [9,50,10,770], [14,19,59,20] ], /todegree )
sun = jpl_eph(t,jpl_body(/sun),center=jpl_body(/earth),/to_sphere,/deg,/location)
print, 'Source RA/dec: ',los
print, 'Sun    RA/dec: ',sun

; Convert to heliographic coordinates

unit = CvSky(t,fromequator=los,/toheliographic,/degrees)
sun  = CvSky(t,fromequator=sun,/toheliographic,/degrees)
print, 'Source heliographic: ',unit
print, 'Sun    heliographic: ',sun

; Unit vector along line of sight in heliographic coordinates

unit = cv_coord( from_sphere=[unit,1], /to_rect, /degrees )

; Magnetic field component along line of sight.

bpar = scalarproduct(bb,unit)

; Heliographic coordinates of segments along line of sight

nrr = 200
drr = 0.01

r = earthsky3dloc(t,nra=0,ra=los[0],dec=los[1],/degrees,/equator,nrr=nrr,drr=drr,/cv2carrington)

f = InterpolateHeliosphere(r,rho*bpar,xcgrid=xcgrid,xlgrid=xlgrid,rrgrid=rrgrid)

; Rotation measure in radians/m^2

v = 3.9269908*integratelos(r,f*drr)

;v = !physics.chargeofelectron^3/(2*!dpi*!physics.massofelectron^2*!physics.speedoflight^4)*1e69*v

print,v

RETURN  &  END