C+
C NAME:
C	ThomsonLOS
C PURPOSE:
C	Determines the integrated intensity along a line of sight for
C	electron Thomson scattering in a 1/r^2 density using Romberg integration.
C CALLING SEQUENCE:
	function ThomsonLOS(SLower,SUpper,ScSun,Elo,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	ScSun	real	heliocentric distance of observer (solar radii)
C	Elo	real	elongation of s/c line of sight (l.o.s.) in degrees.
C			(Elo=0 is the direction to the Sun)
C	U	real	limb darkening constant
C OUTPUTS:
C	F	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	ThomsonSetupIntegrand, nrQRomb, nrQRomo
C INCLUDE:
	include		'sun.h'
	include		'phys.h'
	include		'math.h'
C EXTERNAL:
	external	ThomsonTang
	external	ThomsonTangMRad
	external	MidInf
C RESTRICTIONS:
C >	For a 1/r^2 density the observer is specified by a heliocentric 
C	distance and the orientation of the line of sight by an elongation.
C >	For a user-specified density function (see Thomson3DDensity) the 
C	observer is specified in heliocentric spherical coordinates (usually
C	ecliptic coordinates): longitude, latitude and distance. The line of
C	sight orientation is specified in topocentric spherical coordinates
C	(again usually ecliptic coordinates): longitude and latitude only.
C	The topocentric longitude is defined relative to the Sun-observer
C	direction (see subroutine Point_On_LOS for more information).
C >	ThomsonLOSStep, ThomsonLOS3DStep: BE CAREFUL !!!
C	if a finite integration range [0,SUpper] is used to represent an
C	integration [0,infinity], the value of SUpper should be selected
C	carefully to avoid that the missing mass (between [SUpper,infinity])
C	introduces large errors. E.g. for ScSun < 1 AU, SUpper=2 AU introduces
C	errors of the order of 6-7% on some lines of light. For ScSun > 1 AU
C	this gets worse very rapidly.
C PROCEDURE:
C >	ThomsonLOSStep, ThomsonLOS3DStep:
C	calculates the integral by straightforward summing of contributions
C	from nStep intervals between SLower and SUpper.
C >	ThomsonLOS, ThomsonLOS3D:
C	calculates the integral using Romberg integration between SLower and
C	SUpper. If SUpper<0 then the integration is extended along the line of
C	sight until the integral has converged sufficiently. This effectively
C	represents an integration along the entire line of sight.
C >	ThomsonLOSFar:
C	calculates an analytical expression for the integral along the
C	entire line of sight for a 1/r^2 density in the limit of small angular
C	size of the Sun.
C
C >	The integrated intensity incident on the observer has cgs units of
C	erg/s/cm^2/sterad. The flux from the Sun incident on the observer has
C	cgs units of erg/s/cm^2, so the ratio will have units 1/sterad.
C >	denAU	electron density in the solar wind at 1 AU (cm^-3)
C >	In principle, the calculation of the density in the DO-loop can
C	be replaced by a call to a function calculating any desired radial
C	density profile.
C >	Coronal density		n(s) = denAU*(r0/r(s))^2,  r0 = 1 AU
C	Intensity/electron	I(s) = 0.5*Pi*Sigma*I0*Func(Omega,Chi)
C				(r and s are in units of AU)
C	The integral has the form
C		B = Integral[0,inf]{n(s)I(s)ds}
C		  = 0.5*Pi*Sigma*I0*Integral[0,inf]{n(s) Func(Omega,Chi) ds}
C	Func(Omega,Chi) is returned by ThomsonBase. The integral is calculated
C	by taking nStep steps from 0 to SUpper AU.
C >	The result is expressed in units of the flux received from the solar
C	disk by the observer:
C		F = (Pi*RSun^2)*I0*(1-U/3) / ScSun^2 (erg/s/cm^2)
C		  = Pi*I0*(1-U/3)* (RSun/ScSun)^2
C >	10^-16 = 10^-26*10^10
C	Factor 10^-26 is from the Thomson cross section, Sigma
C	Factor 10^10  is from the integration step size (in solar radii)
C MODIFICATION HISTORY:
C	1996, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
C-
	    real	SLower
	    real	SUpper
	    real	ScSun
	    real	Elo
	    real	U
	    real	P

	real		nrQRomb
	real		nrQRomo

	Elo = max(Elo,MATH__TINY)

	S0 = ThomsonSetupIntegrand(ScSun,Elo,U)

	!-------
	! Integral split the integration in two section: from SLower to SRomb, using
	! QRomb, and an integration of SRomb to SUpper, using nrQRomo with the MidInf
	! function.

	SRomb = 2.*max(1./sngl(SUN__RAU),ScSun*cosd(Elo))

	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(ThomsonTang    ,S0,S1)
	    rItr = rItr+nrQRomb(ThomsonTangMRad,S0,S1)
	    S0 = S1			! Limits for integration with QRomo
	    S1 = SHi
	end if

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

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

	return
	end