C+
C NAME:
C	READ_SMEI8
C PURPOSE:
C	Read Helios data for synoptic map
C CALLING SEQUENCE:
C	call READ_SMEI8(cWildCard,nCar,JDCar,XCbeg,ICnr,iYr,MJDref,
C    &		XElow,XEhigh,iEorW,XRlim,
C    &		iSMEI,cSMEI,iDay,nDev,dSig,
C    &		NLMAX,NS,NL,NLOC,
C    &		cgbfiles,SRCB,DIST,XC,XL,ZZ,XE,XLS,XCE,XPE,XLE,IYRS,DOYS8)
C INPUTS:
C	cWildCard
C	nCar
C	JDCar
C	XCbeg
C	ICnr,
C	iYr
C	XElow
C	XEhigh
C	iEorW
C	XRlim
C	iSMEI
C	cSMEI
C	iDay
C	nDev
C	dSig,
C	NS			! starting N value
C
C OUTPUTS:
C	NL
C	iSMEI
C	cSMEI
C	iDay
C	nDev
C	dSig,
C	NLMAX
C	NS			! starting N value
C	NL
C	SRCB
C	DIST
C	XC
C	XL
C	ZZ
C	XE
C	XLS
C	XCE
C	XPE
C	XLE
C	IYRS
C	DOYS8
C FUNCTIONS/SUBROUTINES:
C	CWRITE, ECLIPTIC_HELIOGRAPHIC, JULIAN8, 
C	SC_ECLIP, THOM_PNT_SRCE, T_FILTER, XMAP_SC_POS, XMAP_OBS_POS
C EXTERNAL:
C	external	EARTH8
C INCLUDE:
C	include		'sun.h'
C PROCEDURE:
C MODIFICATION HISTORY:
C
C	11/14/03  B. Jackson (UCSD): Provided to read SMEI data
C-
	subroutine READ_SMEI8(cWildCard,nCar,JDCar,XCbeg,ICnr,iYr,MJDref,
     &		Radius, XElow,XEhigh,iEorW,XRlim,
     &		iSMEI,cSMEI,iDay,nDev,dSig,
     &		ICIC,RA,DEC,ORBITV,ZSMEI,ISS,
     &		NORBITST,MX,NLOC,NLMAX,NS,NL,
     &		cgbfiles,SRCB,DIST,XC,XL,ZZ,XE,XLS,XCE,XPE,XLE,IYRS,DOYS8)
     
	include		'sun.h'
	real*8		JDCar(nCar)
	character	      cWildCard*(*),
     &			cSMEI(iSMEI)*(*),
     &			cgbfiles(*)*32,
     &			cgbfilesL*32,
     &			SRCB(*)*34,SRCBS*34
     
     	character	ctxt*4	/'.txt'/
     
	real		DIST(*),		! SMEI Distance from the Sun
     &			XC(*),			! Modified Carrington variable Point P
     &			XL(*),			! Heliographic latitude Point P (deg)
     &			ZZ(*),			! Intensity (S10)
     &			XE(*)			! Topocentric elongation Point P (deg)
						! (+=East;-=West)
	real		XLS(*),			! Topocentric ecliptic longitude Sun
     &			XCE(*),			! Modified Carrington variable of SMEI
     &			XPE(*),			! Topocentric ecliptic lng(P)-lng(Sun)
     &			XLE(*)			! Topocentric ecliptic lat(P)

	integer		IYRS(*)

	real*8		DOYS8(*)
C------
C These arrays are used to store data read from the SMEI data files

C	parameter	(NLOC = 990)		! Total number of RA and Dec locations possible (May, 2003)
C	parameter	(NLOC = 1108)		! Total number of RA and Dec locations possible (Oct., 2003)
C	parameter	(NLOC = 5000)		! Total number of RA and Dec locations possible per orbit
	
	integer		ICIC(*)

	real		RA(*),
     &			DEC(*),
     &			ORBITV(*),
     &			ZSMEI(*)		! Intensities (S10)
     
	integer		ISS(*)
     
	external	EARTH8
	real*8		JEPOCH,dTDoy,Doy8,xTT,MJDref
	real*8		dLngSun, dLatSun, dDisSun
	character	cName*100
	logical		bDataAvailable	/.FALSE./,
     &			bMore
	save		XElowOld , XEhighOld,
     &			iEorWOld , XRlimOld,
     &			iDayOld  , nDevOld , dSigOld
C-------
C Check whether new data have to be read.
	if (bDataAvailable			.and.
     &		XCbegOld	.eq. XCbeg	.and.
     &		iCnrOld	.eq. iCnr	.and.
     &		XElowOld	.eq. XElow	.and.
     &		XEhighOld	.eq. XEhigh	.and.
     &		iEorWOld	.eq. iEorW	.and.
     &		XRlimOld	.eq. XRlim	.and.
     &		iDayOld	.eq. iDay	.and.
     &		nDevOld	.eq. nDev	.and.
     &		dSigOld	.eq. dSig) return
C-------
C Save all parameters involved int the data selection.
C
	XCbegOld	= XCbeg
	iCnrOld	= iCnr
	XElowOld	= XElow
	XEhighOld	= XEhigh
	iEorWOld	= iEorW
	XRlimOld	= XRlim
	iDayOld	= iDay
	nDevOld	= nDev
	dSigOld	= dSig
C-------
C Get relative heliocentric and topocentric ecliptic coordinates of point P
C
C	TCONST	= 101.534/(60.*24.)
	iVel	= 400
	bMore	= .TRUE.
	IC	= XCbeg
	XClast	= XCbeg+ICnr-1
	print *, 'IC, XClast, bMore, MX =', IC, XClast, bMore, MX

	nGood1 = 0
	NLLL = 0
	NLLLo = 0
	NLLLoo = 0
	NL = NS
	do M=1,MX
	    if(bMore) then
C	print *, bMore, NL, NLMax
C*******************
C	The below was put in for the 20 January 2005 shock reconstruction
C
C	ZZZZE   = 0.0		! done to integrate brightness at 90 deg.
C	AZZZZE  = 0.0
C	rEEloP  = 0.0
C	arEElop = 0.0
C	ZZZZW   = 0.0
C	AZZZZW  = 0.0
C	rWEloP  = 0.0
C	arWElop = 0.0 
C	Ddoy8E  = 0.0
C	Adoy8E  = 0.0
C	Ddoy8W  = 0.0
C	Adoy8W  = 0.0
C*******************

		    cName = cWildCard
C		    print *, cName
		    if(NLLLoo.eq.0) iYrrr = iYr
		    call SMEIRead(cName,cgbfilesL,NORBITST,M,MX,NLOC,LOC,iYr,dTDoy,ISSD,ICIC,RA,DEC,ORBITV,ISS,ZSMEI,nGood1) ! Read a SMEI file of data in the interval
		
C		    do I=1,10
C		    print *, 'iYr, RA, DEC, ISS, ZSMEI', iYr, RA(I), DEC(I), ISS(I), ZSMEI(I)
C		    end do
 
		    do L=1,LOC
		      if (NL .lt. NLmax) then
			write(SRCBS(1:1),  '(I1)')   ICIC(L)
		    	write(SRCBS(2:8),  '(F7.2)') RA (L)
	    		write(SRCBS(9:15), '(F7.2)') DEC(L)
	    		write(SRCBS(16:22),'(F7.0)') Orbitv(L)
	    		write(SRCBS(23:26),'(I4.4)') ISS(L)					
			write(SRCBS(27:34),'(F8.2)') ZSMEI(L)
			
			write(cgbfilesL(29:32),'(A)') ctxt

C			print *, L,SRCBS

			Doy8 = dTDoy + dble(float(ISS(L)+ISSD)/86400.)
			
C			if(Doy8.gt.190.) print *, L,ISS(L),Doy8

			if(iyr.lt.2003.or.iyr.gt.2015.and.NLLLoo.le.10) then

			  print *, '******* Error ******** year < 2003 or > 2015 (iyr = ',iyr,')'
			  NLLLoo = NLLLoo + 1

			end if

C			if(NL.lt.10) print *, 'L, M, Doy8, ISS(L), ZSMEI(L)', L, M, Doy8, ISS(L), ZSMEI(L)

			call SunNewcomb8(0,iYr,Doy8,dLngSun,dLatSun,dDisSun)

						! Geocentric eclip. coord(Sun)
			rLng = sngl(dLngSun)
			rLngSun = rLng
			rLat = sngl(dLatSun)				! Geocentric ecliptic lat(Sun)
			call ECLIPTIC_EQUATOR8(0,iYr,Doy8,rLng,rLat)

			dRA = rLng				! RA (Sun)
			rLng  = RA(L)
			rLat  = DEC(L)
			dRA = rLng-dRA				! RA(P)-RA(Sun)
			if (abs(dRA) .gt. 180.) dRA = dRA-sign(1.,dRA)*360.
			call ECLIPTIC_EQUATOR8(1,iYr,Doy8,rLng,rLat)
						
								! Geocentric eclip. coord(P)
			rLng  = rLng-rLngSun
			rLngP = rLng
			rLatP = rLat
			XCsc = XMAP_SC_POS8(EARTH8,iYr,Doy8,-nCar,JDCar)

			rP = max(0.25, cosd(rLng)*cosd(rLat) )	! Geocentric distance P (at least 0.25 AU)
			call POINT_ON_LOS(rLng,rLat,rP,rEloP,EloSun,i)

						! Heliocentric coord(P) rel. Sun
c			print *, ELO,rEloP, ELO-rEloP
			rP    = rP*dDisSun			! Dist. Sun-P (AU)
			rEloP = i*rEloP				! rEloP>0 : East of Sun
								!      <0 : West of Sun					
			rLng = rLngSun+180.+rLng		! Heliocentric eclip. lng(P)
			call ECLIPTIC_HELIOGRAPHIC8(0,iYr,Doy8,rLng,rLat)

								! Heliographic coord(P)
			xLat = rLat				! Heliographic lat(P)
			XCobs = XMAP_OBS_POS(XCsc,rLng)		! Point-P at rP
			
C			if(NL.le.10.or.NL.ge.17244) print *, XCobs, Radius, rP, iVel, Sun__Spiral

			XCobs = XCobs+(Radius-rP)/iVEL*Sun__Spiral ! Point-P traced to Radius
			
		        call Julian8(10,iYr,Doy8,xTT,JEpoch)
			
						! All XE (East+West), only positive XE (East) or only negative XE (West)
		        Rdum = (1-abs(iEorW))*abs(rEloP)+iEorW*rEloP
		
			NLLL = NLLL + 1

			if(xTT .gt. MJDref) then
			  if(NLLLo.le.10.and.NLLL.lt.10) then
			    write( *,'(A,/,2F11.2,2I7)') 'xTT limit exceeded at xTT, MJDref, NLLL, NL =', xTT,MJDref,NLLL,NL
			    if(NLLLo.eq.10) write(*,'(A)')'If no files are read, most probably you have specified the wrong directory,'
			    if(NLLLo.eq.10) write(*,'(A,I6,A,F7.2,A)') 'or there are no files for iYr = ', iYrrr,' Doy = ', 
     &				Doy8, ' in the directory specified.'
			    if(NLLLo.eq.10) write(*,'(A,A)') 'This is especially true if the value of the year or Doy above looks strange.'
			    if(NLLLo.eq.10) write(*,'(A,A)') 'The directory specified is ',cWildCard
			  end if
			  NLLLo = NLLLo + 1
			end if
			if(NLLLo.ge.1000) return				! Limit on the number of bad reads possible

     		        if  ((XElow .le. Rdum .and. Rdum .le. XEhigh) .and.	! Select location according to criterion
     &			  ( abs(dRA) .le. XRlim ) .and.(ZSMEI(L) .gt. -999.) .and.
     &			  (MJDref .eq. 0 .or. xTT .le. MJDref)
     &		          ) then

     			  NL = NL + 1

			  if(NL.le.NLmax) then
			    cgbfiles(NL) = cgbfilesL			! File name line of sight came from
			    SRCB(NL) = SRCBS				! LOS name (Camera, RA & DEC and Orbit value)
			    IYRS(NL) = iYr				! Year of observation
			    DOYS8(NL)= Doy8				! DOY of observation
			    DIST(NL) = sngl(dDisSun)			! Sun - SMEI distance
			    XLS (NL) = sngl(dLngSun)			! Geocentric ecliptic longitude of Sun (deg)
			    XPE (NL) = rLngP				! Geocentric ecliptic longitude of point P (deg) (relative to Sun)		  	
			    XLE (NL) = rLatP				! Geocentric ecliptic latitude  of point P (deg)
			    XCE (NL) = XCsc				! Carrington variable sub-Earth point 			
			    XE  (NL) = rEloP				! Elongation of LOS
			    XC  (NL) = XCobs				! Point-P traced to Radius
			    XL  (NL) = xLat				! Heliographic lat(P)
C			    ZZ  (NL) = ZSMEI(L)/0.55			! Brightness in S10 (1 S10 = 0.55 ADU)			    
C			    ZZ  (NL) = ZSMEI(L)/0.40			! Brightness in S10 (1 S10 = 0.40 ADU)
C			    ZZ  (NL) = ZSMEI(L)				! Brightness in S10 (1 S10 = 1.0 ADU)
C			    ZZ  (NL) = ZSMEI(L)/0.70			! Brightness in S10 (1 S10 = 0.55 ADU + ~20% more e-)

C			    ZZ  (NL) = ZSMEI(L)/0.48			! Brightness in S10 (1 S10 = 0.48ADU + no extra e- for mass) best guess June, 2007

C			    ZZ  (NL) = ZSMEI(L)/0.58			! Brightness in S10 (1 S10 = 0.48 ADU + ~20% more e-) best guess Jul, 2007
C			    ZZ  (NL) = ZSMEI(L)/0.60			! Brightness in S10 (1 S10 = 0.50 ADU + ~20% more e-) best guess May 19, 2007
C			    ZZ  (NL) = ZSMEI(L)/0.65			! Brightness in S10 (1 S10 = 0.52 ADU + ~20% more e-)
C			    ZZ  (NL) = ZSMEI(L)/0.50			! Brightness in S10 (1 S10 = 0.50 ADU --- real brightness!)

C			    ZZ  (NL) = ZSMEI(L)/0.46			! Brightness in S10 (1 S10 = 0.46 ADU) SPIE Aug, 2007
			    ZZ  (NL) = ZSMEI(L)/0.552			! Brightness in S10 (1 S10 = 0.46 ADU + ~20% more e-) SPIE Aug, 2007

C***************************************************************************
C	  The below was put in for the 20 January 2005 shock reconstruction
C
C	  print *, doy8, rElop, ZZ(NL)
C	  if(rEloP.gt.85.0.and.rEloP.lt.95.0) then	! east of the Sun
C	    if(ZZ(NL).gt.-100.) then
C	      rEElop  = rEElop  + rEloP			! count all eastern elongations
C	      arEEloP = arEElop + 1.0
C	      ZZZZE   = ZZZZE   + ZZ(NL)
C	      AZZZZE  = AZZZZE  + 1.0
C	      Ddoy8E  = Ddoy8E  + doy8
C	      Adoy8E  = Adoy8E  + 1.0
C	    end if
C	  end if
C	  if(rEloP.lt.-85.0.and.rEloP.gt.-95.0) then	! west of the Sun
C	    if(ZZ(NL).gt.-100.) then
C	      rWElop  = rWElop  + rEloP			! count all eastern elongations
C	      arWEloP = arWElop + 1.0
C	      ZZZZW   = ZZZZW   + ZZ(NL)
C	      AZZZZW  = AZZZZW  + 1.0
C	      Ddoy8W  = Ddoy8W  + doy8
C	      Adoy8W  = Adoy8W  + 1.0
C	    end if
C	  end if
C***************************************************************************

			    if(NL.le.10) then
C			      print *, 'Into if ', SRCB(NL),DOYS8(NL),DIST(NL),XLS(NL),XPE(NL),XLE(NL),XCE(NL),
C     &						 XE(NL),XC(NL),XL(NL),ZZ(NL)
     			    end if
			  end if
			end if
			
		      end if
		    end do

C*******************
C	The below was put in for the 20 January 2005 shock reconstruction
C
C	rEEloP  = rEEloP/arEElop
C	ZZZZE   = ZZZZE/AZZZZE
C	rWEloP  = rWEloP/arWElop
C	ZZZZW   = ZZZZW/AZZZZW
C	Ddoy8E  = Ddoy8E/Adoy8E
C	Ddoy8W  = Ddoy8W/Adoy8W
C	print *, M, Ddoy8E, rEEloP, ZZZZE, AZZZZE, Ddoy8W, rWEloP, ZZZZW, AZZZZW
C
C*******************
	            if(NL.ge.NLMAX) bMore = .FALSE.
	    end if
	end do
	if(.not.bMore) then
	    write (*,'(1X,A,I7/)') 
     &	    'The number of SMEI brightness values maxed out at', NL
            NL = NLMAX
	end if
	write (*,'(1X,A,I8,A,I8,A/)') 
     &	    'There were',NL,' SMEI brightness values accepted, of a total',NLLL,' read'

	return
	end
