C+
C NAME:
C	ThomsonLOS3D
C PURPOSE:
C	Determines the integrated intensity along a line of sight for
C	electron Thomson scattering for a user-defined density using Romberg integration
C CALLING SEQUENCE:
	function ThomsonLOS3D(SLower,SUpper,ScLng,ScLat,ScSun,SLng,SLat,U, P)
C INPUTS:	(the input values are retained on return):
C	SLower	real	lower limit of integration (solar radii)
C	SUpper	real	upper limit of integration (solar radii)
C			If SUpper<0 the upper limit is effectively infinity:
C			      SUpper is set to the parameter MATH__PINF (=1.e30) solar radii
C	ScLng	real	heliocentric longitude of observer (deg)
C	ScLat	real	heliocentric latitude of observer (deg)
C	ScSun	real	heliocentric distance of observer (solar radii)
C	SLng	real	topocentric longitude line of sight relative to
C			Sun-observer direction
C	SLat	real	topocentric longitude line of sight relative to
C			Sun-observer direction
C	U	real	limb darkening constant
C OUTPUTS:
C	ThomsonLOS3D
C		real	Integrated Thomson scattering intensity received
C			per sterad in units of 10^-16 times the flux received
C			from the solar disk (at the observer location).
C	P	real	polarization
C CALLS:
C	ThomsonSetupLOS, ThomsonSetupIntegrand, nrQRomb, nrQRomo
C INCLUDE:
	include		'sun.h'
	include		'phys.h'
	include		'math.h'
C EXTERNAL:
	external	ThomsonTang3D
	external	ThomsonTangMRad3D
	external	MidInf
C RESTRICTIONS:
C PROCEDURE:
C	See href=ThomsonLOS=
C MODIFICATION HISTORY:
C	1996, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
C-
	    real	SLower
	    real	SUpper
	    real	ScLng
	    real	ScLat
	    real	ScSun
	    real	SLng
	    real	SLat
	    real	U
	    real	P

	real		nrQRomb
	real		nrQRomo

	Elo = acosd(cosd(SLng)*cosd(SLat))
	S0 = ThomsonSetupLOS(ScLng,ScLat,ScSun,SLng,SLat)
	S0 = ThomsonSetupIntegrand(ScSun,Elo,U)

	SRomb = 2.*max(1./sngl(SUN__RAU),ScSun*cosd(Elo))
					! Switch from QRomb to QRomo
	SHi = SUpper			! Upper limit
	if (SHi .lt. 0.) SHi = MATH__PINF! Improper integral

	S0 = max(SLower,0.)		! Lower limit
	S1 = min(SHi,SRomb)		! Don't integrate past SRomb with QRomb

	rIt  = 0.
	rItr = 0.
	if (S0 .lt. S1) then	 	! Near part on line of sight
	    rIt  = rIt +nrQRomb(ThomsonTang3D    ,S0,S1)
	    rItr = rItr+nrQRomb(ThomsonTangMRad3D,S0,S1)
	    S0 = S1			! Update integration limits
	    S1 = SHi
	end if

	if (S0 .lt. S1) then		! Far part on line of sight
	    rIt  = rIt +nrQRomo(ThomsonTang3D    ,S0,S1,MidInf)
	    rItr = rItr+nrQRomo(ThomsonTangMRad3D,S0,S1,MidInf)
	end if

	rI = rIt+rIt-rItr		! Total intensity
	P  = 0.
	if (rI .ne. 0.) P  = rItr/rI	! Polarization
					! Intensity/(Pi*I0)
	ThomsonLOS3D = ScSun*ScSun/(1-U/3)*0.5*PHYS__THOMSON*SUN__R*rI

	return
	end