C+
C NAME:
C	MkDMaptdn0n
C PURPOSE:
C	Change the density map from the line of sight (L.O.S.) values of G-level and/or
C	Thomson-scattering weights by the ratio of observed to model values for each line 
C	of sight.  Originally the value of PW was included in this to include in the weighting
C	for g-level data.  This option is not now used, and the version of this subroutine for
C	g-level and for Thomson-scattering brightness is now identical. 
C	A big step taken here was to make real and smooth the numbers of points that 
C	constitude convergence. This smoothing is now done using gridsphere and timesmooth.
C	The reasoning behind this is that if the values of a map are smoothed over space and time
C	shouldn't the criteria that determines whether convergence is possible also be smoothed
C	in a similar way.  Thus, 1.5 lines of sight crossing can constitute convergence just as 
C	2 or 3 lines of sight crossing can.  In reality many lines of sight contribute to a 
C	given point on the map (maybe not 100%, each), so it stands to reason that the convergence
C	criterion need not be an integer 2 or 3 or whatever.  This also allows "high resolution"
C	cadences and spatial scales to be explored digitally while still retaining the same 
C	approximate spatial and temporal structure resolution.  It works!  It makes the whole 
C	iteration procedure more stable and gives better answers when there are few data points. 
C	I also tried to smooth the weights and changes using this idea, but I got worse answers,
C	at least in the test cases I tried.   Since this did not seem to help significantly 
C	because of the computational increase, I removed the weight and changes smoothing portion 
C	of the code from the analysis.  Another innovation was to increase the (now real) convergence 
C	criterion with latitude according to 1/cos.  This seems appropriate to me since the spatial 
C	map scales hung together with gridsphere are actually smaller at high latitudes.  Thus, fewer
C	lines of sight crossings are necessary at large latitudes. 
C   
C CATEGORY:
C	Data processing
C CALLING SEQUENCE:
C	call MkDMaptdn0n(LON,LAT,ITIM,GWTij,PW,FIX,NL,NLOS,NLOSP1,
C		nLng,nLat,nT,DMAP,ANLIMIT,GSSIG,SIGB,NS,NBADSG,SIGMAP,
C		ANMAP,AWTMAP,DMA,FIXM,WTSM,DMEANF,DFIXM2,WTP,Alng,ConstL,CONR,CONT,aNday,Ztmp)
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	ITIM(NLOSP1,NL)		integer		Projected times at reference surface
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       nT                      integer         Number of time intervals
C	DMAP(nLng,nLat,nT)	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 = don't
C	NBADSG(NL)		integer		Bad source G-value 0 - bad, 1 - good
C OUTPUTS:
C	DMAP(nLng,nLat,nT)	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	AWTMAP(nLng,nLat,nT)	real		Combination Wt and LOS crossing number map
C
C SCRATCH ARRAYS:
C	ANMAP(nLng,nLat,nT)	real
C	DMA(nLng,nLat,nT)	real
C	FIXM(nLng,nLat,nT)	real
C	WTSM(nLng,nLat,nT)	real
C	DMEANF(nLng,nLat,nT)	real
C	DFIXM2(nLng,nLat,nT)	real
C	WTP(NLOS,NL)		real
C	Ztmp(nLng,nLat,nT)	real
C
C FUNCTIONS/SUBROUTINES:
C	ArrR4Bad, ArrR4Zero, ArrR4getminmax, GridSphere2D, Timesmooth, Say
C PROCEDURE:
C MODIFICATION HISTORY:
C	NOV, 1995 B. Jackson (STEL,UCSD)
C-
	subroutine MkDMaptdn0n(LON,LAT,ITIM,GWTij,PW,FIX,NL,NLOS,NLOSP1,
     &		nLng,nLat,nT,DMAP,ANLIMIT,GSSIG,SIGB,NS,NBADSG,SIGMAP,
     &		ANMAP,AWTMAP,DMA,FIXM,WTSM,DMEANF,DFIXM2,WTP,Alng,ConstL,CONR,CONT,aNday,Ztmp)
     
	integer		LON(NLOSP1,NL),		! Projected Carr. rot. value of point on L.O.S.
     &			LAT(NLOSP1,NL),		! Heliographic lat. (in deg.) point on L.O.S.
     &                  ITIM(NLOSP1,NL),        ! Projected time values
     &			NBADSG(NL)		! Bad source velocity 0 - bad, 1 - good
	real		GWTij(NLOS,NL),		! m weights of each point along the L.O.S.
     &			FIX(NL),		! Ratio observed/model G for total L.O.S.
     &			GSSIG(NL),		! Deviation of source from mean in sigma
     &			DMAP(nLng,nLat,nT)	! Density map
 
C	Scratch arrays
	real		ANMAP(nLng,nLat,nT),	! Number of points per bin
     &			AWTMAP(nLng,nLat,nT),
     &			DMA(nLng,nLat,nT),
     &			FIXM(nLng,nLat,nT),
     &			WTSM(nLng,nLat,nT),
     &			DMEANF(nLng,nLat,nT),
     &			DFIXM2(nLng,nLat,nT),
     &			WTP(NLOS,NL),
     &			Ztmp(nLng,nLat,nT)
     
	character	cStr*150
	include		'MAPANGLES.H'
	nLngLat	= nLng*nLat
	nLngLatnT = nLng*nLat*nT
	nLOSnL	  = NLOS*NL
	Bad = BadR4()
	nLat1 = nLat - 1
	daLat2 = 90./nLat
	call ArrR4Bad(nLngLatnT,ANMAP)
	call ArrR4Bad(nLngLatnT,AWTMAP)
	call ArrR4Bad(nLngLatnT,DMA)
	call ArrR4Zero(nLngLatnT,FIXM)
	call ArrR4Bad(nLngLatnT,WTSM)
	call ArrR4Zero(nLngLatnT,DMEANF)
	call ArrR4Zero(nLngLatnT,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
		ITlast = 0
		do J=1,NLOS
		    LN = LON(J,I)			! map longitude
		    LT = LAT(J,I)			! map latitude
                    IT = ITIM(J,I)                      ! integer time
		    WTSj = abs(GWTij(J,I) )
C		    FWS = WTSj*(FIX(I)**(1.0/(2.0*PW)))
C		    if(FIXM(LN,LT,IT).eq.Bad) FIXM(LN,LT,IT) = 0.
		    FIXM(LN,LT,IT) = FIXM(LN,LT,IT)+WTSj*FIX(I)
		    WTP(J,I) = WTSj
		    if(WTSM(LN,LT,IT).eq.Bad) WTSM(LN,LT,IT) = 0.
		    WTSM(LN,LT,IT) = WTSM(LN,LT,IT)+WTP(J,I)
C	print *, LN, LT, IT, FIXM(LN,LT,IT), WTSj, FIX(I)
		    if (LN.ne.LNlast.or.LT.ne.LTlast.or.IT.ne.ITlast) then
			if(ANMAP(LN,LT,IT).eq.Bad) ANMAP(LN,LT,IT) = 0.
     			ANMAP(LN,LT,IT) = ANMAP(LN,LT,IT)+1.
		    end if
		    LNlast = LN
		    LTlast = LT
		    ITlast = IT
		end do
	    end if
	end do
	do N = 1, nT
	     do I=1,nLng
	     	 do J=1,nLat
		     if(ANMAP(I,J,N).eq.0.0) ANMAP(I,J,N) = badR4()
		 end do
	     end do
	     call ArrR4getminmax(nLngLat,ANMAP(1,1,N),amin,amax)
C	print *, 'In MKDMAPTDN N, amax', N, amax
	     if(amax.ne.Bad) then
		call GridSphere2D(ALng,nLng,nLat,1,ANMAP(1,1,N),CONR,3,0.0,0.0)
C
C  Now modify the latitude map count to approximate the cos latitude filter effect
C
		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,N).ne.Bad) then
			    ANMAP(I,J,N) = ANMAP(I,J,N)*sqonebycos
			end if
		    end do
		end do
C
	    end if
	end do
	call Timesmooth(nLng,nLat,nT,ANMAP,0,1.0*CONT/aNday,0.,Ztmp)
	
C	print *, nLng,nLat,nT

	do LN=1,nLng
	    do LT=1,nLat
		do IT=1,nT
C		    if(LN.eq.100.and.LT.eq.18) print *, LN,LT,IT,FIXM(LN,LT,IT),WTSM(LN,LT,IT)
		    if (FIXM(LN,LT,IT) .gt. 0. .and. WTSM(LN,LT,IT) .gt. 0.) then
			DMEANF(LN,LT,IT) = FIXM(LN,LT,IT)/WTSM(LN,LT,IT)
			DMA(LN,LT,IT) = DMEANF(LN,LT,IT)*DMAP(LN,LT,IT)
C	print *, LN,LT,IT,FIXM(LN,LT,IT),WTSM(LN,LT,IT),DMEANF(LN,LT,IT),DMAP(LN,LT,IT),DMA(LN,LT,IT)
		    else
			ANMAP(LN,LT,IT) = 0.
		    end if
		end do
	    end do
	end do
C  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
		IT = ITIM(J,I)			! map time
		if (DMA(LN,LT,IT) .ne. Bad) DFIXM2(LN,LT,IT)=DFIXM2(LN,LT,IT)
     &		    +WTP(J,I)*((FIX(I)-DMEANF(LN,LT,IT))**2)
	    end do
	    end if
	end do
C  calculation of the "total" Chi square distribution for all latitudes and longitudes
	CONVER  = 0
	NCONVER = 0
	do I=1,nLng
	    do J=1,nLat
		do K=1,nT
		    if (ANMAP(I,J,K) .lt. ANLIMIT) DMA(I,J,K) = Bad ! data not deconvolved
		    if (ANMAP(I,J,K).ne.Bad.and.WTSM(I,J,K).ne.0.0.and.WTSM(I,J,K).ne.Bad) then
		    	AWTMAP(I,J,K) = WTSM(I,J,K)
		    end if		    
		    if (DMA(I,J,K) .ne. Bad) then
			CONVER = CONVER+DFIXM2(I,J,K)
			NCONVER = NCONVER+1
		    end if
		    DMAP(I,J,K) = DMA(I,J,K)
C	if(DMAP(I,J,K).gt.1.0e7.and.DMAP(I,J,K).lt.1.0e22) print *, I, J, DMAP(I,J,K)
		end do
	    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,F9.4,A,F8.3,A)')
     &		'#D map convergence =',SIGMAP,', from last =',abs(1-SIGMAP/SIGOLD)*100,'%'
	call Say('MkDMapn','I','Info',cStr)
	return
	end