C+
C NAME:
C	MkGModeltdn
C PURPOSE:
C	Make model line of sight G-levels from a density map at a
C	given height. Following, a latex file containing the scintillation
C	index value m and the weighting function.
C	Line-of-sight integration
.  
C	This is the same g-level model-making program as used in ipsdtd just
C	in case the ipshdt program is to be used for g-levels.  However, the 
C	main ipshtdt program has never been tested in this mode of operation. 
C
C	Scintillation index \\
C	   $$ m^2 = \int dz W(z) $$ 
C	Weighting function \\
C	   $$ W(z) = (\Delta N_e(z))^2 \int dq 
C	          \sin^2(\frac{q^2\lambda z}{4\pi})
C	          \exp(-\frac{\theta_o^2 q^2 z^2}{2})
C	          q^{-3} $$
C CATEGORY:
C	Data processing
C CALLING SEQUENCE:
C	call MkGModeltd(XCbe,XCtbeg,XCtend,XCtpr,XLON,XLAT,RP,
C		WTS2,DFAC,PW,DEN1AU,NL,NLOS,NLOSP1,DMAP,
C		nLng,nLat,nT,RRS,GM2,GWTij,GMWTij,GMWTi)
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	PW		real		Density power
C	DEN1AU		real		Average density at 1 AU
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	DMAP(nLng,nLat,nT)	real	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	GM2(NL)		real		Model G-levels for a given source
C	GWTij(NLOS ,NL)	real		Value weights
C SCRATCH ARRAYS:
C	GMWTij(NLOS,NL) real
C	GMWTi(NL)	real
C FUNCTIONS/SUBROUTINES:
C	Get4Dval
C PROCEDURE:
C MODIFICATION HISTORY:
C	NOV, 1995 B. Jackson (STEL,UCSD)
C-
	subroutine MkGModeltdn(XCbe,XCtbeg,XCtend,XCtpr,XLON,XLAT,RP,
     &		WTS2,DFAC,PW,DEN1AU,NL,NLOS,NLOSP1,DMAP,
     &		nLng,nLat,nT,RRS,GM2,GWTij,GMWTij,GMWTi)
     
	real		XLON(NLOSP1,NL),	! Projected Carr. rot. value of point on L.O.S.  
     &			XCbe(2,nT),		! Beginning and ending Carr. variables
     &			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
     &			DMAP(nLng,nLat,nT),	! Density map
     &			GM2(NL),		! Model G-levels^2 for a given source
     &			GWTij(NLOS ,NL)
	real		GMWTij(NLOS,NL), 	 ! Scratch arrays
     &			GMWTi(NL)
     
	real*8		XCtpr(NLOSP1,NL)	! Projected time value on surface
	
	real*8		XCtbeg,XCtend		! Fixed 6/11/04 BVJ

	BADD = BadR4()
	R1AU = 1.
	NDENj = 0
	do I=1,NL
	    GMWTi(I) = 0
            do J=1,NLOS
                call Get3DTval(XCbe,XCtbeg,XCtend,nLng,nLat,nT,DMAP,1,
     &                  XLON(J,I),XLAT(J,I),XCtpr(J,I),DENj)
		if (DENj .eq. BADD) stop 'Bad density in MkGModeltdn'
		DENj = DENj*DFAC(J,I)*((RRS/RP(J,I))**2)
C	if(DENj.le.0.) print *, 'Bad density in MkGmodeltdn', J,I,XLON(J,I),XLAT(J,I),XCtpr(J,I),DENj,DFAC(J,I),RRS,RP(J,I),WTS2(J,I),PW
	if(DENj.le.0.) NDENj = NDENj + 1
	if(DENj.le.0.) DENj = 0.001
C		GWTij(J,I)  = (DENj**(2.0*PW))*WTS2(J,I)/(RP(J,I)**2) ! Asai's
		GWTij(J,I)  = (DENj**(2.0*PW))*WTS2(J,I)
C		GWTij(J,I)  = (DENj**PW)*WTS2(J,I)*cosd(XLAT(J,I))
		DENMj = DEN1AU*((R1AU/RP(J,I))**2)
		GMWTij(J,I) = (DENMj**(2.0*PW))*WTS2(J,I)
C		GMWTij(J,I) = (DENMj**PW)*WTS2(J,I)*cosd(XLAT(J,I))
		GMWTi(I)    = GMWTi(I) + GMWTij(J,I)
	    end do
	if(GMWTi(I).eq.0) print *, I, DENMj, PW, WTS2(1,I), WTS2(NLOS,I)
	end do
	if(NDENj.ne.0) print *, NDENj, ' L.O.S. density values in MkGmodeltdn are below zero'
	do I=1,NL
	    GM2(I) = 0.
C	Scintillation index line of sight integration for both the modeled
C	density weighting (GWTij) and the mean density weighting (GMWTij).
	    do J=1,NLOS
		GWTij(J,I)  = GWTij(J,I)/GMWTi(I)
		GM2(I) 	    = GM2(I) + GWTij(J,I)
	    end do
	end do
	do I=1,NL
	    do J=1,NLOS
		GWTij(J,I)  = GWTij(J,I)/GM2(I)
	    end do
	end do
	return
	end