C+
C NAME:
C	READ_HOS
C PURPOSE:
C	Read Helios data for synoptic map
C CALLING SEQUENCE:
	subroutine READ_HOS(cWildCard,nCar,JDCar,XCbeg,ICnr,
     &		Radius,Speed,Eref,PowV, XElow,XEhigh,iEorW,XRlim,
     &		iSc,iHOS,cHOS,iSectBeg,iSectEnd,CC,iDay,nDev,dSig,
     &		NL,XC,YL,ZZ,
     &		cVel,LngV,LatV,VSol,VSol2,NSIDE)
C INPUTS:
C	cWildCard	character*(*)	wildcard specification for Helios files
C					'%%%%%_%%%_####.HOS', where #### is the
C					position of the rotation number
C	nCar		integer		# rotation start times
C	JDCar(nCar)	real*8		rotation start times (Julian days)
C	XCbeg		real		Carrington variable
C	ICnr		integer
C	Radius		real		reference distance (AU)
C	Speed		real		traceback speed (km/s)
C					Speed .eq. 0: IPS velocities are used for traceback
C					Speed .ne. 0: traceback at constant speed
C	Eref		real		reference elongation used for rescaling intensities
C					Eref .eq. 0: no rescaling
C	PowV		real		(currently effective only if Speed .eq. 0)
C					output ZZ is multiplied by V^PowV, where
C					V is the UCSD IPS speed (i.e. allows
C					calculation of mass flux, etc.)
C	XElow		real		minimum elongation limit
C	XEhigh		real		maximum elongation limit
C	iEorW		integer		-1,0,1 for west only/east+west/east only
C	XRlim		real		limit on ecliptic longitude rel. to Sun
C	iSc		integer		1 or 2 for Helios A or B
C	iHOS		integer		max. permitted # Helios files to be read
C	iSectBeg	integer		modified start sector number
C	iSectEnd	integer		modified end   sector number
C
C	The following three parameters are used by href=T_FILTER= to subtract a
C	running mean from the time series.
C
C	iDay		integer		time window
C	nDev		integer		glitch level (# standard deviations)
C	dSig		real		extra fraction of standard deviation added to
C					running mean
C
C	cVel		character*(*)	wildcard '$IPS:V1234_8.005' for
C					locating the IPS V-files used for
C					traceback if Speed=0
C OUTPUTS:
C	iHOS		integer		actual # Helios files read
C	cHOS(iHOS)	character*(*)	names of Helios files
C	CC		integer		1,2,3 for U,B,V
C
C	The following three may have been updated by interactive prompts
C	if the input value of iDay was less than zero (iDay < 0):
C
C	iDay		integer
C	nDev		integer
C	dSig		real
C
C	NL		integer		# data points
C	XC		real		modified Carrington variable Point P after traceback
C	YL		real		heliographic latitude
C	ZZ		real		intensity (S10) units (scaled to Eref 
C					if Eref .ne. 0; otherwise unscaled)
C					multiplied by V^PowV, where V is the
C					UCSD IPS speed.
C	cVel		character*(*)	set to IPS velocity file
C					'$IPS:V1234_8.005'
C CALLS:
C	BadR4, iGetLogical, Say, ECLIPTIC_HELIOGRAPHIC, GAL_CNTR, Julian, HOSRead
C	iSearch, ArrR4Copy
C	SC_ECLIP, SC_ECLIP90, ThomsonUBVConst, ThomsonLOSS10Far, T_FILTER, XMAP_SC_POS
C	XMAP_OBS_POS, HOSOrbID, iGridScan, iFltArr, ArrR4TimesConstant, MKVTRACE
C INCLUDE:
	include		'hosucsd.h'		! Array size parameters
	include		'sun.h'
	include		'dirspec.h'
	include		'hos_e9.h'
C EXTERNAL:
	external	HELIOS_1
	external	HELIOS_2
	external	fncWZ			! not used
C PROCEDURE:
C MODIFICATION HISTORY:
C	FEB-1996, Paul Hick (UCSD); extracted from SANHEL.FOR
C-
	    character	cWildCard*(*)
	    integer	nCar
	    real*8	JDCar(nCar)
	    real	XCbeg
	    integer	ICnr
	    real	Radius
	    real	Speed
	    real	Eref
	    real	PowV
	    real	XElow
	    real	XEhigh
	    integer	iEorW
	    real	XRlim
	    integer	iSc
	    integer	iHOS
	    character	cHOS(iHOS)*(*)
	    integer	iSectBeg
	    integer	iSectEnd
	    integer	CC
	    integer	iDay
	    integer	nDev
	    real	dSig
	    integer	NL

	    !-------
	    ! NMAX is defined in hosucsd.h

	    real	XC(NMAX)		! Modified Carrington variable Point P
	    real	YL(NMAX)		! Heliographic latitude Point P (deg)
	    real	ZZ(NMAX)		! Intensity (S10)

c	    real	XLS(NMAX)
c	    real	XCH(NMAX)
c	    real	XOP(NMAX)
c	    real	XAP(NMAX)
c	    real	XRP(NMAX)
c	    real	XEP(NMAX)
c	    real	DOY(NMAX)

c	    integer	NYR(NMAX)

	    character	cVel*(*)
	    integer	LngV
	    integer	LatV
	    real	VSol (LngV,LatV,2)
	    real	VSol2(LngV,LatV,2)
	    integer*1	NSIDE(LngV,LatV,-1:1,-1:1)
						! Values not yet subjected to traceback:
						! (saved between calls; refilled if necessary)
	real		XRN   (NMAX)		! Heliocentric distance point P
	real		XCN   (NMAX)		! Carrington variable point P
	real		XEN   (NMAX)		! Topocentric elongation point P
	real		ZZN   (NMAX)		! Intensity
	real		VTrace(NMAX)		! IPS traceback speeds

	integer		NTrace(NMAX)		! Scratch array for MKVTRACE
	logical		bTrace			! Controls recalculation VTrace

	!-------
	! MXC, NSECT, NBUF are defined in hosucsd.h
	! *HOS arrays are use to store data read from the Helios data files by HOSRead

	integer		PHOS(MXC)		! Photometer ID
	integer		NHOS(MXC)		! Record numbers		NOT USED
	integer		KHOS(MXC)		! Colors			NOT USED
	integer		FHOS(MXC)		! Filters			NOT USED
	real		THOS(MXC)		! Times (doy of year)
	real		RHOS(MXC)		! Distance s/c-Sun (AU)
	real		LHOS(MXC)		! Ecliptic longitude Sun (deg)
	real		ZHOS(NSECT*MXC)		! Intensities (S10)

	real		SX (NSECT*NBUF)		! Scratch arrays for T_FILTER
	real		SNX(NSECT*NBUF)
	real		SXX(NSECT*NBUF)

	!-------
	! Filled by SC_ECLIP: relative heliocentric (H* arrays) and topocentric
	! (T* arrays) ecliptic coordinates of point P for each line of sight

	real		HLNG(NBSECT:NESECT,2)	! Helioc. ecl. lng(P) (deg)
	real		HLAT(NBSECT:NESECT,2)	! Helioc. ecl. lat(P) (deg)
	real		HR  (NBSECT:NESECT,2)	! Distance P-Sun (unit Sun-S/C dist)
	real		HELO(NBSECT:NESECT,2)	! Angle P-Sun-S/C (deg)

	real		TLNG(NBSECT:NESECT,2)	! Longitude (deg)
	real		TLAT(NBSECT:NESECT,2)	! Latitude (deg)
	real		TR  (NBSECT:NESECT,2)	! Distance P-S/C (unit Sun-S/C dist)
	real		TELO(NBSECT:NESECT,2)	! Angle P-S/C-Sun (elongation; deg)

	real*8		JEpoch
	character	cFile*80
	character	cSc(2)		/'A','B'/
	character	cSay*7		/'READ_HOS'/

	integer		IPSD		/5/
	integer		nPad(4)		/4*0./

	!integer	kPP(3)		/HOS__P_1,HOS__P_2,HOS__P_3/
	integer		kCC(3)		/HOS__C_1,HOS__C_2,HOS__C_3/
	!integer	kFF(5)		/HOS__F_1,HOS__F_2,HOS__F_3,HOS__F_4,HOS__F_5/

	logical		bDataAvailable	/.FALSE./
	logical		bMore
	logical		bRead
	logical		bOneCC
	logical		bOneFF

	save		XRN, XCN, XEN, ZZN, NTrace, VTrace, bTrace

	save		RadiusOld, SpeedOld, ErefOld, XElowOld	 , XEhighOld,
     &			iEorWOld , XRlimOld, iScOld , iSectBegOld, iSectEndOld,
     &			iDayOld  , nDevOld , dSigOld

	bOneCC(kC) = kC .eq. HOS__C_1 .or. kC .eq. HOS__C_2 .or. kC .eq. HOS__C_3
	bOneFF(kF) = kF .eq. HOS__F_1 .or. kF .eq. HOS__F_2 .or. kF .eq. HOS__F_3 .or. kF .eq. HOS__F_4 .or. kF .eq. HOS__F_5

	Bad = BadR4()

	!-------
	! New data have to be read if there are no data available from a previous call
	! or if any of the parameters involved in the data selection changed

	bRead = .not. bDataAvailable		.or.
     &		XCbegOld	.ne. XCbeg	.or.
     &		iCnrOld		.ne. iCnr	.or.
     &		XElowOld	.ne. XElow	.or.
     &		XEhighOld	.ne. XEhigh	.or.
     &		iEorWOld	.ne. iEorW	.or.
     &		XRlimOld	.ne. XRlim	.or.
     &		iScOld		.ne. iSc	.or.
     &		iSectBegOld	.ne. iSectBeg	.or.
     &		iSectEndOld	.ne. iSectEnd	.or.
     &		iDayOld		.ne. iDay	.or.
     &		nDevOld		.ne. nDev	.or.
     &		dSigOld		.ne. dSig

	!-------
	! If Speed .ne. 0 (traceback at constant speed) don't multiply intensity
	! with V^P (would be multiplication by a constant).

	if (Speed .ne. 0) PowV = 0.

	!-------
	! Changes in Radius, Speed, Eref and PowV do not require reading of data files
	! and are treated separately.

	if (.not. bRead 			.and.
     &		RadiusOld	.eq. Radius	.and.
     &		SpeedOld	.eq. Speed	.and.
     &		ErefOld		.eq. Eref	.and.
     &		PowVOld		.eq. PowV) return

	iBeg90 = 1
	iEnd90 = 1					! Read only 90 deg intensities
	i90 = iEnd90-iBeg90+1

	if (bRead) then					! If new data are read from file

	    !-------
	    ! Get relative heliocentric and topocentric ecliptic coordinates of point P
	    ! for all lines of sight (photometer-sector combinations)

	    call SC_ECLIP  ( iSc,NBSECT,NESECT,HLNG,HLAT,HR,HELO)	! Heliocentric
	    call SC_ECLIP  (-iSc,NBSECT,NESECT,TLNG,TLAT,TR,TELO)	! Topocentric
	    call SC_ECLIP90( iSc,HLNG90,HLAT90,HR90,HELO90)		! Heliocentric
	    call SC_ECLIP90(-iSc,TLNG90,TLAT90,TR90,TELO90)		! Topocentric

	    iHOSMax = iHOS
	    iHOS   = 0
	    NL	   = 0					! Count data points

	    IC	   = XCbeg
	    XClast = XCbeg+ICnr-1

	    bDataAvailable = .TRUE.
	    bTrace 	   = .TRUE.			! VTrace array has to rebuild

	end if

	bMore = bRead					! Read new data (skipped if bMore = .FALSE.)
	do while (IC .le. XClast .and. bMore)

	    I = iGetLogical(cEnvi//'HOS',cFile)
	    cFile(I+1:)    = cWildCard
	    cFile(I+1:I+1) = cSc(iSc)
	    write (cFile(I+11:I+14),'(I4.4)') IC

	    if (iSearch(1,cFile,cFile) .eq. 1) then
		call Julian(1,IAYr,RADoy,JDCar(IC  ),JEpoch)
	        call Julian(1,IBYr,RBDoy,JDCar(IC+1),JEpoch)

		iHOS = iHOS+1
		if (iHOS .le. iHOSMax) cHOS(iHOS) = cFile

		write (*,'(3X,2A,I4.4,A,2(I3.3,A))') cFile(:itrim(cFile)),'  (',IAYr,'  ',int(RADoy),'-',int(RBDoy),')'

		kPCF = HOS__P_12+HOS__C_ALL+HOS__F_ALL+HOS__SWAP_SECT	! Phot 1&2, swap sectors
		call HOSRead(kPCF,cFile,MXC,iCnt,NHOS,THOS,PHOS,KHOS,FHOS,RHOS,LHOS,NSECT,ZHOS)

		if (.not. bOneCC(kPCF) .or. .not. bOneFF(kPCF)) call Say(cSay,'E','OneCF',
     &			'only one color and filter allowed for 16/31 deg photometers')

		!--------
		! It's possible to switch from one color or filter to another when reading a new file.
		! Is this supposed to be this way???

		do I=1,3				! Find color
		    if (iand(kPCF,HOS__C_ALL) .eq. kCC(I)) CC = I
		end do

		!do I=1,5				! Find filter
		!    if (iand(kPCF,HOS__F_ALL) .eq. kFF(I)) iFilt = I
		!end do

		call HOSOrbID(THOS(1),RHOS(1),LHOS(1)+180.,iSc,iYr,dRMin)
		call GAL_CNTR(iSc,iCnt,NBSECT,NESECT,LHOS,PHOS,ZHOS)
		call T_FILTER(iCnt,NSECT,THOS,ZHOS,PHOS,iDay,nDev,dSig,NBUF,SNX,SX,SXX)

		ICoff = IC-XCbeg			! Map all points to rotation XCbeg

		I = 1
		do while (I .le. iCnt .and. bMore)		! Loop over all MC's
								! Carrington variable for S/C
		    if (iSc .eq. 1) XCHOS = XMAP_SC_POS(HELIOS_1,iYr,THOS(I),nCar,JDCar)
		    if (iSc .eq. 2) XCHOS = XMAP_SC_POS(HELIOS_2,iYr,THOS(I),nCar,JDCar)
		    XCHOS = XCHOS-ICoff

		    J = iSectBeg
		    do while (J .le. iSectEnd .and. bMore)	! Loop over sectors

			if (iEorW .eq. 0) then			! All XE (East+West)
			    rElo = abs(TELO(J,PHOS(I)))
			else
			    rElo = iEorW*TELO(J,PHOS(I))	! Only positive XE (East) or
			end if					! only negative XE (West)

			dRA = TLNG(J,PHOS(I))			! Topocentric lng(P)-lng(Sun)
			if (abs(dRA) .gt. 180.) dRA = dRA-sign(1.,dRA)*360.

			JHOS = (I-1)*NSECT+J-NBSECT+1

			if (ZHOS(JHOS) .ne. Bad .and.			! Skip flagged data
     &				XElow .le. rElo .and. rElo .le. XEhigh .and.
     &				abs(dRA) .le. XRlim) then		! Check elongation and RA range

			    bMore = NL .lt. NMAX
			    if (bMore) then
				NL = NL+1
									! Don't do traceback yet.
				XRN(NL) = HR(J,PHOS(I))*RHOS(I)		! Heliocentric dist(P) (AU)
				XEN(NL) = abs(TELO(J,PHOS(I)))		! Topocentric elongation P (deg)
				ZZN(NL) = ZHOS(JHOS)			! Intensity (S10)

				rLng = LHOS(I)+180.+HLNG(J,PHOS(I))	! Heliocentric ecliptic lng(P)
				rLat = HLAT(J,PHOS(I))			! Heliocentric ecliptic lat(P)
				call ECLIPTIC_HELIOGRAPHIC(0,iYr,THOS(I),rLng,rLat)
									! rLng: Heliographic lng(P)
									! rLat: Heliographic lat(P)
				XCN(NL) = XMAP_OBS_POS(XCHOS,rLng)	! Carrington variable point P
				YL (NL) = rLat				! Heliographic lat(P) (deg)

				!-------
				! May be needed for the deconvolution program (arrays haven't been declared yet)
				!
				! NYR(NL) = iYr				! Year of observation
				! DOY(NL) = THOS(I)			! Day of year
				! XLS(NL) = LHOS(I)			! Topocentric ecliptic lng(Sun)
				! XCH(NL) = XCHOS			! Carrington variable s/c
				! XOP(NL) = TLNG(J,PHOS(I))		! Topocentric ecliptic lng(P)
				! XAP(NL) = TLAT(J,PHOS(I))		! Topocentric ecliptic lat(P)
				! XRP(NL) = HR(J,PHOS(I))*RHOS(I)	! Heliocentric dist(P) (AU)
				! XEP(NL) = abs(TELO(J,PHOS(I)))	! Topocentric elongation P (deg)

			    else
				call Say(cSay,'W','RDALL','Not all data read; NMAX too small')

			    end if

			    ! print *,nl, xc(nl),yl(nl),xen(nl)

			end if

			J = J+1		! Next sector
		    end do

		    I = I+1		! Next MC
		end do

		!-------
		! Read the 90 degree photometer data.

		kPCF = HOS__P_3+kCC(CC)+HOS__F_ALL		! Phot 3, color CC
		call HOSRead(kPCF,cFile,MXC,iCnt,NHOS,THOS,PHOS,KHOS,FHOS,RHOS,LHOS,i90,ZHOS)
		call T_FILTER(iCnt,i90,THOS,ZHOS,PHOS,iDay,nDev,dSig,NBUF,SNX,SX,SXX)

		I = 1
		do while (I .le. iCnt .and. bMore)		! Loop over all MC's
								! Carrington variable for S/C
		    if (iSc .eq. 1) XCHOS = XMAP_SC_POS(HELIOS_1,iYr,THOS(I),nCar,JDCar)
		    if (iSc .eq. 2) XCHOS = XMAP_SC_POS(HELIOS_2,iYr,THOS(I),nCar,JDCar)
		    XCHOS = XCHOS-ICoff

		    J = iBeg90
		    do while (J .le. iEnd90 .and. bMore)	! Loop over int,pB

			rElo = TELO90				! Topocentric elongation (=90 deg)

			dRA = TLNG90				! Topocentric lng(P)-lng(Sun)
			if (abs(dRA) .gt. 180.) dRA = dRA-sign(1.,dRA)*360.

			JHOS = (I-1)*i90+J-iBeg90+1

			if (ZHOS(JHOS) .ne. Bad .and.			! Skip flagged data
     &				XElow .le. rElo .and. rElo .le. XEhigh .and.
     &				abs(dRA) .le. XRlim) then		! Check elongation and RA range

			    bMore = NL .lt. NMAX
			    if (bMore) then
				NL = NL+1
									! Don't do traceback yet.
				XRN(NL) = HR90*RHOS(I)			! Heliocentric dist(P) (AU)
				XEN(NL) = TELO90			! Topocentric elongation P (deg)
				ZZN(NL) = ZHOS(JHOS)			! Intensity (S10)

				rLng = LHOS(I)+180.+HLNG90		! Heliocentric ecliptic lng(P)
				rLat = HLAT90				! Heliocentric ecliptic lat(P)
				call ECLIPTIC_HELIOGRAPHIC(0,iYr,THOS(I),rLng,rLat)
									! rLng: Heliographic lng(P)
									! rLat: Heliographic lat(P)
				XCN(NL) = XMAP_OBS_POS(XCHOS,rLng)	! Carrington variable point P
				YL (NL) = rLat				! Heliographic lat(P) (deg)

				!-------
				! May be needed for the deconvolution program (arrays haven't been declared yet)
				!
				! NYR(NL) = iYr				! Year of observation
				! DOY(NL) = THOS(I)			! Day of year
				! XLS(NL) = LHOS(I)			! Topocentric ecliptic lng(Sun)
				! XCH(NL) = XCHOS			! Carrington variable s/c
				! XOP(NL) = TLNG90			! Topocentric ecliptic lng(P)
				! XAP(NL) = TLAT90			! Topocentric ecliptic lat(P)
				! XRP(NL) = HR90*RHOS(I)		! Heliocentric dist(P) (AU)
				! XEP(NL) = TELO90			! Topocentric elongation P (deg)

			    else
				call Say(cSay,'W','RDALL','Not all data read; NMAX too small')

			    end if

			    ! print *,nl, xc(nl),yl(nl),xen(nl)

			end if

			J = J+1		! Next sector
		    end do

		    I = I+1		! Next MC
		end do

	    end if

	    IC = IC+1			! Next file
	end do

	if (NL .eq. 0) return		! No data read

	call ArrR4Copy(NL,ZZN,ZZ)	! Copy unscaled intenisities to ZZ

	if (Eref .ne. 0) then		! Scale intensities to elongation Eref
	    call ThomsonUBVConst(CC,UU,APM)
	    Tmp = ThomsonLOSS10Far(0.,-1.,sngl(Radius/sind(Eref)/SUN__RAU),Eref,APM,PP)
	    do I=1,NL
		ZZ(I) = ZZ(I)*Tmp/ThomsonLOSS10Far(0.,-1.,sngl(XRN(I)/sind(XEN(I))/SUN__RAU),XEN(I),APM,PP)
	    end do
  	end if

	if (Speed .eq. 0) then		! IPS velocity traceback
	    if (bTrace) then		! Refresh VTrace array
		I = XCbeg
		call MKVTRACE(0.5d0*(JDCar(I+ICnr)+JDCar(I)),cVel,IPSD,
     &			NL,XRN,YL,XCN,NTrace,VTrace, LngV,LatV,VSol,VSol2,NSIDE)
		bTrace = .FALSE.	! VTrace valid until new data are read
	    end if

	    do I=1,NL			! Perform traceback to Radius
		XC(I) = XCN(I)+SUN__SPIRAL*(Radius-XRN(I))/VTrace(I)
	    end do

	    if (PowV .ne. 0.) then	! Flux calculation: I*VTrace^PowV
		do I=1,NL
		    ZZ(I) = ZZ(I)*VTrace(I)**PowV
		end do
	    end if

	else				! Traceback at constant speed
	    do I=1,NL
		XC(I) = XCN(I)+SUN__SPIRAL*(Radius-XRN(I))/Speed
	    end do
					! Flux calculation: I*Speed^PowV
	    if (PowV .ne. 0.) call ArrR4TimesConstant(NL,ZZ,Speed**PowV,ZZ)

	end if

	XCbegOld	= XCbeg		! Save all parameters involved
	iCnrOld		= iCnr		! .. in the data selection.
	XElowOld	= XElow
	XEhighOld	= XEhigh
	iEorWOld	= iEorW
	XRlimOld	= XRlim
	iScOld		= iSc
	iSectBegOld	= iSectBeg
	iSectEndOld	= iSectEnd
	iDayOld		= iDay
	nDevOld		= nDev
	dSigOld		= dSig

	RadiusOld	= Radius
	SpeedOld	= Speed
	ErefOld		= Eref
	PowVOld		= PowV

	if (Speed .eq. 0) then		! In case of IPS Traceback

	    !-------
	    ! Read velocity map near Radius. Needed for blanking out areas in the map
	    ! where no velocity data are available (VSol and NSIDE are used in calling 
	    ! program only if Speed .eq. 0)

	    nIPS = int(Radius*100/IPSD)*IPSD
	    nIPS = max(1,nIPS)

	    write (cVel(15:15+2),'(I3.3)') nIPS
	    cFile = cVel
	    I = iFltArr(cFile,40,0.,1.,fncWZ,BadR4(),LngV*LatV,LngV,LatV,nPad,VSol)
	    I = iGridScan(LngV,LatV,VSol,BadR4(),NSIDE)

	end if

	return
	end