;+
; NAME:
;	ThomsonLOSRomb
; PURPOSE:
;	Calculates integrated line-of-sight intensity for
;	Thomson scattering using Romberg integration.
; CATEGORY:
;	Physics: Thomson scattering
; CALLING SEQUENCE:
	FUNCTION ThomsonLOSRomb, pos, los, P, $
		lower	= lower 	, $
		upper	= upper 	, $
		density = density	, $
		degrees = degrees	, $
		au		= au		, $
		udark	= udark 	, $
		apm 	= apm		, $
		s10 	= S10
; INPUTS:
;	pos			array[3,n]; type: float
;					heliocentric location of observer:
;					(ecliptic) longitude and latitude; and heliocentric distance
;					(a scalar 'pos' is interpreted as [0,0,pos])
;	los			array[2,m]; type: float
;					topocentric (ecliptic) longitude and latitude of line of sight
;					relative to Sun-observer direction
;					(a scalar 'los' is interpreted as [los,0])
; OPTIONAL INPUT PARAMETERS:
;	/s10		if set intensities are returned in S10 units
;					In this case the apparent magnitude APM MUST BE SPECIFIED
;	udark=udark	scalar; type: float; default: 0
;					limb darkening constant
;	apm=apm		scalar; type: float
;					apparent magnitude of Sun
;	density=density
;				scalar; type: string; default: undefined
;					name of function used to calculate the electron density.
;					Keyword is passed to href=ThomsonSetupLOS=.
;	lower=lower
;				scalar; type: float; default: 0
;					Lower limit of los integration (solar radii)
;					If Lower <=0 then Lower = 0 is used.
;	upper=upper
;				scalar; type: float; default: 9000 solar radii
;					Upper limit of los integration (solar radii)
;					If Lower <=0 then Upper = 9000.0 solar radii (~40 AU) is used.
;
;	/degrees	if set all angles are in degrees (default: radians)
;	/au			if set all distances are in AU (default: solar radii)
;
; OUTPUTS:
;	Result		scalar[n,m]; type: float
;					Integrated Thomson scattering intensity;
;					units depend on setting of /S10:
;					/S10 NOT set: per sterad in units of 10^-16 times the flux received
;						from the solar disk (at the observer location
;					/S10 set: S10 units
;	P[n,m]		scalar; type: float
;					polarization
; INCLUDE:
	@compile_opt.pro		; On error, return to caller
; CALLS:
;	InitVar, ToRadians, sphere_distance, ThomsonS10, ThomsonSetupLOS
;	ThomsonSetupRomb, SyncArgs, ThomsonSolarFlux, ToSolarRadii
; EXTERNAL:
;	ThomsonTang, ThomsonTangMRad
; SEE ALSO:
;	ThomsonLOSStep
; PROCEDURE:
; >	The lower and upper limits can actually be 1-dim arrays, but that is
;	probably not useful.
; >	If no 'density' function is specified then a 1/r^2 density is used with
;	density of 1 at 1 AU.
; MODIFICATION HISTORY:
;	SEP-1999, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;-

InitVar, udark, 0
InitVar, S10, /key


IF S10 THEN $
	IF IsType(APM, /undefined) THEN $
		message, 'no apparent magnitude for the Sun specified'
		
IF n_elements(pos) EQ 1 THEN pos = [0,0,pos]
IF n_elements(los) EQ 1 THEN los = [los,0]

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

npos = szpos[szpos[0]+2]/szpos[1]
nlos = szlos[szlos[0]+2]/szlos[1]

IF npos*nlos GT 1 THEN BEGIN

	pos = reform(pos,szpos[1],npos)
	los = reform(los,szlos[1],nlos)

	B = fltarr(nlos,npos)
	P = fltarr(nlos,npos)

	FOR ipos=0,npos-1 DO BEGIN
		FOR ilos=0,nlos-1 DO BEGIN

			B[ilos,ipos] = ThomsonLOSRomb(pos[*,ipos], los[*,ilos], tmp, $
				lower	= lower 	, $
				upper	= upper 	, $
				density = density	, $
				degrees = degrees	, $
				au		= au		, $
				udark	= udark 	, $
				apm 	= apm		, $
				s10 	= s10)
			P[ilos,ipos] = tmp
		ENDFOR
	ENDFOR

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

	CASE IsType(tmp,/defined) OF
	0: BEGIN
		B = B[0]
		P = P[0]
	END
	1: BEGIN
		B = reform(B, tmp, /overwrite)
		P = reform(p, tmp, /overwrite)
	END
	ENDCASE

	SyncDims, pos, size=szpos
	SyncDims, los, size=szlos

	RETURN, B

ENDIF

Tiny = 1.0E-6
jmax = 100

rpm  = ToRadians(degrees=degrees)
rpau = ToSolarRadii(au=au)

; Direction to sun (in same sun-centered coordinate system as los,
; i.e. the sun is at lng=0 at the opposite latitude)

tmp = pos[0:1]
tmp[0]  =  0.0d0
tmp[1] *= -1.0d0

elo  = sphere_distance(tmp*rpm, los*rpm)		; Elongation (radians)
RSun = pos[2]*rpau								; Observer-Sun distance (solar radii)

ThomsonSetupLOS , pos*[rpm,rpm,rpau], los*rpm, density=density
ThomsonSetupRomb, RSun, elo, udark

SRomb = 2.0d0*((RSun*cos(elo)) > (1.0d0/!sun.RAu))

InitVar, lower, 0.0
InitVar, upper, 9000.0/rpau

Hi = upper*rpau
L  = where(upper LT 0)
IF L[0] NE -1 THEN Hi[L] = infinitevalue(Hi)

S0 = lower*rpau > 0.0		; Lower limit
S1 = Hi < SRomb				; Don't integrate past SRomb with QRomb

It  = 0.0
Itr = 0.0
SyncArgs, S0, S1, It, Itr

L = where(S0 LT S1)
IF L[0] NE -1 THEN BEGIN	; Near part on line of sight
	It [L] += QRomb('ThomsonTang'    ,S0[L],S1[L],eps=Tiny,jmax=jmax)
	Itr[L] += QRomb('ThomsonTangMRad',S0[L],S1[L],eps=Tiny,jmax=jmax)
	S0 [L]  = S1 [L]		; Update integration limits
	S1 [L]  = Hi [L]
ENDIF

L = where(S0 LT S1)
IF L[0] NE -1 THEN BEGIN	; Far part on line of sight
	It [L] += QRomo('ThomsonTang'	 ,S0[L],S1[L],eps=Tiny,jmax=jmax,/midinf)
	Itr[L] += QRomo('ThomsonTangMRad',S0[L],S1[L],eps=Tiny,jmax=jmax,/midinf)
ENDIF

It = It+It-Itr				; Total intensity
P  = Itr/It					; Polarization




cv = 0.5d0*!physics.thomsoncrosssection

; Replacing RSun_ by r0 if S10=1 doesn't change the result, but
; is faster when RSun is an array.

r0 = 1.0d0/!sun.RAu

CASE S10 OF
0: cv /= 							ThomsonSolarFlux(RSun,udark)	; Intensity/Solar Flux
1: cv *= 1.0d-16*ThomsonS10(r0,apm)/ThomsonSolarFlux(r0  ,udark)	; S10 units
ENDCASE


RETURN,	cv*!sun.R*It  &  END