;+
; NAME:
;	ThomsonLOSDensity
; PURPOSE:
; CATEGORY:
;	sat/idl/toolbox/thomson
; CALLING SEQUENCE:
	FUNCTION ThomsonLOSDensity, S
; INPUTS:
;	S		array[k]; type: float
;				topocentric distance electron-observer (in solar radii)
; OUTPUTS:
;	F		array; type: float
;				density at electron location (in electrons/cm^-3)
; INCLUDE:
	@compile_opt.pro		; On error, return to caller
	@thomson_common.pro		; Dummy comment
; CALLS:
;	InitVar, SuperArray, SubArray, SyncArgs, CvPointOnLOS, IsType
; PROCEDURE:
; >	The observer locations and directions of the lines of sight are
;	first setup in a common block by href=ThomsonSetupLOS=.
; >	The distance along the lines of sight 'S', in combination with
;	the common block variables are used to get the heliocentric coordinates
;	of all line-of-sight locations.
; >	The function 'density' is a user-defined function, which returns the
;	density for a given heliocentric location.
;	The function is called using the IDL call_function routine.
; >	If no density function is specified than a 1/r^2 density with a
;	1 electrons/cm^3 density at 1 AU is used.
; MODIFICATION HISTORY:
;	JAN-1998, Paul Hick (UCSD)
;	JUN-2001, Paul Hick (UCSD/CASS)
;		Added /grid keyword
;	JUL-2008, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;		Removed /grid keyword again
;		Added some code to allow calculation to work for observers
;		out of the ecliptic (with i.e. with ecliptic latitude
;		not zero).
;-

; Copy Pos and Dir from common block in @thomson_common

pos = Pos_						; [3,npos], Observer location (lng,lat,dist)
los = Dir_						; [3,nlos], LOS direction (sun-centered lng,lat)

szpos = size(pos)
szlos = size(los)
szseg = size(S  )

npos = szpos[szpos[0]+2]/szpos[1]
nlos = szlos[szlos[0]+2]/szlos[1]
nseg = szseg[szseg[0]+2]
ntot = nseg*nlos*npos

abg = dblarr(3,npos)
abg[1,*] = pos[1,*]

pos = SuperArray(pos,nseg*nlos,after=1) 	 		; [3,nseg*nlos,npos]
pos = reform(pos,szpos[1],ntot,/overwrite)			; [3,nseg*nlos*npos]

abg = SuperArray(abg,nseg*nlos,after=1)
abg = reform(abg,3,ntot,/overwrite) 				; [3,nseg*nlos*npos]

los = SuperArray(los,nseg,after=1)					; [3,nseg,nlos]
los[2,*,*] = SuperArray(S,nlos,/trail)	 			; Line-of-sight distance
los = SuperArray(los,npos,/trail)					; [3,nseg,nlos,npos]

los = reform(los,szlos[1],ntot,/overwrite) 			; [3,nseg*nlos*npos]


; los contains spherical coordinates in a topocentric
; coordinate system 1 with
;	x-axis in ecliptic plane toward sun
;	z-axis to ecliptic north

; Define coordinate system 2 with
;	y-axis along cross product pos and ecliptic north (as in 1)
;	z-axis along cross product y and pos
;	x-axis points away from observer (opposite to pos).

; The Euler angles rotating from 1 to 2:
;	(0,lat,0) where lat = pos[1,*]

tmp = where(abg[1,*] NE 0)
IF tmp[0] NE -1 THEN los[*,tmp] = EulerRotate(abg[*,tmp],los[*,tmp])

; Convert LOS distance in units of observer-Sun distance

los[2,*] /= pos[2,*]

; Convert from topocentric to heliocentric coordinates

los = CvPointOnLOS(los)

; Convert distance back to physical units

los[2,*] *= pos[2,*]

; los now are heliocentric lng,lat in coordinate
; system 3 with
;	x-axis toward observer (pos)
;	z-axis as in system 2.

; Rotate from system 3 to heliocentric ecliptic coordinates
; with lng relative to observer pos.

IF tmp[0] NE -1 THEN los[*,tmp] = EulerRotate(-reverse(abg[*,tmp]),los[*,tmp])

; los[0,*] is a heliocentric longitude relative to the
; observer. Convert to absolute longitude by
; adding the longitude of the observer.

los[0,*] += pos[0,*]

; If no density function is specified a 1/r^2 density with
; density of 1 electron/cm^3 at 1 AU is used.

CASE IsType(density, /string) OF
0: los = 1.0/(reform(los[2,*])*!sun.RAu)^2	; [nseg*nlos*npos]
1: los = call_function(density, los)
ENDCASE

destroyvar, sz
IF szseg[0] GT 0 THEN boost, sz, szseg[1:szseg[0]]
IF szlos[0] GT 1 THEN boost, sz, szlos[2:szlos[0]]
IF szpos[0] GT 1 THEN boost, sz, szpos[2:szpos[0]]

CASE IsType(sz,/defined) OF
0: los = los[0]
1: los = reform(los, sz, /overwrite)
ENDCASE

RETURN, los  &  END