C+
C NAME:
C	MkDHModeltd
C PURPOSE:
C	Make model line of sight brightnesses from a density map at a
C	given height. 
Both model base brightnesses and model variation
C	brightnesses are output. Many different weighting schemes were tried,
C	and the weights currently used that weight simply as the brightness
C	may not be optimal for reconstructions of the heliospheric model at 
C	Earth since the heliosphere is less bright at 90 degrees, where 
C	Helios and SMEI might look.
C
C CATEGORY:
C	Data processing
C CALLING SEQUENCE:
C	call MkDHModeltd(XCbe,XCtbeg,XCtend,XCtpr,XLON,XLAT,RP,
C		WTS2,DFAC,PWR,NL,NLOS,DIST,NLOSP1,DMAP, DBASE,OBSB,
C		nLng,nLat,nT,RRS,GM,GWTij)
C INPUTS:
C       XCbe(2,nT)      real    	Beginning and ending Carrington variables
C	XCtbeg		real		Beginning time interval
C	XCtend		real		Ending time variable
C	XCtpr(NLOSP1,NL)real*8		Projected time
C	XLON(NLOSP1,NL)	integer		Projected Carr. rot. value of point on L.O.S.
C	XLAT(NLOSP1,NL)	integer		Heliographic lat. (in deg.) point on L.O.S.
C	RP  (NLOSP1,NL)	real		Distance above Sun of point on L.O.S.
C	WTS2(NLOS  ,NL)	real		Weights of each point along the L.O.S.
C	DFAC(NLOSP1,NL)	real		Density factors for each L.O.S. point
C	PWR		real		Density falloff with power
C	NL		integer		Number of G-level data points
C	NLOS		integer		Number of L.O.S. distance segments
C	DIST(NL)	real		Helios distance from the Sun
C	NLOSP1		integer		Number of L.O.S. distance segments + 1
C	DMAP (nLng,nLat,nT)	real	Density map  (Total = Base plus density)
C	DBASE(nLng,nLat,nT)	real	Base density map
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	dRR		real		Interval between height maps
C	RRS		real		Height of deconvolution surfaces
C OUTPUTS:
C	OBSB  (NL)	real		Observed base brightness (calculated here from base density)
C	GM(NL)		real		Model G-levels for a given source
C	GWTij(NLOS ,NL)	real		Value weights
C FUNCTIONS/SUBROUTINES:
C	Get3DTval
C PROCEDURE:
C MODIFICATION HISTORY:
C	NOV, 1995 B. Jackson (STEL,UCSD)
C-
	subroutine MkDHModeltd(XCbe,XCtbeg,XCtend,XCtpr,XLON,XLAT,RP,
     &		WTS2,DFAC,PWR,NL,NLOS,DIST,NLOSP1,DMAP, DBASE,OBSB,
     &		nLng,nLat,nT,RRS,GM,GWTij)
	real		XCbe  (2,nT),		! Beginning and ending Carr. variables
     &			XLON  (NLOSP1,NL),	! Projected Carr. rot. value of point on L.O.S.  
     &			XLAT  (NLOSP1,NL),	! Heliographic lat. (in deg.) point on L.O.S.
     &			RP    (NLOSP1,NL),	! Distance above Sun of point on L.O.S.
     &			WTS2  (NLOS  ,NL),	! Weights of each point along the L.O.S.
     &			DFAC  (NLOSP1,NL),	! Density factors for each L.O.S. point
     &			DIST  (NL),		! Helios distance from the Sun
     &			DMAP  (nLng,nLat,nT),	! Density map  (Total = Base + Density)
     &			DBASE (nLng,nLat,nT),	! Base density
     &			OBSB   (NL),		! Observed base brightness (calculated here from base density)
     &			GM    (NL),		! Model G-levels for a given source
     &			GWTij (NLOS  ,NL)
     
	real*8		XCtpr(NLOSP1,NL)	! Projected time value on surface
	
	real*8		XCtbeg,XCtend		! Fixed 6/11/04 BVJ

	BADD = BadR4()
	CONSTWT = 5.0/( RRS**2)
	R1AU = 1.
	NDENj = 0
	do I=1,NL
	    GWTj = 0.0
	    AM   = 0.0
	    OB	 = 0.0
	    WTj  = 0.0
	    if(WTS2(1,I).eq.0.0) print *, ' WTS2(1,I) = 0.0, I =', I
            do J=1,NLOS
                call Get3DTval(XCbe,XCtbeg,XCtend,nLng,nLat,nT,DMAP,1,
     &                  XLON(J,I),XLAT(J,I),XCtpr(J,I),DENjB)
		if (DENjB .eq. BADD) stop 'Bad density in MkDHModeltd'
		DENj = DENjB*DFAC(J,I)*((RRS/RP(J,I))**PWR)
                call Get3DTval(XCbe,XCtbeg,XCtend,nLng,nLat,nT,DBASE,1,			! Allows a base different from a constant
     &                  XLON(J,I),XLAT(J,I),XCtpr(J,I),DENBB)
		if (DENBB .eq. BADD) stop 'Bad Base density found in MkDHModeltd'
		if (DENBB .le. 0.0)  stop 'Negative Base density found in MkDHModeltd'
		DENB = DENBB*((RRS/RP(J,I))**PWR)
		if(I.eq.1.and.DENj.le.0.) print *, 'Density below zero in MkDHmodeltd', J,I
C	J,I,XLON(J,I),XLAT(J,I),XCtpr(J,I),DENj,DFAC(J,I),RRS,RP(J,I),WTS2(J,I)
		if(DENj.le.0.) NDENj = NDENj + 1
		if(DENj.le.0.) DENj = 0.001
C The way it was:
C		GWTij(J,I)  = WTS2(J,I)
C		GWTij(J,I)  = GWTij(J,I)*cosd(XLAT(J,I))
C		DTW = DENj*WTS2(J,I)
C		DTW = DTW*cosd(XLAT(J,I))
C		GWTj = GWTj + GWTij(J,I)
C		AM  = AM + DTW
C New

C*******************************************************************
C		DENj = 20.0		! Put in for 20 Jan 2005 shock reconstruction
C*******************************************************************

		DTW  = DENj*WTS2(J,I)
		AM   = AM + DTW

C*******************************************************************
C		print *, I,J, DTW, AM	! Put in for 20 Jan 2005 shock reconstruction
C*******************************************************************		
		
		OB   = OB + DENB*WTS2(J,I)
		
C		GWTij(J,I) = DENjB*DENj*WTS2(J,I)
C		GWTij(J,I) = DENjB
C		GWTij(J,I) = DENjB*DENj

		GWTij(J,I) = DTW
		WTj 	   = WTj + WTS2(J,I)
C	if (I.eq.7000.and.J.eq.1) print *, 'I, J, DENB, RRS, RP(J,I), PWR', I, J, DENB, RRS, RP(J,I), PWR
C	if (I.eq.7000.and.J.le.2) print *, 'I, J, DENj, DENjB, DFAC(J,I), WTS2(J,I), AM, OB',
C     &				 	    I, J, DENj, DENjB, DFAC(J,I), WTS2(J,I), AM, OB
C	if (I.eq.7000.and.J.ge.NLOS-1) print *, 'I, J, DENj, DENjB, DFAC(J,I), WTS2(J,I), AM, OB',
C     &				 	    	 I, J, DENj, DENjB, DFAC(J,I), WTS2(J,I), AM, OB
	    end do
	    GM(i)   = AM
	    OBSB(i) = OB
	    do j=1,NLOS
C		GWTij(j,i) = GWTij(j,i)/GWTj 		 ! makes every line of sight the same weight (the way it was)
C		GWTij(j,i) = GWTij(j,i)/WTj/CONSTWT	! make every line of sight weight according to base density 
							! and relative line of sight weight (session 4)(DENjB*WTS2(J,I)/Wtj)
C		GWTij(j,i) = GWTij(j,i)/CONSTWT		! base density times the l.o.s. weight at that position
							! (session 5) (DENjB*WTS2(J,I))
C		GWTij(j,i) = GWTij(j,i)/(CONSTWT*1000.)	! base density times the l.o.s. weight including density 
							! at that position (session 10) (DENjB*DENj*WTS2(J,I))
C		GWTij(j,i) = GWTij(j,i)*AM/(CONSTWT*1000.) ! base density times the total l.o.s. weight including density 
							! in that direction (session 11) (DENjB*AM)
C		GWTij(j,i) = GWTij(j,i)/(CONSTWT*1000.) ! base density times the deconvolved density at that position 
							! (session 12) (DENjB*DENj)
		GWTij(j,i) = GWTij(j,i)			! the deconvolved density times the L.O.S. weight at that position 
							! (session 13) (DTW) (because the base density ratio fix is relative to 1)
	    end do
	end do
        if(NDENj.ne.0) print *, NDENj,' L.O.S. density values in MkDHmodeltd are below zero'
	return
	end