C+
C NAME:
C	FillMapOLS
C PURPOSE:
C	To fill maps at opposite the region of interest with a smooth average of the 
C	data from the other nearby longitudes. 
C CATEGORY:
C	Data processing
C CALLING SEQUENCE:
C	call FillMapOLS(Mo,XCbe,XCint,nLng,nLat,nT,N,nTmax,RRS,Speed,Const,AEmax,AEmin,AEang,DMap)
C INPUTS:
C	Mo		   integer	1 - Smooth/interpolate region most distant from Earth
C					2 - Limit values farthest from Earth
C	XCbe(2,nT)	   real		Boundary values of time maps
C	XCint(NTmax)	   real		Start and end times of intervals
C	nLng		   integer	# Longitude points
C	nLat		   integer	# Latitude points
C	nT		   integer	# Time points
C	N		   integer	Map time value
C	nTmax		   integer	Total # of map time values
C	RRS		   real		Deconvolution source surface distance
C	Speed		   real		Average speed used for data analysis
C	Const		   real		Longitude filter constant (in terms of the numbers of intervals in one rotation)
C	AEmax		   real		Anti-Earth maximum value
C	AEmin		   real		Anti-Earth minimum value
C	AEang		   real		Angle from Earth to begin longitude limit
C	DMap(nLng,nLat,nT) real		Input map
C OUTPUTS:
C	DMap(nLng,nLat,nT) real		Output map
C CALLS:
C PROCEDURE:
C	Bad values (indicated by BadR4()) are processed
C MODIFICATION HISTORY:
C	MAY, 1999 B. Jackson (STEL,UCSD)
C-
	subroutine FillMapOLS(Mo,XCbe,XCint,nLng,nLat,nT,N,NTmax,RRS,Speed,Const,AEmax,AEmin,AEang,DMap)

	real	XCbe(2,nT)			! Boundary values of time maps
	real	XCint(nTmax)			! Start and end times of intervals		
	real	DMap(nLng,nLat,nT)		! Map

	! Average the other non-Bad longitude values.
	
	XCmid = (XCint(N)+XCint(N+1))/2.0
	if (Speed .ne. 0.0) XCmid = XCmid-(0.8-RRS)*((400.0/Speed)*4.0)/25.38

	XCbeg = XCbe(1,N)
	XCend = XCbe(2,N)

	Bad = BadR4()
	ExpPWC = 50.0

	NXC  = XCend-XCbeg + 0.1		! Usually 3
	Lmod = nLng/NXC				! Usually 18 for the IPS
	IXCopp = (XCend-XCmid)*Lmod+0.5		! Middle I of region of interest

	! This part filters or smooths the data

	if (Mo .eq. 1) then

	    nLngB1 = IXCopp-Lmod/2
	    nLngE1 = nLngB1+1
	    nLngE2 = IXCopp+Lmod/2
	    nLngB2 = nLngE2-1
	    NX = 2

	    !print *,'NXC,Lmod,XCbeg,XCend,XCmid,IXCopp,nLngB1,nLngE1'
	    !print *, NXC,Lmod,XCbeg,XCend,XCmid,IXCopp,nLngB1,nLngE1,nLngB2,nLngE2

	    do I=nLngB1,nLngE1

		do J=1,nLat
		    DMNT   = 0.0
		    ANDMNT = 0.0
		    do II=1,NX
			IP = I + II
			IM = I - II
			if (1 .le. IP .and. IP .le. nLng) then
			    if (DMap(IP,J,N) .ne. Bad) then
				ExpPW = (II/Const)**2
				Expot = 999999.0
				if (ExpPW .gt. ExpPWC) Expot = 0.0
				if (Expot .ne. 0.0) Expot = exp(-ExpPW)
				DMNT   = DMNT  +Expot*DMap(IP,J,N)
				ANDMNT = ANDMNT+Expot
			    end if
			end if
			if (1 .le. IM .and. IM .le. nLng) then
			    if (DMap(IM,J,N) .ne. Bad) then
				ExpPW = (II/Const)**2
				Expot = 999999.0
				if (ExpPW .gt. ExpPWC) Expot = 0.0
				if (Expot .ne. 0.0) Expot = exp(-ExpPW)
				DMNT   = DMNT  +Expot*DMap(IM,J,N)
				ANDMNT = ANDMNT+Expot
			    end if
			end if
		    end do

		    if (1 .le. I .and. I .le. nLng) then
			if (ANDMNT .ne. 0.0 .and. DMap(I,J,N) .eq. Bad) DMap(I,J,N) = DMNT/ANDMNT
		    end if

		end do

	    end do

	    do I=nLngB2,nLngE2

		do J=1,nLat
		    DMNT   = 0.0
		    ANDMNT = 0.0
		    do II=1,NX
			IP = I + II
			IM = I - II
			if (1 .le. IP .and. IP .le. nLng) then
			    if (DMap(IP,J,N) .ne. Bad) then
				ExpPW = (II/Const)**2
				Expot = 999999.0
				if (ExpPW .gt. ExpPWC) Expot = 0.0
				if (Expot .ne. 0.0) Expot = exp(-ExpPW)
				DMNT   = DMNT  +Expot*DMap(IP,J,N)
				ANDMNT = ANDMNT+Expot
			    end if
			end if
			if (1 .le. IM .and. IM .le. nLng) then
			    if (DMap(IM,J,N) .ne. Bad) then
				ExpPW = (II/Const)**2
				Expot = 999999.0
				if (ExpPW .gt. ExpPWC) Expot = 0.0
				if (Expot .ne. 0.0) Expot = exp(-ExpPW)
				DMNT   = DMNT  +Expot*DMap(IM,J,N)
				ANDMNT = ANDMNT+Expot
			    end if
			end if
		    end do

		    if (1 .le. I .and. I .le. nLng) then
			if (ANDMNT .ne. 0.0 .and. DMap(I,J,N) .eq. Bad) DMap(I,J,N) = DMNT/ANDMNT
		    end if

		end do

	    end do

	end if

	! This part limits the Anti-Earth maps beyond an angle AEang

	if (Mo .eq. 2) then

	    Idist = Lmod*AEang/360.0+0.5		! Integer distance of anti-Earth angle
	    Ibeg = IXCopp - Idist
	    Iend = IXCopp + Idist

	    !print *, 'XCmid, XCbeg, IXCopp, Idist, Ibeg, Iend', XCmid, XCbeg, IXCopp, Idist, Ibeg, Iend

	    do I=1,nLng
		if (I .lt. Ibeg .or. I .gt. Iend) then
		    do J=1,nLat
			if (DMap(I,J,N) .ne. Bad) DMap(I,J,N) = max(AEmin,min(DMap(I,J,N),AEmax))
		    end do
		end if
	    end do

	end if

	return
	end