C+
C NAME:
C	ReadG
C PURPOSE:
C	Read data from unedited G*. file or edited E*.DAT file,
C	containing Cambridge IPS data
C CATEGORY:
C	Data processing
C CALLING SEQUENCE:
	subroutine ReadG(cDat,nCar,JDCar,MJD,iEdt,Radius,Speed,GPow, iOK,XCmin,XCmax,
     &			 IYRS,DOYS, XDS,XLS,XLL,XDL,XCE,XC,YL,XE,XR,XV,XG,cStr,bMessage)
C INPUTS:
C	cDat		character*(*)	name of logical for dir containing data 
C					files
C	nCar		integer		(read-only
C					# start times for Carrington rotations
C	JDCar(nCar)	real*8		(read-only)
C					Carrington start times (Julian days)
C	MJD		integer		(read-only)
C					modified Julian day (JD-2400000.5) for
C					the requested file
C	iEdt		integer		(read-only)
C					0 = Read unedited data file
C					1 = Read edited data file
C	Radius		real		Reference distance in AU
C					(not used if Speed=0)
C	Speed		real		trace back velocity (km/s)
C					(if Speed=0 then no traceback is done)
C	GPow		real		if non-zero used to scale the measured
C					G-value: G = G(measured)/RP^GPow
C					where RP is the heliocentric distance 
C					of the point P.
C					(i.e. if GPow=0 the unscaled value is 
C					returned)
C	bMessage	logical		if .FALSE. suppresses Say call before return
C OUTPUTS:
C	(NTIME=72,NDEC=14)
C
C	iOK		integer		= 0 if file was not opened succesfully
C					= 1 if file was opened
C					= 3 if file contained no valid points
C	XCmin		real		Minimum XC value (BadR4() if file empty)
C	XCmax		real		Maximum XC value (BadR4() if file empty)
C
C	IYRS(NTIME,NDEC)integer		time: year
C	DOYS(NTIME,NDEC)real		time: doy of year
C	XDS(NTIME,NDEC)	real		Earth-Sun distance (AU)
C	XLS(NTIME,NDEC) real		Geocentric ecliptic longitude Sun (degrees)
C	XLL(NTIME,NDEC) real		Difference between geocentric ecliptic longitudes of
C					line of sight and Sun
C	XDL(NTIME,NDEC) real		Geocentric ecliptic latitude of line of sight
C	XCE(NTIME,NDEC) real		Modified Carrington variable of Earth
C
C	XC(NTIME,NDEC)	real		Modified Carrington variable of 'Point P', except for
C					an offset of NCoff, i.e. the conventional Carrington
C					numbers are NCoff+XC(I)
C	YL(NTIME,NDEC)	real		Heliographic latitudes (degrees) of 'Point P'
C	XE(NTIME,NDEC)	real		Elongation of the line of sight used
C					If XE<0	the l.o.s. points west of the Sun
C					If XE>0 the l.o.s. points east of the Sun
C	XR(NTIME,NDEC)	real		Right Ascension relative to Sun
C	XV(NTIME,NDEC)	real		velocities (currently set to zero)
C	XG(NTIME,NDEC)	real		G-factors
C		Each observation XG refers to a line of sight extending from Earth.
C		The coordinates XC,YL define the heliocentric position of the point
C		P along the line of sight closest to the Sun. abs(XE) is the angle
C		between line of sight and the direction to the center of the solar disk.
C	cStr		character*(*)
C CALLS:
C	ECLIPTIC_EQUATOR, ECLIPTIC_HELIOGRAPHIC, SunNewcomb, Julian, DATE_DOY
C	iSearch, iFilePath, bOpenFile, iFreeLun, ReadGHD, Say, POINT_ON_LOS, T_GST,
C	ArrR4Bad, ArrI4Bad, ArrR4GetMinMax
C EXTERNAL:
	external	EARTH
C INCLUDE:
	include		'openfile.h'
	include		'sun.h'
C RESTRICTIONS:
C >	Input files are searched for in the directory assigned to logical $DATA
C >	Make sure to test iOK before using the output arrays (the arrays are
C	not always cleared before returning; you may end up using the values
C	from a previous call to ReadG)
C PROCEDURE:
C	Unedited data files: G12345.   where 12345 is the modified Julian day
C
C    >	The data files contain for each day 72 records (i.e. 72 observation
C	times, evenly spaced over the whole day) of 14 declinations.
C    >	The original G-factors are coded as hexidecimal numbers. The 
C	translation is done using the Z format specifier
C    >  The final G-factors are calculated from
C		G = GMAX**(IG/G0-1.)
C	where GMAX and G0 are constants given in the data file and IG is the
C	translated hexidecimal number
C    >  Given the UT and the equatorial coordinates (HA=0, declination 
C	specified in data file), the heliocentric coordinates of the point P
C	are calculated.
C    >  If no G-factor is available for a line of sight the G-factor and
C	heliocentric coordinates are set to BadR4() (real*4) or BadI4 (int*4)
C
C	Edited data files: E12345.DAT  where 12345 is the modified Julian day
C
C    >  Each ASCII file contains 73 records. The first repeats the modified
C	Julian day number. The other 72 records each contain 14 G-factors
C	in F8.3 format. The organization is the same as for the unedited files.
C MODIFICATION HISTORY:
C	MAR/APR-1992, Paul Hick (UCSD)
C	JUL-1993, Paul Hick (UCSD); fixed bug in equatorial-to-ecliptic
C		coordinate conversion
C	FEB-1998, Paul Hick (UCSD); added arguments IYRS,..,XCE for support of
C		IPS deconvolution program.
C-
	parameter	(NTIME	= 72)
	parameter	(NDEC	= 14)
	parameter	(NTD    = NTIME*NDEC)

	    character	cDat*(*)
	    integer	nCar
	    real*8	JDCar(nCar)
	    integer	MJD
	    integer	iEdt
	    real	Radius
	    real	Speed
	    real	GPow
	    integer	iOK
	    real	XCmin
	    real	XCmax

	    integer	IYRS(NTIME,NDEC)	! Output arrays
	    real	DOYS(NTIME,NDEC)
	    real	XDS (NTIME,NDEC)	! Earth-Sun distance
	    real	XLS (NTIME,NDEC)	! Ecliptic longitude Sun
	    real	XLL (NTIME,NDEC)	! Ecliptic longitude difference los and Sun-Earth direction
	    real	XDL (NTIME,NDEC)	! Ecliptic latitude los
	    real	XCE (NTIME,NDEC)	! Carrington variable of Earth

	    real	XC  (NTIME,NDEC)	! Modified Carrington variable of Point P
	    real	YL  (NTIME,NDEC)	! Heliographic latitude of Point P
	    real	XE  (NTIME,NDEC)	! Elongation
	    real	XR  (NTIME,NDEC)	! Right ascension relative to Sun
	    real	XV  (NTIME,NDEC)	! Velocities
	    real	XG  (NTIME,NDEC)	! G-factors

	    character	cStr*(*)
	    logical	bMessage

	character	cFile*80		! Input file name
	character	cMon*3
	character	cZEdt(0:1)*8	/'(14Z2.0)','(14F8.3)'/
	character	cEdt(0:1)	/'g','e'/
	character	cExt(0:1)*4	/'.   ','.dat'/
	character	cDum*5

	integer		IG(NDEC)
	real		GG(NDEC)
	equivalence	(IG,GG)
	
	real		DecRef(14) /-7.8, 1.9, 9.8,16.9,23.4,29.5,35.4,
     &				     41.0,46.6,52.1,57.7,63.2,68.9,74.8/
	real		CAMBRIDGE	/-.095/	! Geographic long Cambridge (degrees)

	real*8		JEpoch
	real*8		dLngSun
	real*8		dLatSun
	real*8		dDisSun			 ! Args to SunNewcomb		
	logical		bEOH
	logical		bOpenFile

	if (iEdt .ne. 0 .and. iEdt .ne. 1) call Say('ReadG','E','INVTYP','invalid type: must be 0/1')

	write (cDum,'(I5.5)') MJD
	I = iFilePath(cDat,0,' ',cEdt(iEdt)//cDum//cExt(iEdt),cFile)

	iOK = 0
	if (.not. bOpenFile(OPN__REOPEN+OPN__TEXT+OPN__READONLY,iU,cFile,iRecl)) return! Open data file

	iOK = 3				! File opened, but may be empty

	if (iEdt .eq. 1) then
	    read (iU,'(A)',iostat=I) cDum
	    bEOH = I .eq. 0
	else				! Read file header section
	    call ReadGHD(MJD,iU,cFile, bEOH,GMAX,G0)
	end if

	if (.not. bEOH) then		! Error reading header: probably empty file
	    iU = iFreeLun(iU)
	    call Say('ReadG','W',cFile,'is an empty file')
	    return			! Return with iOK = 3
	end if
	
	call Julian(1,iYr,Doy,SUN__MJDtoJD+MJD,JEpoch)
	iDoy = nint(Doy)
	call DATE_DOY(1,iYr,cMon,iMon,iDay,iDoy)

	call ArrI4Bad(NTD, IYRS)	! If no G-value available
	call ArrR4Bad(NTD, DOYS)
	call ArrR4Bad(NTD, XDS )
	call ArrR4Bad(NTD, XLS )
	call ArrR4Bad(NTD, XLL )
	call ArrR4Bad(NTD, XDL )
	call ArrR4Bad(NTD, XCE )
	call ArrR4Bad(NTD, XC  )
	call ArrR4Bad(NTD, YL  )
	call ArrR4Bad(NTD, XE  )
	call ArrR4Bad(NTD, XR  )
	call ArrR4Bad(NTD, XV  )
	call ArrR4Bad(NTD, XG  )

	!-------
	! Read NTIME records from data file. The UT times assigned to each record
	! are evenly spaced over the 24 hours of the day. Each record is an average
	! over a fraction 1/NTIME of the day. The center time of each interval is
	! used to calculate various coordinates.

	write (cZEdt(iEdt)(2:3),'(I2.2)') NDEC
	NP = 0
	do I=1,NTIME

	    dUT = (I-.5)/NTIME		! Center time (fraction of day)
	    Doy = iDoy+dUT		! UT time of day (fraction for time of day)
	    Sec = 86400.*dUT		! UT time of day (seconds)

	    call SunNewcomb(0,iYr,Doy,dLngSun,dLatSun,dDisSun)
	    rLngSun = dLngSun

	    XCEarth = XMAP_SC_POS(EARTH,iYr,Doy,nCar,JDCar) ! Carrington var Earth
	    RALOS = T_GST(iYr,iDoy,Sec)+CAMBRIDGE ! RA (deg) = LST (Hour Angle is zero)

	    RASun = rLngSun		! Get RA of Sun
	    rLat  = dLatSun
	    call ECLIPTIC_EQUATOR(0,iYr,Doy,RASun,rLat)
	    RASun = RALOS-RASun
	    if (abs(RASun) .gt. 180) RASun = RASun-sign(1.0,RASun)*360.

	    if (iEdt .eq. 0) read (iU,cZEdt(iEdt)) IG	!(IG(J),J=1,NDEC) ! Read rec of hex nrs from file
	    if (iEdt .eq. 1) read (iU,cZEdt(iEdt)) GG	!(GG(J),J=1,NDEC)

	    !-------
	    ! Proces the NDEC declinations in each record:
	    ! -	Calculate G-factor
	    ! -	Calculate heliographic coordinates of the point P (nearest to the Sun
	    !	along the line of sight

	    do J=1,NDEC

		if ( (iEdt .eq. 0 .and. IG(J) .gt. 0)   .or.
     &		     (iEdt .eq. 1 .and. GG(J) .gt. 0) ) then	! ????????

		    rLng = RALOS
		    rLat = DecRef(J)
		    call ECLIPTIC_EQUATOR(1,iYr,Doy,rLng,rLat)	! Ecliptic coord. los
		    rLng = rLng-rLngSun		! Longitude diff rel Earth-Sun

		    IYRS(I,J) = iYr
		    DOYS(I,J) = Doy
		    XDS	(I,J) = dDisSun
		    XLS (I,J) = rLngSun
		    XLL (I,J) = rLng
		    XDL (I,J) = rLat
		    XCE (I,J) = XCEarth
				 		! Dist. Earth-PointP (not closer than .25 AU)
		    RP	 = max(.25, cosd(rLat)*cosd(rLng) )
		    call POINT_ON_LOS(rLng,rLat,RP,ELO,ELOSUN,iEorW)
						! iEorW: -1=West; +1=East
		    rLng = rLng+rLngSun-180	! Ecliptic long point P
		    call ECLIPTIC_HELIOGRAPHIC(0,iYr,Doy,rLng,rLat)
						! Heliographic coord point P
		    XC(I,J) = XMAP_OBS_POS(XCEarth,rLng)
		    if (Speed .ne. 0) XC(I,J) = XC(I,J)+SUN__SPIRAL*(Radius-RP)/Speed

		    YL(I,J) = rLat		! Heliographic lat point P
		    XE(I,J) = iEorW*ELO		! Elongation (sign for E/W)
		    XR(I,J) = RASun		! RA relative to Sun

		    XV(I,J) = 0.
		    if (iEdt .eq. 0 ) XG(I,J) = GMAX**(IG(J)/G0-1.) ! G-factor
		    if (iEdt .eq. 1 ) XG(I,J) = GG(J)
		    if (GPow .ne. 0.) XG(I,J) = XG(I,J)/RP**GPow
		    NP = NP+1

		end if

	    end do

	end do

	iU = iFreeLun(iU)			! Close data file

	call ArrR4GetMinMax(-NTD,XC,XCmin,XCmax)! Get min,max (excl. bad values)
						! cStr is output argument
	write (cStr,'(A,I4,3A,I4,I6,A)') cFile(:itrim(cFile)),iDay,'-',cMon,'-',iYr,NP,' pnts'
	if (bMessage) call Say(' ',' ',' ',cStr)
	if (NP .ne. 0) iOK = 1			! File contains valid data

	return
	end