C+
C NAME:
C	MkFModeltd
C PURPOSE:
C	Make model line of sight Faraday Rotation values from a density map and a magnetic field determined at a
C	given height. Following is shown the way this is done at a given frequency and height (A LaTex file 
C	describing same).
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 MkFModeltd(XCbe,XCtbeg,XCtend,XCtpr,XLON,XLONP,XCEG,ELONG,DIST,XLAT,XLATP,RP,IYRS,DOYS8,
C		DFAC,VFAC,DEN1AU,dLOSF,FRCONST,NL,NLMAXF,NLOS,NLOSP1,DMAP,FMAP,
C		nLng,nLat,nT,RRS,GM2,GWTij,GMWTij,SIGNPM,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)		real	Projected Carr. rot. value of point on L.O.S.
C	XLONP(NLOSP1,NL)	real	Carr. rot. value of L.O.S. point
C	XCEG (NL),		real	Carr. rot. value of Earth
C	ELONG(NL)		real	Elongation of L.O.S.
C	DIST (NL)		real	Distance of Earth at the time of the observation 
C	XLAT(NLOSP1,NL)		real	Projected Heliographic lat. (in deg.) point on L.O.S.
C	XLATP(NLOSP1,NL)	real	Heliographic lat. (in deg.) of L.O.S. point
C	RP  (NLOSP1,NL)		real	Distance above Sun of point on L.O.S.
C	IYRS (NL),		integer	Year source was observed
C	DOYS8(NL),		real*8	Day of year source was observed
C	DFAC(NLOSP1,NL)		real	Density factors for each L.O.S. point
C	DEN1AU			real	Average density at 1 AU
C	dLOSF			real	Distance of a line of sight element in AU
C	FRCONST			real	Faraday Rotation constant amount according to frequency
C	NL			integer	Number of FR data points
C	NLMAXF			integer Total number of possible FR 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	FMAP(nLng,nLat,nT,3)	real	Magnetic field 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	FM2(NL,3)		real	Model Faraday rotation values for a given source
C	FWTij(NLOS,NL,3)	real	Value weights (plus and minus)
C	SIGNPM(NLOS,NL)		real	Sign of the positive model (+1.0 or -1.0)
C SCRATCH ARRAYS:
C	FMWTi(NL)		real	An absolute value
C FUNCTIONS/SUBROUTINES:
C	Get3DTval
C PROCEDURE:
C MODIFICATION HISTORY:
C	NOV, 1995 B. Jackson (STEL,UCSD)
C-
	subroutine MkfModeltd(XCbe,XCtbeg,XCtend,XCtpr,XLON,XLONP,XCEG,ELONG,DIST,XLAT,XLATP,RP,IYRS,DOYS8,
     &		DFAC,VFAC,DEN1AU,dLOSF,FRCONST,NL,NLMAXF,NLOS,NLOSP1,DMAP,VMAP,FMAP,
     &		nLng,nLat,nT,nTmax,RRS,FM2,FWTij,SIGNPM,FMWTi)

	real		XLON (NLOSP1,NL),		! Projected Carr. rot. value of point on L.O.S.
     &			XLONP(NLOSP1,NL),		! Carr. rot. value of L.O.S. point
     &			XCEG (NL),			! Carr. rot. value of Earth
     &			ELONG(NL),			! Elongation of L.O.S.
     &			DIST (NL),			! Distance of Earth at the time of the observation 
     &			XCbe (2,nT),			! Beginning and ending Carr. variables
     &			XLAT (NLOSP1,NL),		! Projected heliographic lat. (in deg.) point on L.O.S.
     &			XLATP(NLOSP1,NL),		! Heliographic lat. (in deg.) of L.O.S. point
     &			RP   (NLOSP1,NL),		! Distance above Sun of point on L.O.S.
     &			DFAC (NLOSP1,NL),		! Density factors for each L.O.S. point
     &			VFAC (NLOSP1,NL),		! Velocity factors for each L.O.S. point
     &			DMAP (nLng,nLat,nT),		! Density map
     &			VMAP (nLng,nLat,nT),		! Velocity map
     &			FMAP (nLng,nLat,nTmax,3),	! Magnetic field map (3 component, longitude, latitude, R)
     &			FM2  (NLMAXF,3),		! Model Faraday rotation values for a given source
     &			FWTij(NLOS,NLMAXF,3),		! Magnetic value weights (plus and minus)
     &			SIGNPM(NLOS,NLMAXF,3)		! Sign of the positive model

	real*8		XCtpr(NLOSP1,NL)		! Projected time value on surface
	
	real*8		XCtbeg,XCtend,DOYS8(NL),Doy8

	real		FMWTi(NL,3)			! Scratch arrays

        real 		FIELD(3)
	
	integer		IYRS(NL)

	include		'sun.h'		! Contains executable

	BADD = BadR4()

	R1AU   = 1.
	NDENj  = 0
	FMINWT = 0.01
	FRCON  = FRCONST*dLOSF
	dLOSF2 = dLOSF/2.
	OMEGA  = 100.*Sun__Omega*Sun__AU
	ROTC   = 360.
	Icnt = 0
	do I=1,NL
	    do K=1,3
	      FM2(I,K)   = 0.
	    end do
	    FMWTi(I,1) = 0.
	    FMWTi(I,2) = 0.
	    FMWTi(I,3) = 0.
            do J=1,NLOS
	        PD = dLOSF*float(J) - dLOSF2
                call Get3DTval(XCbe,XCtbeg,XCtend,nLng,nLat,nT,DMAP,1,
     &                  XLON(J,I),XLAT(J,I),XCtpr(J,I),DENj)
     		if (DENj .eq. BADD) print *, I,J,XCtbeg,XCtend,nLng,nLat,nT,XLON(J,I),XLAT(J,I),XCtpr(J,I),DENj
		if (DENj .eq. BADD) stop 'Bad density in MkFModeltd'
		DENj = DENj*DFAC(J,I)*((RRS/RP(J,I))**2)		! Density at L.O.S. location
	        if(DENj.le.0.) NDENj = NDENj + 1
		if(DENj.le.0.) DENj = 0.001

                call Get3DTval(XCbe,XCtbeg,XCtend,nLng,nLat,nT,VMAP,1,
     &                  XLON(J,I),XLAT(J,I),XCtpr(J,I),Vj)
		if (Vj .eq. BADD) stop 'Bad velocity in MkFModeltd'
     		Vj = Vj*VFAC(J,I)					! Velocity at L.O.S. location

		do K=1,3
		  call Get3DTval(XCbe,XCtbeg,XCtend,nLng,nLat,nT,FMAP(1,1,1,K),1,
     &              XLON(J,I),XLAT(J,I),XCtpr(J,I),FIELD(K))
		  if (FIELD(K) .eq. BADD) stop 'Bad magnetic field in MkFModeltd'
		end do

		FIELDLN1 = ((RRS/RP(J,I))/VFAC(J,I))*
     &				(OMEGA*RRS*cosd(XLATP(J,I))/Vj)		! Mean Longitudinal heliographic field on L.O.S. point
		FIELDLT1 = (RRS/RP(J,I))/VFAC(J,I)			! Mean Latitudinal heliographic field on L.O.S. point
		FIELDR1  = (((RRS/RP(J,I))/VFAC(J,I))**2)		! Mean Radial field on L.O.S. point

		FIELDLN = FIELD(1)*((RRS/RP(J,I))/VFAC(J,I))*
     &				(OMEGA*RRS*cosd(XLATP(J,I))/Vj)		! Longitudinal heliographic field on L.O.S. point
		FIELDLT = FIELD(2)*(RRS/RP(J,I))/VFAC(J,I)		! Latitudinal heliographic field on L.O.S. point
		FIELDR  = FIELD(3)*(((RRS/RP(J,I))/VFAC(J,I))**2)	! Radial field on L.O.S. point

		Ierr = 0
		if(Ierr.eq.0) go to 222
		if(I.eq.245.and.J.eq.20) 
     &		print *, 'mkfmodeltd I,J,DIST(I),ELONG(I),XCEG(I),PD',I,J,
     &								DIST(I),ELONG(I),XCEG(I),PD
     		if(I.eq.800.and.J.eq.20) 
     &		print *, 'mkfmodeltd I,J,DIST(I),ELONG(I),XCEG(I),PD',I,J,
     &								DIST(I),ELONG(I),XCEG(I),PD

		if(I.eq.245.and.J.eq.1) 
     &		print *, 'mkfmodeltd I,J,OMEGA,XLONP(J,I),XLATP(J,I),RP(J,I),VFAC(J,I),DFAC(J,I),Vj,',I,J,
     &								OMEGA,XLONP(J,I),XLATP(J,I),RP(J,I),VFAC(J,I),DFAC(J,I),Vj

		if(I.eq.245.and.J.eq.20) 
     &		print *, 'mkfmodeltd I,J,OMEGA,XLONP(J,I),XLATP(J,I),RP(J,I),VFAC(J,I),DFAC(J,I),Vj,',I,J,
     &								OMEGA,XLONP(J,I),XLATP(J,I),RP(J,I),VFAC(J,I),DFAC(J,I),Vj
     
     		if(I.eq.800.and.J.eq.20) 
     &		print *, 'mkfmodeltd I,J,OMEGA,XLONP(J,I),XLATP(J,I),RP(J,I),VFAC(J,I),DFAC(J,I),Vj,',I,J,
     &								OMEGA,XLONP(J,I),XLATP(J,I),RP(J,I),VFAC(J,I),DFAC(J,I),Vj
		
		if(I.eq.245.and.J.eq.20) 
     &		print *, 'mkfmodeltd I,J,FIELD(1),FIELD(2),FIELD(3),FIELDLN,FIELDLT,FIELDR,DENj',I,J,
     &					    FIELD(1),FIELD(2),FIELD(3),FIELDLN,FIELDLT,FIELDR,DENj
     
		if(I.eq.800.and.J.eq.20) 
     &		print *, 'mkfmodeltd I,J,FIELD(1),FIELD(2),FIELD(3),FIELDLN,FIELDLT,FIELDR,DENj',I,J,
     &					    FIELD(1),FIELD(2),FIELD(3),FIELDLN,FIELDLT,FIELDR,DENj
 222  continue
C		For this system the magnetic field positive is defined as + longitude, + latitude, and + (outward) radial heliographic

		iYr = IYRS(I)
		Doy8 = DOYS8(I)
		PLN = (XCEG(I) - XLONP(J,I))*ROTC			! Heliographic angular distance to the L.O.S. point
		PLN = PLN + EARTH8(iYr,Doy8)				! Heliographic longitude of the L.O.S. point
		PLT = XLATP(J,I)					! Heliographic latitude of L.O.S. point 

		call ECLIPTIC_HELIOGRAPHIC8(1,iYr,Doy8,PLN,PLT)		! Convert heliographic to ecliptic coordinates

		EP    = EARTH8(-iYr,Doy8) - PLN				! Minus the ecliptic angular distance to the L.O.S. point
		PP    = RP(J,I)*cosd(PLT)				! Projected ecliptic equatorial radial distance of L.O.S. point		
		SQEP  = SQRT(PP*PP + DIST(I)*DIST(I) - 2.*PP*DIST(I)*cosd(EP))		! Law of cosines
		COSLN =	DIST(I)*sind(EP)/SQEP				! Direction cosine along L.O.S. in longitude at L.O.S. point
		XLTD  = -RP(J,I)*sind(PLT)			! Distance of L.O.S. point above the equator
		COSLT = XLTD/PD						! Direction cosine along L.O.S. in latitude at L.O.S. point
		COSR  =	(DIST(I)*cosd(ELONG(I)) - PD)/RP(J,I)		! Direction cosine along L.O.S. in radial at L.O.S. point
		
C		if(I.eq.245.and.J.eq.20) print *, 'mkfmodeltd I,J,EP,PP,PLN,PLT',I,J,EP,PP,PLN,PLT
C		if(I.eq.245.and.J.eq.20) print *, 'mkfmodeltd I,J,COSLN,COSLT,COSR',I,J,COSLN,COSLT,COSR
C		if(I.eq.800.and.J.eq.20) print *, 'mkfmodeltd I,J,EP,PP,PLN,PLT',I,J,EP,PP,PLN,PLT
C		if(I.eq.800.and.J.eq.20) print *, 'mkfmodeltd I,J,COSLN,COSLT,COSR',I,J,COSLN,COSLT,COSR
				
		BLN1 = COSLN*FIELDLN1*DENj				! Mean B*Ne parallel from longitude at L.O.S. point
		BLT1 = COSLT*FIELDLT1*DENj				! Mean B*Ne parallel from latitude at L.O.S. point
		BR1  = COSR*FIELDR1*DENj				! Mean B*Ne parallel from radial at L.O.S. point

		BLN = COSLN*FIELDLN*DENj				! B*Ne parallel from longitude at L.O.S. point
		BLT = COSLT*FIELDLT*DENj				! B*Ne parallel from latitude at L.O.S. point
		BR  = COSR*FIELDR*DENj					! B*Ne parallel from radial at L.O.S. point

C		if(I.eq.245.and.J.eq.20) print *, 'mkfmodeltd I,J,BLN,BLT,BR,BLN1,BLT1,BR1',I,J,BLN,BLT,BR,BLN1,BLT1,BR1
		
C		if(I.eq.800.and.J.eq.20) print *, 'mkfmodeltd I,J,BLN,BLT,BR,BLN1,BLT1,BR1',I,J,BLN,BLT,BR,BLN1,BLT1,BR1
C
		FM2(I,1) = FM2(I,1) + BLN*FRCON
		FM2(I,2) = FM2(I,2) + BLT*FRCON
		FM2(I,3) = FM2(I,3) + BR*FRCON
		
C		FWTij(J,I,1)  = DENj*FIELDLN*abs(COSLN)			! Ne*B FR wt from longitude in degrees (+- according to field)
C		FWTij(J,I,2)  = DENj*FIELDLT*abs(COSLT)			! Ne*B FR wt from latitude in degrees (+- according to field)
C		FWTij(J,I,3)  = DENj*FIELDR*abs(COSR)			! Ne*B FR wt from radial in degrees (+- according to field)

C		FWTij(J,I,1)  = 2.0**(abs(BLN1))			! Ne*B FR wt from longitude in degrees (+- according to cos*field)
C		FWTij(J,I,2)  = 2.0**(abs(BLT1))			! Ne*B FR wt from latitude in degrees (+- according to cos*field)
C		FWTij(J,I,3)  = 2.0**(abs(BR1))				! Ne*B FR wt from radial in degrees (+- according to cos*field)
		
C		FWTij(J,I,1)  = abs(BLN1)				! Ne*B FR wt from longitude in degrees (+- according to cos*field)
C		FWTij(J,I,2)  = abs(BLT1)				! Ne*B FR wt from latitude in degrees (+- according to cos*field)
C		FWTij(J,I,3)  = abs(BR1)				! Ne*B FR wt from radial in degrees (+- according to cos*field)	

C		FWTij(J,I,1)  = BLN1	(was this 9/14/08)		! Ne*B FR wt from longitude in degrees (+- according to cos*field)
C		FWTij(J,I,2)  = BLT1					! Ne*B FR wt from latitude in degrees (+- according to cos*field)
C		FWTij(J,I,3)  = BR1					! Ne*B FR wt from radial in degrees (+- according to cos*field)	

		FWTij(J,I,1)  = abs(BLN*FRCON)				! Ne*B FR wt from longitude in degrees (+- according to cos*field)
		FWTij(J,I,2)  = abs(BLT*FRCON)				! Ne*B FR wt from latitude in degrees (+- according to cos*field)
		FWTij(J,I,3)  = abs(BR*FRCON)				! Ne*B FR wt from radial in degrees (+- according to cos*field)	

C		FMWTi(I,1) = FMWTi(I,1) + FWTij(J,I,1)
C		FMWTi(I,2) = FMWTi(I,2) + FWTij(J,I,2)
C		FMWTi(I,3) = FMWTi(I,3) + FWTij(J,I,2)

	    	SIGNPM(J,I,1) = BLN1/abs(BLN1)
	    	SIGNPM(J,I,2) = BLT1/abs(BLT1)
	    	SIGNPM(J,I,3) = BR1/abs(BR1)
		
		do K=1,3
C		  FMWTi(I)	= FMWTi(I) + abs(FWTij(J,I,K))		! Total absolute wt applied for each line of sight
		  FMWTi(I,K)	= FMWTi(I,K) + FWTij(J,I,K)		! Total wt applied for each line of sight (was tried
C		  FMWTi(I,1)	= FMWTi(I,1) + FWTij(J,I,K)		! Total wt applied for each line of sight
		end do
	    end do

	    do K=1,3
	      if(FMWTi(I,K).eq.0.) then					! Will never happen with an abs function above.
C	        print *, 'Zero L.O.S. model wt',I,K, DENj, FMWTi(I,K)
		Icnt = Icnt + 1
	        FMWTi(I,K) = FMINWT
	      end if
	    end do
	    
	end do
	if(Icnt.ne.0) print *, 'Zero L.O.S. model wt for',Icnt,' LOS: all set to small value'
	if(NDENj.ne.0) print *, NDENj, ' L.O.S. density values in MkGmodeltdn are below zero'

C	Normalized Faraday rotation line of sight weighting (FWTij). This is the mean FR component wt attributed 
C	to each L.O.S. point and this might not be optimal, or correct)

	do I=1,NL
	  do J=1,NLOS
	    do K=1,3
C	      FWTij(J,I,K) = FWTij(J,I,K)/FMWTi(I,K)			! Normalize each line of sight weight to one.
	    end do
C		if(I.eq.245.and.J.eq.20) print *, 
C     &		'mkfmodeltd I,J,FWTij(J,I,1),FWTij(J,I,2),FWTij(J,I,3)',I,J,FWTij(J,I,1),FWTij(J,I,2),FWTij(J,I,3)
C		if(I.eq.800.and.J.eq.20) print *, 
C     &		'mkfmodeltd I,J,FWTij(J,I,1),FWTij(J,I,2),FWTij(J,I,3)',I,J,FWTij(J,I,1),FWTij(J,I,2),FWTij(J,I,3)
	  end do
	end do

	return
	end