C+
C NAME:
C	BuildSourceSurface
C PURPOSE:
C	Construct a new density/velocity map from line of sight values and weights based
C	on the ratio or difference of observed to model values for each line of sight. 
C CATEGORY:
C	Tomography
C CALLING SEQUENCE:
	subroutine BuildSourceSurface(I__VD,bRemove,bBound,BINXmin,NBAD,FIX,FIXMEAN,FIXSTDV,
     &		PPprj,Weight,XCbegMAT,XCendMAT,ZMAP,SIGB,Tmp1,Tmp2,BINX,BIN)
C INPUTS:
C	I__VD			integer		= TOM__V: new V map
C						= TOM__G: new G^2 map
C	bRemove			logical		activates sigma cut-off if .TRUE.
C	bBound                  logical
C
C	BINXmin			real		min. # los crossings required for deconvolution  
C	NBAD(NL)		integer		flag for bad lines of sight data (0-bad,1-good)
C	FIX (NL)		real		ratio observed/model los observation (from LOSTweak)
C	FIXMEAN			real
C	FIXSTDV			real
C
C	PPprj (3,NLOS,NL)	real		Lng/lat/time of los segment on source surface
C	Weight(NLOS,NL)		real		Weights of los segments (from href=LOSIntegralV= or href=LOSIntegralG=)
C
C	XCbegMAT(nTim)
C	XCendMAT(nTim)
C
C	ZMAP(nLng,nLat,nTim)	real		Velocity or g^2 map
C						(the input map is assumed to have no holes)
C	SIGB			real		sigma cut-off (applied to SIG)
C	bBound			logical		limit the D-map to its boundaries?
C
C	Tmp1(nLng,nLat,nTim)	real		scratch array
C	Tmp2(nLng,nLat,nTim)	real		scratch array
C	BIN (PRJ__N,NLOS,NL)	integer		scratch array
C OUTPUTS:
C	ZMAP(nLng,nLat,nTim)	real		new map; holes in the map are indicated by BadR4()
C	NBAD(NL)		integer		(Updated only if NS=1)
C						flag for bad line of sight data (0-bad,1-good)
C	BINX(nLng,nLat,nTim)	real		# projected los segments in each bin at the source surface
C
C	Not used anywhere as yet:
C
C	BIN (PRJ__N,NLOS,NL)	integer		longitude index of source surface map bin where segment was projected
C						(NOMOD only: will be zero if segment projected outside map)
C				integer		latitude index of bin where segment was projected
C				integer		time index of bin where segment was projected
C CALLS:
C	BadR4, ArrR4Zero, iArrR4ValuePresent, ArrR4Total, Say, Int2Str, Flt2Str, Str2Str
C	FLINT
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 RESTRICTIONS:
C	There should be no holes (BadR4()) values in the input array ZMAP
C PROCEDURE:
C	MODE	= iand(MODE,TOM__DIFF) = 0: corrections based on ratios
C		= iand(MODE,TOM__DIFF) > 0: corrections based on differences
C		= iand(MODE,TOM__MOD ) = 0: no remapping
C		= iand(MODE,TOM__MOD ) > 0: remap XC using mod(XCend-XCbeg)
C
C	iand(MODE,TOM__MOD)
C	> 0: segments with XCprj values outside [XCbeg,XCend] are mapped back into this range
C	= 0: segments with XCprj values outside [XCbeg,XCend] are ignored
C MODIFICATION HISTORY:
C	NOV-1995, B. Jackson (STEL,UCSD)
C	MAY-1999, Paul Hick (UCSD/CASS); extensive revisions
C	AUG-2001, Paul Hick (UCSD/CASS; pphick@ucsd.edu); modified for time-dependent tomography
C-
	    integer	I__VD
	    logical	bRemove
	    logical	bBound
	    real	BINXmin
	    integer	NBAD    (*)		! Bad source velocity 0 - bad, 1 - good
	    real	FIX     (*)		! Ratio observed/model G for total los
	    real	FIXMEAN
	    real	FIXSTDV

	    real	PPprj   (*)		! Lng/lat/time on source surface of point on los
	    real	Weight  (*)		! Weights of each point on los

	    real	XCbegMAT(*)
	    real	XCendMAT(*)

	    real	ZMAP    (*)		! Velocity or density map
	    real	SIGB
 	    real	Tmp1    (*)		! Scratch arrays
	    real	Tmp2    (*)
	    real	BINX    (*)		! # los crossings in each bin of the map

	    integer	BIN     (*)		! Scratch array

	character	cSay*18			/'BuildSourceSurface'/
	character	cStr*200
	integer		Str2Str
	integer		Flt2Str

	real		SIGOLD(2)		/-1.,-1./
	save		SIGOLD

	real		Tiny			/0.0/

	logical		bV
	logical		bMod
	logical		bDiff

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


	call T3D_iget(T3D__MODE,0,MODE)

	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)

	call T3D_iget(T3D__NLOS_V,0,NLOS)

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


	bV	= I__VD .eq. TOM__V
	bMod	= iand(MODE,TOM__MOD ) .ne. 0
	bDiff	= iand(MODE,TOM__DIFF) .ne. 0

	nSST	= nLng*nLat*nTim
	Bad	= BadR4()

        call ArrR4Zero(nSST,Tmp1)
        call ArrR4Zero(nSST,Tmp2)
        call ArrR4Zero(nSST,BINX)

	NSS = 0
	do I=1,NL					! Loop over all lines of sight

	    if (NBAD(I) .eq. 1) then
		if (abs(FIX(I)-FIXMEAN) .gt. SIGB*FIXSTDV) then	! Source above limit
		    NSS = NSS+1
		    if (bRemove) NBAD(I) = 0		! Mark source as bad
	    	end if
	    end if

	    if (NBAD(I) .eq. 1) then
		FixI = FIX(I)

		LastLON = 0
		LastLAT = 0
		LastTIM = 0


		do J=1,NLOS				! Step along line of sight

		    locPP = locLOS(J,I)			! Loc of los segment in array Weight

		    K = locPRJ(0,J,I)

		    locXC = K+PRJ__XC			! Lng  loc at source surface
		    locXL = K+PRJ__XL			! Lat  loc at source surface
		    locTT = K+PRJ__TT			! Time loc at source surface

		    F = PPprj( locXC )			! Projected longitude
		    T = PPprj( locTT )			! Projected time

		    XCbeg   = FLINT(1,nTim,XCbegMAT,T,Tiny)! Needed for statement fncs XCindx, XCmod
		    XCend   = FLINT(1,nTim,XCendMAT,T,Tiny)
		    XCrange = XCend-XCbeg

		    if (bMod) F = XCmod(F)		! Rescale to range [XCbeg,XCend]
 
		    NearLON = nint(XCindx(F))		! Nearest longitude index
		    NearTIM = nint(TTindx(T))		! Nearest time index

		    !-------
		    ! Skip segments outside map

		    if (NearLON .lt. 1 .or. NearLON .gt. nLng .or. NearTIM .lt. 1 .or. NearTIM .gt. nTim) then
			BIN   ( locXC ) = 0
			Weight( locPP ) = 0.		! Mark skipped segment

		    else
			NearLAT = nint(XLindx(PPprj( locXL ))) ! Nearest latitude index

			BIN( locXC ) = NearLON
			BIN( locXL ) = NearLAT
			BIN( locTT ) = NearTIM

			K = locSS(NearLON,NearLAT,NearTIM)

			!-------
			! Several neighboring segments on a los may end up in the same bin
			! at the source surface. Count only the first of these (i.e. when the
			! los moves into another bin)

			if (J .eq. 1 .or. NearLON .ne. LastLON .or.
     &				NearLAT .ne. LastLAT .or. NearTIM .ne. LastTIM) BINX(K) = BINX(K)+1.0

			LastLON = NearLON		! Line of sight crosses bin NearLON, NearLAT, NearTIM
			LastLAT = NearLAT
			LastTIM = NearTIM

			if (bBound .and. (PPprj( locXC ) .lt. XCbeg .or.
     &					  PPprj( locXC ) .gt. XCend)) Weight( locPP ) = 0.

			W = Weight( locPP )
			F = FixI

			Tmp1(K) = Tmp1(K)+W*F	 	! Sum of correction factors
			Tmp2(K) = Tmp2(K)+W		! Sum of weights

		    end if
 		end do
	    end if
	end do

	if (bMod) then
	    do K=1,nTim
		do J=1,nLat
		    I1 = locSS(1,J,K)
		    I2 = I1+nLng1			!locSS(nLng,J,K)

		    Tmp1(I1) = Tmp1(I1)+Tmp1(I2)
		    Tmp1(I2) = Tmp1(I1)
		    Tmp2(I1) = Tmp2(I1)+Tmp2(I2)
		    Tmp2(I2) = Tmp2(I1)
		    BINX(I1) = BINX(I1)+BINX(I2)
		    BINX(I2) = BINX(I1)
		end do
	    end do
	end if

	do I=1,nSST					! Calculate new map
	    if (BINX(I) .lt. BINXmin) then		! Data not deconvolved: bad value in map
		Tmp1(I) = Bad
		ZMAP(I) = Bad
	    else 
		Tmp1(I) = Tmp1(I)/Tmp2(I)		! Weighted average correction
		if (bDiff) then				! New good value
		    ZMAP(I) = max(0.,ZMAP(I)+Tmp1(I))
		else
		    ZMAP(I) = Tmp1(I)*ZMAP(I)
		end if
	    end if
	end do

	nGood = nSST-iArrR4ValuePresent(nSST,ZMAP,Bad)! # good values
	if (nGood .eq. 0) call Say(cSay,'E','Empty','deconvolved '//char(I__VD)//' map')

	!-------
	! Nothing of what follows leaves this subroutine. So it doesn't matter what is done exactly.
	! In particular SIGMAP is not (yet?) used to terminate the iteration process.

	call ArrR4Zero(nSST,Tmp2)			! Calculate Chi square across map

	W = 0.
	do I=1,NL
	    if (NBAD(I) .ne. 0) then
		FixI = FIX(I)

		do J=1,NLOS

		    locXC = locPRJ(PRJ__XC,J,I)		! Lng  loc at source surface
		    NearLON = BIN( locXC )

		    if (NearLON .ne. 0) then		! Segment inside map

			locXL = locXC-PRJ__XC+PRJ__XL	! Lat  loc at source surface
			locTT = locXC-PRJ__XC+PRJ__TT	! Time loc at source surface

			NearLAT = BIN( locXL )
			NearTIM = BIN( locTT )

			locNear = locSS(NearLON,NearLAT,NearTIM)

			if (ZMAP( locNear ) .ne. Bad) then
			    F = FixI-Tmp1( locNear )

			    locPP = locLOS(J,I)		! Loc of los segment in array Weight

			    W = W+Weight( locPP )
			    Tmp2( locNear ) = Tmp2( locNear )+Weight( locPP )*(F*F)
			end if
		    end if
		end do
	    end if
	end do

	SIGMAP = sqrt(ArrR4Total(nSST,Tmp2,I)/W)	! Calculate average Chi square for whole map
							! For bad ZMAP values Tmp2 is zero
	!-------
	! Write message with convergence information to screen

	IVD = 1
	if (bV) IVD = 2

	I = 0
	I = I+Int2Str(NSS,			cStr(I+1:)) ! Build string for Say call
	I = I+Str2Str('/',			cStr(I+1:))
	I = I+Int2Str(NL ,			cStr(I+1:))
	I = I+Str2Str(' (',			cStr(I+1:))
	I = I+Flt2Str(100.*NSS/NL,2,		cStr(I+1:))
	I = I+Str2Str(' %) sources above',	cStr(I+1:))+1
	I = I+Flt2Str(SIGB,1,			cStr(I+1:))
	I = I+Str2Str(' sigma',			cStr(I+1:))
	if (bRemove) I = I+Str2Str(' REMOVED',	cStr(I+1:))
	I = I+Str2Str('#map convergence:',	cStr(I+1:))+1
	I = I+Flt2Str(SIGMAP,3,			cStr(I+1:))
	if (SIGOLD(IVD) .ne. -1.) then
	    I = I+Str2Str(' (',			cStr(I+1:))
	    I = I+Flt2Str(SIGMAP/SIGOLD(IVD),3, cStr(I+1:))
	    I = I+Str2Str(' x previous value)',	cStr(I+1:))
	end if

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

	SIGOLD(IVD) = SIGMAP				! Saved for comparison on next pass

	return
	end