C+
C NAME:
C	MkDMap
C PURPOSE:
C	Change the density map from the line of sight (L.O.S.) values of G-level and
C	weights by the ratio of observed to model values for each line of sight.
C CATEGORY:
C	Data processing
C CALLING SEQUENCE:
	subroutine MkDMap(LON,LAT,GWTij,PW,FIX,NL,NLOS,NLOSP1,
     &		nLng,nLat,DMAP,ANLIMIT,GSSIG,SIGB,NS,
     &		NBADSG,SIGMAP,ANMAP,DMA,FIXM,WTSM,DMEANF,DFIXM2,WTP,CONR)  
C INPUTS:
C	LON(NLOSP1,NL)		integer		Projected Carr. rot. value of point on L.O.S.
C	LAT(NLOSP1,NL)		integer		Heliographic lat. point on L.O.S.
C	GWTij(NLOS,NL)		real		m Weights of each point along the L.O.S.
C	PW			real		Power of G to N
C	FIX(NL)			real		Ratio observed/model G for total L.O.S.
C	NL			integer		Number of G-level data points
C	NLOS			integer		Number of L.O.S. distance segments
C	NLOSP1			integer		Number of L.O.S. distance segments + 1
C	nLng			integer		Number of longitude points in DMAP
C	nLat			integer		Number of latitude points in DMAP
C	DMAP(nLng,nLat)		real		Density map
C	SIGMAP			real		Density map convergence criteria
C	ANLIMIT			real		Number of points constituting deconvolution  
C	GSSIG(NL)		real		Deviation of source from mean in sigma
C	SIGB			real		G-level sigma cut-off
C	NS			integer		1 = Activate sigma cut-off, 0 = do not
C	NBADSG(NL)		integer		Bad source G-value 0 - bad, 1 - good
C	CONR			real		Gridsphere filter factor
C SCRATCH ARRAYS:
C	ANMAP(nLng,nLat)	real		Scratch array
C	DMA  (nLng,nLat)	real		Scratch array
C	FIXM (nLng,nLat)	real		Scratch array
C	WTSM (nLng,nLat)	real		Scratch array
C	DMEANF(nLng,nLat)	real		Scratch array
C	DFIXM2(nLng,nLat)	real		Scratch array
C	WTP(NLOS,NL)		real		Scratch array
C OUTPUTS:
C	DMAP(nLng,nLat)		real		Density map
C	BADD			real		Value of map where there is no projection
C	NBADSG(NL)		integer		Bad source G-value 0 - bad, 1 - good
C	SIGMAP			real		Density map convergence criteria
C CALLS:
C	BadR4, ArrR4Bad, ArrR4Zero, GridSphere2D, Say
C PROCEDURE:
C MODIFICATION HISTORY:
C	NOV, 1995 B. Jackson (STEL,UCSD)
C-
	    integer	LON  (NLOSP1,NL)	! Projected Carr. rot. value of point on L.O.S.
	    integer	LAT  (NLOSP1,NL)	! Heliographic lat. (in deg.) point on L.O.S.
	    real	GWTij(NLOS  ,NL)	! m weights of each point along the L.O.S.
	    real	FIX  (       NL)	! Ratio observed/model G for total L.O.S.
	    real	GSSIG(       NL)	! Deviation of source from mean in sigma
	    real	DMAP  (nLng,nLat)	! Density map
	    integer	NBADSG(      NL)	! Bad source velocity 0 - bad, 1 - good
	    real	ANMAP (nLng,nLat)	! Number of points per bin
	    real	DMA   (nLng,nLat)
	    real	FIXM  (nLng,nLat)
	    real	WTSM  (nLng,nLat)
	    real	DMEANF(nLng,nLat)
	    real	DFIXM2(nLng,nLat)
	    real	WTP  (NLOS,NL)

	character	cStr*150

	include		'mapangles.h'

	Bad = BadR4()

	nLngLat = nLng*nLat
	nLOSnL  = NLOS*NL
	nLat1   = nLat - 1
	daLat2  = 90./nLat1

	call ArrR4Bad	(nLngLat,ANMAP)
	call ArrR4Bad	(nLngLat,DMA)
	call ArrR4Zero	(nLngLat,FIXM)
	call ArrR4Bad	(nLngLat,WTSM)
	call ArrR4Zero	(nLngLat,DMEANF)
	call ArrR4Zero	(nLngLat,DFIXM2)
	call ArrR4Zero	(nLOSnL,WTP)

	NSS = 0

	do I=1,NL
	    if (abs(GSSIG(I)) .gt. SIGB) then		! Source above limit
		NSS = NSS+1
		if (NS .eq. 1) NBADSG(I) = 0		! Mark source as bad
	    end if
	    if (NBADSG(I) .ne. 0) then			! Only use valid sources
		LNlast = 0
		LTlast = 0
		do J=1,NLOS
		    LN = LON(J,I)			! map longitude
		    LT = LAT(J,I)			! map latitude
		    WTSj = GWTij(J,I)
		    FIXM(LN,LT) = FIXM(LN,LT)+WTSj*FIX(I)
		    WTP(J,I) = WTSj
		    if(WTSM(LN,LT) .eq. Bad) WTSM(LN,LT) = 0.
		    WTSM(LN,LT) = WTSM(LN,LT)+WTP(J,I)
		    !print *, LN, LT, FIXM(LN,LT), WTSj, FIX(I)
		    if (LN .ne. LNlast .or. LT .ne. LTlast) then
			if(ANMAP(LN,LT).eq.Bad) ANMAP(LN,LT) = 0.
			ANMAP(LN,LT) = ANMAP(LN,LT)+1.
		    end if		
		    LNlast = LN
		    LTlast = LT
		end do
	    end if
	end do

	call GridSphere2D(ALng,nLng,nLat,1,ANMAP,CONR,3,0.0,90.0)

	!-------
	! Now modify the latitude map count to approximate the cos latitude filter effect

	do J=1,nLat
	    alat = abs(XLdeg(J))
	    if(alat.eq.90.) alat = alat - dalat2
	    sqonebycos = SQRT(1./cosd(alat))
	    do I=1,nLng
		if(ANMAP(I,J).ne.Bad) then
		    ANMAP(I,J) = ANMAP(I,J)*sqonebycos
		end if
	    end do		    
	end do

	do LN=1,nLng
	    do LT=1,nLat
		if (FIXM(LN,LT) .gt. 0. .and. WTSM(LN,LT) .gt. 0.) then
		    DMEANF(LN,LT) = FIXM(LN,LT)/WTSM(LN,LT)
		    DMA(LN,LT) = DMEANF(LN,LT)*DMAP(LN,LT)
		    !print *, LN, LT, FIXM(LN,LT), WTSM(LN,LT), DMEANF(LN,LT), DMAP(LN,LT), DMA(LN,LT)
		else
		    ANMAP(LN,LT) = 0.
		end if
	    end do
	end do

	!-------
	! calculation of the Chi square distribution for each latitude and longitude

	do I=1,NL
	    if (NBADSG(I) .ne. 0) then			! Only use valid sources
	    do J=1,NLOS
		LN = LON(J,I)			! map longitude
		LT = LAT(J,I)			! map latitude
		if (DMA(LN,LT) .ne. Bad) DFIXM2(LN,LT) = DFIXM2(LN,LT)+WTP(J,I)*((FIX(I)-DMEANF(LN,LT))**2)
	    end do
	    end if
	end do

	!-------
	! calculation of the "total" Chi square distribution for all latitudes and longitudes

	CONVER  = 0
	NCONVER = 0
	do I=1,nLng
	    do J=1,nLat
		if (ANMAP(I,J) .lt. ANLIMIT) DMA(I,J) = Bad ! data not deconvolved
		if (DMA(I,J) .ne. Bad) then
		    CONVER = CONVER+DFIXM2(I,J)
		    NCONVER = NCONVER+1
		end if
		DMAP(I,J) = DMA(I,J)
		!if(DMAP(I,J).gt.1.0e7.and.DMAP(I,J).lt.1.0e22) print *, I, J, DMAP(I,J)
	    end do
	end do

	SIGOLD = SIGMAP
	SIGMAP = CONVER/NCONVER

	write (cStr,'(A,F4.1,A,F5.2,A)') 'density sources above',SIGB,' sigma:',100.*NSS/NL,'%'
	if (NS .eq. 1) write (cStr(itrim(cStr)+1:),'(A,I4,A)') '#you just removed',NSS,' sources'
	write (cStr(itrim(cStr)+1:),'(A,F7.3,A,F7.3,A)')
     &		'#D map convergence =',SIGMAP,', from last =',abs(1-SIGMAP/SIGOLD)*100,'%'

	call Say('MkDMap','I','Info',cStr)

	return
	end