C+
C NAME:
C	LOSReach
C PURPOSE:
C	Does a sanity check on the los segment locations output by LOSProjection
C CATEGORY:
C	Data processing
C CALLING SEQUENCE:
	subroutine LOSReach(I__VD,XCbegMAT,XCendMAT,PPlos,NBS)
C INPUTS:
C	MODE		integer		= iand(MODE,TOM__MOD) > 0: remap XC using mod(XCend-XCbeg)
C					= iand(MODE,TOM__MOD) = 0: no remapping
C	NCoff		integer		Carrington rotation offset to be applied to
C					XCbeg, XCend and XClos
C
C	TT		real		start time as Carrington variable (- NCoff)
C	dTT		real		time resolution as Carrington variable
C					(if nTim = 1 then dTT should be set to zero)
C	RR		real		heliocentric distance of source surface (AU) 
C	dRR		real		radial distance resolutions (AU)
C
C	nLng		integer		# longitudes
C	nLat		integer		# latitudes
C	nRad		integer		# heliocentric distances
C	nTim
C
C	XCbegMAT(nTim)	real		start Carrington variable
C	XCendMAT(nTim)	real		end Carrington variable
C
C	NL		integer		# IPS G-level observations
C	NLOS		integer		# segments of length dLOS along los
C	PPlos(4,NLOS,NL)real		Lng/lat/rad/time of point on los
C OUTPUTS:
C	NBS(NL)		integer		iand(MODE,TOM__MOD) > 0: all set to 1
C					iand(MODE,TOM__MOD) = 0: set to 1 if whole line of sight
C					insided range [XCbeg,XCend]; otherwise set to 0.
C					(los with NBS=0 will be discarded later by href=LOSClean=)
C	iMap		integer		Lowest heliocentric distance index (<=nRad) above
C					all los segments
C CALLS:
C	T3D_iget, T3D_get_grid, T3D_iset, pInfR4, ArrR4Total, cFlt2Str, cInt2Str
C	iwhitespace, Say, ArrI4Constant, iArrI4Total
C INCLUDE:
	include		't3d_param.h'
	include		't3d_array.h'
	include		't3d_index.h'
C	include		't3d_grid_fnc.h'
C	include		't3d_loc_fnc.h'
C PROCEDURE:
C >	XClos and XCprj are on the same scale as XCbeg and XCend.
C >	XClos, XLlos, RRlos are output from href=LOSPosition=.
C >	iand(MODE,TOM__MOD) > 0:
C		Program is terminated if XLlos or RRlos values lie outside the volume
C		covered by the reconstruction. XClos is allowed to fall outside [XCbeg,XCend]
C		(these will be mapped inside this range by GetLOSValue, when necessary).
C >	iand(MODE,TOM__MOD) = 0:
C		Program is terminated if XClos, XLlos or RRlos values lie outside the volume.
C >	Effectively LOSReach makes sure that all segments lie inside the reconstruction volume
C >	Note that segments can still end up outside [XCbeg,XCend] after traceback to the
C	source surface.
C MODIFICATION HISTORY:
C	SEP-1999, P. Hick (UCSD/CASS; pphick@ucsd.edu)
C-
	    integer	I__VD
	    real	XCbegMAT(*)
	    real	XCendMAT(*)
	    real	PPlos   (*)	! Lng/lat/rad/time of point on los
	    integer	NBS     (*)

	character	cSay*8		/'LOSReach'/
	character	cStr*120
	character	cFlt2Str*14
	character	cInt2Str*14

	logical		bMod
	integer		nRadSave(2)	/2*0/
	save		nRadSave

	include		't3d_grid_fnc.h'
	include		't3d_loc_fnc.h'

	if (I__VD .eq. TOM__V) call T3D_iget(T3D__NL_V,0,NL)
	if (I__VD .eq. TOM__G) call T3D_iget(T3D__NL_G,0,NL)

	if (NL .eq. 0) return

	call T3D_iget(T3D__MODE,0,MODE)
	bMod = iand(MODE,TOM__MOD) .ne. 0


	call T3D_iget(T3D__NLOS_V,0,NLOS )
	call T3D_iget(T3D__NCOFF,0,NCoff)

	call T3D_get_grid(TT,dTT,RR,dRR, nLng,nLat,nRad,nTim,dTTi,nLng1,nLat1)


	Big = pInfR4()

	XClosMin =  Big
	XClosMax = -Big
	XLlosMin =  Big
	XLlosMax = -Big
	RRlosMin =  Big
	RRlosMax = -Big
	TTlosMin =  Big
	TTlosMax = -Big

	do K=1,NL
	    do J=1,NLOS
		L = locPOS(0,J,K)

		XClosMin = min(XClosMin,PPlos(L+LOS__XC))
		XClosMax = max(XClosMax,PPlos(L+LOS__XC))
		XLlosMin = min(XLlosMin,PPlos(L+LOS__XL))
		XLlosMax = max(XLlosMax,PPlos(L+LOS__XL))
		RRlosMin = min(RRlosMin,PPlos(L+LOS__RR))
		RRlosMax = max(RRlosMax,PPlos(L+LOS__RR))
		TTlosMin = min(TTlosMin,PPlos(L+LOS__TT))
		TTlosMax = max(TTlosMax,PPlos(L+LOS__TT))
	    end do
	end do


	XCbeg = ArrR4Total(nTim,XCbegMAT,J)/nTim
	XCend = ArrR4Total(nTim,XCendMAT,J)/nTim
	XCrange = XCend-XCbeg

	cStr =	' Carrington variable : range '//cFlt2Str(NCoff+XClosMin   ,-3)//', '//cFlt2Str(NCoff+XClosMax,-3)//
     &		'# ..... reconstruction range '//cFlt2Str(NCoff+XCvar(nLng),-3)//', '//cFlt2Str(NCoff+XCvar(1),-3)
	I = iwhitespace(cStr)

	call Say(cSay,'I',char(I__VD),cStr)

	cStr =	 'heliographic latitude: range '//cFlt2Str(XLlosMin,-2)//', '//cFlt2Str(XLlosMax   ,-2)//
     &		'# ...... reconstruction range '//cFlt2Str(XLdeg(1),-2)//', '//cFlt2Str(XLdeg(nLat),-2)
	I = iwhitespace(cStr)

	call Say(cSay,'I',char(I__VD),cStr)

	cStr =	 'heliocentric distance: range '//cFlt2Str(RRlosMin ,-4)//', '//cFlt2Str(RRlosMax    ,-4)//
     &		'# ...... reconstruction range '//cFlt2Str(RRhght(1),-4)//', '//cFlt2Str(RRhght(nRad),-4)
	I = iwhitespace(cStr)

	call Say(cSay,'I',char(I__VD),cStr)

	cStr =	 'Carrington time      : range '//cFlt2Str(NCoff+TTlosMin       ,-3)//', '//cFlt2Str(NCoff+TTlosMax    ,-3)//
     &		'# ...... reconstruction range '//cFlt2Str(NCoff+TTtime(1),-3)//', '//cFlt2Str(NCoff+TTtime(nTim),-3)
	I = iwhitespace(cStr)

	call Say(cSay,'I',char(I__VD),cStr)


	if (.not. bMod .and. (XClosMin .lt. XCvar(nLng) .or. XClosMax .gt. XCvar(1))) then
	    do I=1,NL
		L = locPOS( LOS__XC,1,I )

		XClosMin =  Big
		XClosMax = -Big

		do J=1,NLOS
		    L = locPOS( LOS__XC,J,I )

		    XClosMin = min(XClosMin,PPlos(L))
		    XClosMax = max(XClosMax,PPlos(L))
		end do

		if (XClosMin .ge. XCvar(nLng) .and. XClosMax .le. XCvar(1)) then
		    NBS(I) = 1
		else
		    NBS(I) = 0
		end if

	    end do
	    cStr = cInt2Str(NL-iArrI4Total(NL,NBS,I))//'/ '//cInt2Str(NL)//'los extend outside Carrington range'
	    I = iwhitespace(cStr)

	    call Say(cSay,'W',char(I__VD),cStr)

	else
	    call ArrI4Constant(NL,1,NBS)

	end if

	if (XLlosMin .lt. XLdeg(1) .or. XLlosMax .gt. XLdeg(nLat))
     &		call Say(cSay,'F',char(I__VD),'los segments outside permitted heliographic latitude range')

	if (RRlosMin .lt. RRhght(1) .or. RRlosMax .gt. RRhght(nRad))
     &		call Say(cSay,'F',char(I__VD),'los segments outside permitted heliocentric distance range')

	if (dTT .ne. 0 .and. (TTlosMin .lt. TTtime(1) .or. TTlosMin .gt. TTtime(nTim)))
     &		call Say(cSay,'F',char(I__VD),'los segments outside permitted time range')

	if (I__VD .eq. TOM__V) nRadSave(1) = min(int(RRindx(RRlosMax))+1,nRad)
	if (I__VD .eq. TOM__G) nRadSave(2) = min(int(RRindx(RRlosMax))+1,nRad)

	if (nRadSave(1) .gt. 0 .and. nRadSave(2) .gt. 0) then
	    nRad = max(nRadSave(1),nRadSave(2))
	    call T3D_iset(T3D__NRAD,0,nRad)
	end if

	return
	end