C+
C NAME:
C	ReadVIPSn8
C	ReadVIPSLOSCheckn
C PURPOSE:
C	Read IPS velocity data into arrays
C CATEGORY:
C	I/O
C CALLING SEQUENCE:
C	call ReadVIPSn8(iReadVIPSn,iProcessVIPSn,cWildVIPS,
C     &		nCar,JDCar,bAuto,XCtest1,XCtest2,NCoff,MJDref,
C     &		Radius,
C     &		XElow,XEhigh,iEorW,XRlim,
C     &		NLmax,NL,cbvfiles,SRCVG,SRC,SRCVGsav,SRCsav,
C     &		IYRF,IREC,IYRS,DOYS8,XDS,XLS,XLL,XDL,XCE,XE,XC,YL,VV,GG,
C     &		NSmax,IYRFsav,IRECsav,IYRSsav,DOYSsav8,XDSsav,XLSsav,XLLsav,
C     &		XDLsav,XCEsav,XEsav,XCsav,YLsav,VVsav,GGsav)
C
C	call ReadVIPSLOSCheck(XCbeg,XCend,NLOS,NINSIDE,DLOS)
C INPUTS:
C   ReadVIPS:
C	iReadVIPSn8	external integer function
C	iProcessVIPSn	external integer function
C					external functions control data file access
C	cWildVIPS	character*(*)	wildcard used to locate data files;
C					the wildcard must contain the substring 
C					'%%'
C	nCar		integer		dimension of JDCar
C	JDCar(nCar)	integer		start times (Julian days) of Carrington rotations
C	bAuto		logical
C	XCtest1		real		modified Carrington variable for start of search
C	XCtest2		real		modified Carrington variable for end of search
C	NCoff		integer		JDCar(I) is the start of rotation NCoff+I
C	MJDref		real*8		if not equal 0, only data prior to
C					modified Julian day MJDref are used
C	Radius		real		reference distance (AU)
C
C	XElow		real
C	XEhigh		real
C	iEorW		integer
C	XRlim		real
C
C	NLmax		integer		max. # points read into arrays 
C					XC,YL,VV,GG, XCsav,YLsav,VVsav,GGsav
C	NSmax		integer
C
C    ReadVIPSLOSCheck:
C	XCbeg		real		modified Carrington variable for start of map
C	XCend		real		modified Carrington variable for end of map
C	NLOS		integer
C	NINSIDE		integer
C	DLOS
C OUTPUTS:
C	NL		integer		# valid data points in arrays XC,YL,VV,GG
C	SRC(NLmax)	character*8	Character Source identifier
C	IYRF(NLmax)	integer		year of file from which observation was extracted (should be same as IYRS)
C	IREC(NLmax)	integer		record number on file IYRF
C	DOYS8(NLmax)	real*8		time of observation: day of year (incl. fraction for time of day)
C	XDS (NLmax)	real		Sun-Earth distance
C	XLS (NLmax)	real		geocentric ecliptic longitude Sun
C	XLL (NLmax)	real		geocentric ecliptic lng(P)-lng(Sun) (deg)
C	XDL (NLmax)	real		geocentric ecliptic lat(P) (deg)
C	XCE (NLmax)	real		modified Carrington variable of sub-Earth point on Sun
C	XE  (NLmax)	real		elongation (deg) (>0: East of Sun; <0: West of Sun)
C	XC  (NLmax)	real		modified Carrington variable point P after traceback to
C					  heliocentric distance Radius at speed VV
C	YL(NLmax)	real		heliographic latitude point P
C	VV(NLmax)	real		IPS velocity (km/s)
C	GG(NLmax)	real		IPS disturbance factor
C	NSmax		integer		# data points saved in XCsav,YLsav,VVsav
C	SRCsav(NSmax)	character	scratch arrays (used internally only),
C	IYRFsav(NSmax)	integer		
C	IRECsav(NSmax)	integer
C	IYRSsav(NSmax)	integer
C	DOYSsav(NSmax)	real*8
C	XDSsav (NSmax)	real
C	XLLsav (NSmax)	real
C	XDLsav (NSmax)	real
C	XCEsav (NSmax)	real
C	XEsav  (NSmax)	real
C	XCsav  (NSmax)	real
C	YLsav  (NSmax)	real
C	VVsav  (NSmax)	real
C	GGsav  (NSmax)	real
C RESTRICTIONS:
C	Data files are read using a logical unit number assigned by iGetLun.
C	Since a data file may be left open on return, there is a potential
C	danger if the same unit is used somewhere. Assign unit numbers
C	with iGetLun to make sure this doesn't happen.
C FUNCTIONS/SUBROUTINES:
C	bOpenFile, itrim, SunNewcomb, iSetFileSpec, iGetFileSpec, Julian
C EXTERNAL:
C INCLUDE:
C	include		'filparts.h'
C	include		'sun.h'
C COMMON BLOCKS:
C SIDE EFFECTS:
C PROCEDURE:
C >	All points in the range [XCtest1,XCtest2] are extracted from the data files
C >!!!!	The modified Carrington variables are in units offset by NCoff from
C	the regular Carrington rotation number (i.e. in units which can be
C	used as index to JDCar)
C >	Two external user-defined are needed to read the data files.
C	External functions exist for reading the Nagoya and the UCSD velocity
C	IPS data: iReadNagoyan8, iProcessNagoyan8
C		  iReadUCSDn8  , iProcessUCSDn8
C		  iReadOotyn8 , iProcessOotyn8
C	These are programmed in the form:
C
C	function iReadNagoyan8(iU,iFirst)
C	... ( read a record from the data file ) ...
C	return
C	entry iProcessNagoyan8(iRec,nCar,JDCar,NCoff,Radius,
C		cSRCVG,cSrc,nYr,Doy8,XCsc,XCobs,xLat,xVV,xGG,dRA,
C		rLngSun,rLngP,rLatP,rEloP )
C	... ( process the record to produce the output values	...
C	...   nYr,Doy8,XCsc,XCobs,xLat,xVV,xGG,dRA,		...
C	...   rLngSun,rLngP,rLatP,rEloP )
C	return
C	end
C
C	Input values to iReadVIPSn
C	iU	integer		logical unit number assigned by bOpenFile
C	iFirst	integer		1 : initializes record counter
C					(used when reading 1st record of file)
C				0 : adds one to record counter
C 	Output values of iProcessVIPSn:
C	iRec	integer		record # on file
C	cSrc	character	Source name
C	nYr	integer		year and ..
C	Doy	real		.. day of year of observation
C	XCsc	real		modified Carington variable for sub-SC point
C				  (NCoff already subtracted)
C	XCobs	real		modified Carrington variable of point P after
C				  trace-back to distance Radius (NCoff already subtracted)
C	xLat	real		heliographic latitude point P (deg)
C	xVV	real		IPS velocity (km/s)
C	xGG	real		IPS disturbance factor
C	dRA	real		RA(point P)-RA(Sun) (deg)
C	rLngSun real		geocentric ecliptic longitude Sun (deg)
C	rLngP	real		geocentric ecliptic longitude of point P (deg)
C					(relative to Sun)
C	rLatP	real		geocentric ecliptic latitude  of point P (deg)
C	rEloP	real		geocentric source elongation (deg)
C				NOTE: the sign is used to indicate East/West of Sun
C					rEloP>0: East of Sun; rEloP<0: West of Sun
C MODIFICATION HISTORY:
C	JUN-1994, Paul Hick (UCSD)
C	OCT-1996, Paul Hick (UCSD); added Ooty routines
C-
	subroutine ReadVIPSn8(iReadVIPSn,iProcessVIPSn,cWildVIPS,
     &		nCar,JDCar,bAuto,XCtest1,XCtest2,NCoff,MJDref,
     &		Radius,
     &		XElow,XEhigh,iEorW,XRlim,
     &		NLmax,NL,cvgfiles,SRCVG,SRC,SRCVGsav,SRCsav,
     &		IYRF,IREC,IYRS,DOYS8,XDS,XLS,XLL,XDL,XCE,XE,XC,YL,VV,GG,
     &		NSmax,IYRFsav,IRECsav,IYRSsav,DOYSsav8,XDSsav,XLSsav,XLLsav,
     &		XDLsav,XCEsav,XEsav,XCsav,YLsav,VVsav,GGsav)

	real*8		JDCar(nCar),
     &			MJDref,
     &			xTT,JEpoch,
     &			Doy8,
     &			DOYS8(NLmax),DOYSsav8(NSmax)		! Time: doy of year

	integer		IYRF(NLmax), IYRFsav(NSmax),
     &			IREC(NLmax), IRECsav(NSmax),
     &			IYRS(NLmax), IYRSsav(NSmax)

	real		XDS (NLmax), XDSsav (NSmax),		! Sun-Earth distance
     &			XLS (NLmax), XLSsav (NSmax),		! Ecliptic longitude Sun
     &			XLL (NLmax), XLLsav (NSmax),
     &			XDL (NLmax), XDLsav (NSmax),
     &			XCE (NLmax), XCEsav (NSmax),
     &			XE  (NLmax), XEsav  (NSmax),		! Elongation point P
     &			XC  (NLmax), XCsav  (NSmax),		! Carrington variable of Point P
     &			YL  (NLmax), YLsav  (NSmax),		! Heliographic latitude of Point P
     &			VV  (NLmax), VVsav  (NSmax),		! Solar wind velocity (km/s)
     &			GG  (NLmax), GGsav  (NSmax)		! G-factor

	character	cWildVIPS*(*),
     &			cvgfiles(*)*11,
     &			SRC(NLmax)*9,SRCsav(NSmax)*9,cSrc*9,SRCVG(NLmax)*115,SRCVGsav(NSmax)*115,cSRCVG*115,
     &			cSay*9			/'ReadVIPSn8'/,
     &			cFile*80,cFiles*11

	integer		iFrstYr			/0/,
     &			NS			/0/

	real		XCmargin		/1./

	logical		bAuto,
     &			bOpen			/.FALSE./,
     &			bGetRecord,
     &			bReading,
     &			bOpenFile,
     &			bMore,
     &			bReadErr,
     &			bSaveOverrun		/.FALSE./

	integer		nYrLen			/4/		! # digits used to identify year in file name
	integer		Str2Str, Flt2Str

	include		'dirspec.h'
	include		'filparts.h'
	include		'openfile.h'

	save		bOpen, bGetRecord, bSaveOverrun,
     &			iU, iYr, NS,
     &			iFrstYr, iFile

	logical		bLOSCheckLocn8
	integer		NLOS			/0/
	save		XCbeg,XCend,NLOS,NINSIDE,DLOS
C-------
C On first call check whether data files are present

	if (iFrstYr .eq. 0) then
	    I = iSetFileSpec(cWildVIPS)
	    I = iGetFileSpec(FIL__NAME,FIL__TYPE,cFile)

	    iFile = index(cFile,cWildChar(:nYrLen))
	    if (iFile .eq. 0) call Say(cSay,'E','INVTMPLTE',
     &			'wildcard template must contain '//cWildChar(:nYrLen))

	    if (iSearch(1,cWildVIPS,cFile) .ne. 1)
     &		call Say(cSay,'E','FILNOTFOUND',
C     &			'no files matching wildcard found :'//cWildVIPS)
     &			'no files matching wildcard found :')

	    I = iSetFileSpec(cFile)
	    I = iGetFileSpec(FIL__NAME,FIL__TYPE,cFile)

	    I = 1
	    call Str2Flt_Int(.TRUE.)			! Save year of first file
	    call Str2Flt(cFile(iFile:iFile+nYrLen-1),I,iFrstYr)
	    iFile = index(cWildVIPS,cWildChar(:nYrLen)) ! Save position of wildcard %%%% in cWildVIPS
	end if
C------
C Check the save scratch arrays for valid data points.
C NOTE: NS can only be unequal zero in AUTO mode. In AUTO mode only XCtest1
C	and XCtest2 change (both increase). All other data selection
C	criteria remain the same.

	NL = 0				! Initialize data point counter
	LOS = 0				! Counter for rejected lines of sight by bLOSCheckLocn8

	if (NS .ne. 0) then		! Only in AUTO mode
C-------
C Scan the `save' arrays. Only points with XC >= XCtest1 are needed for
C subsequent maps. Discard points with XC < XCtest1.

	    N = 0
	    do I=1,NS
		if (XCsav(I) .ge. XCtest1) then
		    N = N+1
		    SRCVGsav(N)= SRCVGsav(I)
		    SRCsav(N)  = SRCsav(I)
		    IYRFsav(N) = IYRFsav(I)
		    IRECsav(N) = IRECsav(I)

		    IYRSsav(N) = IYRSsav(I)
		    DOYSsav8(N)= DOYSsav8(I)
 		    XDSsav (N) = XDSsav (I)
		    XLSsav (N) = XLSsav (I)
		    XLLsav (N) = XLLsav (I)
		    XDLsav (N) = XDLsav (I)
		    XCEsav (N) = XCEsav (I)
		    XEsav  (N) = XEsav  (I)
		    XCsav  (N) = XCsav  (I)
		    YLsav  (N) = YLsav  (I)
		    VVsav  (N) = VVsav  (I)
		    GGsav  (N) = GGsav  (I)
		end if
	    end do
	    NS = N			! Valid # point in save arrays

C--------
C Pick up all points in the save arrays which are needed for the current
C map. Since all points with XC < XCtest1 have just been discarded, the
C points are only tested for XC <= XCtest2.

	    do I=1,NS
		if (XCsav(I) .le. XCtest2) then
		    NL = NL+1
	            IL1 = iSetFileSpec(cFile)
	            IL1 = iGetFileSpec(FIL__NAME,FIL__TYPE,cFiles)
	    	    cvgfiles(NL) = cFiles
		    SRCVG(NL)= SRCVGsav(I)
		    SRC(NL)  = SRCsav(I)
		    IYRF(NL) = IYRFsav(I)
		    IREC(NL) = IRECsav(I)

		    IYRS(NL) = IYRSsav(I)
		    DOYS8(NL)= DOYSsav8(I)
 		    XDS (NL) = XDSsav (I)
 		    XLS (NL) = XLSsav (I)
 		    XLL (NL) = XLLsav (I)
 		    XDL (NL) = XDLsav (I)
 		    XCE (NL) = XCEsav (I)
 		    XE  (NL) = XEsav  (I)

		    XC  (NL) = XCsav  (I)
		    YL  (NL) = YLsav  (I)
		    VV  (NL) = VVsav  (I)
		    GG  (NL) = GGsav  (I)
		endif
	    end do
	end if

C-------
C If no file open then check which file is to be read.
C (bOpen=.TRUE. only happens in AUTO mode. Then also bGetRecord=.FALSE.,
C indicating an unprocessed record is still in memory).
C The data files are assumed to be organized in yearly files: one file
C for each year of data. XCtest1-XCmargin is used to determine which yearly
C file must be opened first.

	if (.not. bOpen) then

	    XCt = XCtest1-XCmargin
	    I = XCt
	    if (I .lt. 1 .or. I .ge. nCar) call Say(cSay,'E','InvXCtest1',
     &		    'start value outside range covered by JDCar')

	    xTT = JDCar(I)+(XCt-I)*(JDCar(I+1)-JDCar(I))
	    call Julian8(1,iYr,Doy8,xTT,JEpoch)
	    iYr = max(iYr,iFrstYr)	! Year of file to be opened

	    iYr = iYr-1			! Will be incremented again in a second
	    bGetRecord = .TRUE.		! Make sure to read record

	end if
C-------
C  Read data until
C  -  end point XCtest2 is reached, or
C  -  EOF is reached, or
C  -  arrays are full
C
C bOpen and bGetRecord are used in subsequent calls.
C If a record has been read succesfully, but is rejected because it's out
C of range (iProc=2) or because the arrays are full (NL>=NLmax), then
C the file is kept open (bOpen=.TRUE.) and bGetRecord is set .FALSE.
C On the next call to ReadVIPSn8, the procedure will skip straight to the
C iProcessVIPS call to try processing the last record read in the previous
C call before reading the next record from the file that is still open.

	bMore	 = .TRUE.
	bReading = .FALSE.		! Only used to control messages to screen

	do while (bMore)

	    if (.not. bOpen) then	! Open next year's file
		iYr = iYr+1
		iFirst = 1		! Force ReadVIPSn8 to initialize record counter
C-------
C Complete the data file name by changing the file extension to the
C proper year. Open data file. If bOpenFile returns 0 the file probably
C didn't exist (exit do while loop). bOpenFile assigns iU if successfully
C opening the file.

		cFile = cWildVIPS	! Replace wildcard by year
		I = Int2Str(iYr,cFile(iFile:iFile+nYrLen-1))
		
		bOpen = bOpenFile(OPN__REOPEN+OPN__TEXT+OPN__READONLY,iU,cFile,iRecl)
		bMore = bOpen		! Open error will exit while loop
	        IL1 = iSetFileSpec(cFile)
	        IL1 = iGetFileSpec(FIL__NAME,FIL__TYPE,cFiles)
	    end if

	    bReadErr = .FALSE.
	    if (bOpen .and. bGetRecord) then
		if (.not. bReading) then
		    call Say(cSay,'I','Reading',cFile)
		    bReading = .TRUE.
		end if

		bMore = iReadVIPSn(iU,iFirst) .eq. 0
		if (bMore) then
		    iFirst = 0		! Record counter has been initialized
		else			! Read error: EOF assumed
		    iU = iFreeLun(iU)
		    bOpen = .FALSE.
		    bReadErr = .TRUE.
		end if

	    end if

	    if (bReadErr) then		! After read error, open new file
		bMore = .TRUE.

	    else if (bMore) then
		iProc = iProcessVIPSn(iR,nCar,JDCar,NCoff,XCtest2,Radius,
     &			cSRCVG,cSrc,nYr,Doy8,XCsc,XCobs,xLat,xVV,xGG,dRA,
     &			rLngSun,rDisSun,rLngP,rLatP,rEloP)

		bMore = iProc .ne. 2		! Check for XCsc past XCtest2
		if (iProc .eq. 1) then		! Record contains valid data with XCsc <= XCtest2

		    call Julian8(10,nYr,Doy8,xTT,JEpoch)
						! All XE (East+West), only positive XE (East) or only negative XE (West)
		    Rdum = (1-abs(iEorW))*abs(rEloP)+iEorW*rEloP
     		    if ((XElow .le. Rdum .and. Rdum .le. XEhigh) .and.
     &			( abs(dRA) .le. XRlim ) .and.
     &			(MJDref .eq. 0 .or. xTT .le. MJDref)
     &		       ) then

			if (XCtest1 .le. XCobs) then
			    if (bAuto) then
			        bSaveOverrun = NS .ge. NSmax
				if (.not. bSaveOverrun) then

				    NS = NS+1
				    SRCVGsav(NS)= cSRCVG
				    SRCsav(NS)	= cSrc
				    IYRFsav(NS) = iYr
				    IRECsav(NS) = iR

				    IYRSsav(NS) = nYr
				    DOYSsav8(NS)= Doy8
				    XDSsav (NS) = rDisSun
				    XLSsav (NS) = rLngSun
				    XLLsav (NS) = rLngP
				    XDLsav (NS) = rLatP
				    XCEsav (NS) = XCsc
				    XEsav  (NS) = rEloP

				    XCsav  (NS) = XCobs
				    YLsav  (NS) = xLat
				    VVsav  (NS) = xVV
				    GGsav  (NS) = xGG
				end if

			    end if
C--------
C Pick up all points needed for the current map.Test only for XC <= XCtest2

			    if (XCobs .le. XCtest2) then
C--------
C Check what part of line of sight is inside the actual boundaries [XCbeg,XCend]

				if (bLOSCheckLocn8(XCbeg,XCend,NLOS,NINSIDE,DLOS,
     &						nYr,Doy8,rLngP,rLatP,XCsc,rLngSun,rDisSun)) then

				    if (NL .lt. NLmax) then

					NL = NL+1
	    	    			cvgfiles(NL) = cFiles
				        SRCVG(NL)= cSRCVG
					SRC(NL)	 = cSrc
					IYRF(NL) = iYr
					IREC(NL) = iR

					IYRS(NL) = nYr		! Time of observation
					DOYS8(NL)= Doy8
					XDS (NL) = rDisSun	! Sun-Earth distance (AU)
					XLS (NL) = rLngSun	! Geocentric ecliptic longitude of Sun (deg)
					XLL (NL) = rLngP	! Geocentric ecliptic lng(P)-lng(Sun)
					XDL (NL) = rLatP	! Geocentric ecliptic lat(P)
					XCE (NL) = XCsc		! Carrington variable sub-Earth point
					XE  (NL) = rEloP

					XC  (NL) = XCobs
					YL  (NL) = xLat
					VV  (NL) = xVV
					GG  (NL) = xGG
c	write (*,'(4F10.3)') XC(NL),YL(NL),VV(NL),GG(NL),TTsav(I)
c	write (*,'(I4,F6.2,4F10.2)') ID,UT,XCtest1,XCobs,XCtest2
				    else
					call Say(cSay,'I','NLmax','Data arrays full: not all data processed')
					bMore = .FALSE.
				    end if
				else
				    LOS = LOS+1

				end if
			    end if
			end if
		    end if
		end if

C-------
C At this point bMore = .FALSE. if iProc=2 (XCsc past XCtest2) or the data arrays are full

		bGetRecord = bMore .or. .not. bGetRecord

	    end if

	end do

	if (.not. bAuto .and. bOpen) then	! If not in AUTO mode, make sure to close open file
	    iU = iFreeLun(iU)
	    bOpen = .FALSE.
	end if

	if (bSaveOverrun) then			! Can only happen in AUTO mode
	    call Say(cSay,'W','NSmax','Save arrays overrun: contents discarded')
	    NS = 0
	    bSaveOverrun = .FALSE.
	end if

	I = 0
	I = I+Int2Str(NL,cFile(I+1:))+1
	I = I+Str2Str('sources in',cFile(I+1:))+1
	I = I+Flt2Str(NCoff+XCtest1,2,cFile(I+1:))
	I = I+Str2Str('-',cFile(I+1:))
	I = I+Flt2Str(NCoff+XCtest2,2,cFile(I+1:))

	if (NLOS .ne. 0) then
	    I = I+2
	    I = I+Int2Str(LOS,cFile(I+1:))+1
	    I = I+Str2Str('; sources rejected',cFile(I+1:))
	end if
	call Say(cSay,'I','# obs',cFile)

	return

	entry ReadVipsLOSCheckn(XCbegIn,XCendIn,NLOSIn,NINSIDEIn,DLOSIn)

	XCbeg 	= XCbegIn
	XCend 	= XCendIn
	NLOS  	= NLOSIn
	NINSIDE = NINSIDEIn
	DLOS  	= DLOSIn

	return
	end

	logical function bLOSCheckLocn8(XCbeg,XCend,NLOS,NINSIDE,DLOS,iYr,Doy8,rLngLOS,rLatLOS,ObsXC,rLngSun,rDisSun)

	real*8		Doy8
	
	bLOSCheckLocn8 = NLOS .eq. 0
	if (bLOSCheckLocn8) return

	ObsLng = rLngSun-180.

	Inside = 0
	do J=1,NLOS			! Loop over all segments along LOS
	    P1 = rLngLOS		! Topocentric ecliptic longitude diff. with Earth-Sun line
	    P2 = rLatLOS		! Topocentric ecliptic latitude
	    P3 = (J-.5)*DLOS		! Distance along LOS (units of ObsR)
					! Convert to heliocentric coordinates
	    call POINT_ON_LOS(P1,P2,P3,ELO,ELOS,iEorW)
	    P3 = P3*rDisSun		! Convert to AU

	    P1 = ObsLng+P1		! Heliocentric longitude of P
	    call ECLIPTIC_HELIOGRAPHIC8(0,iYr,Doy8,P1,P2)
					! Convert to heliographic coordinates
	    P1 = XMAP_OBS_POS(ObsXC,P1)	! Modified Carrington variable of P

	    if (XCbeg .le. P1 .and. P1 .le. XCend) Inside = Inside+1
	end do

	bLOSCheckLocn8 = Inside .ge. NINSIDE

	return
	end
C+
C >	Structure of Nagoya files:
C	cSrce		string		descriptor for radio source
C	iY,iM,iD	integer 	year, month and day of observation
C	UT		real		universal time of day
C	DP		real		heliocentric distance point P (AU)
C	iTHR		integer		?? (only from 1990 - 1993)
C	iHLA		integer		heliocentric ecliptic latitude point P
C	iHLO		integer 	heliocentric longitude point P relative to
C					Sun-Earth direction (i.e. a difference of
C					ecliptic longitudes)
C	iGLA,iGLO	integer		heliographic location point P (degrees)
C	iROT		integer		Carrington rotation number
C	iVEL		integer		IPS velocity (km/s)
C	iDVEL		integer		error in velocity (km/s)
C	SCINDX/G-VALUE	real		scintillation index (arbitrary units);
C					present from 1990 to 1996) or G-value (after 1997)
C Record length of original Nagoya files:
C	1990-1993: 73 chars
C	1994-1996: 68 chars (no THR column)
C	1997-1998: 66 chars (no THR, g-value instead of source index
C In the processed files, iTHR has been omitted, and the source index has been replaced by
C	zero (i.e. results in a zero g-value)
C-
	function iReadNagoyan8(iU,iFirst)
	character	cMon*3,
     &			cSrce*9,cSrc*9,cSRCVGe*66,cSRCVG*66
	real*8		JDCar(nCar),
     &			dLngSun, dLatSun, dDisSun
        real*8		Doy8
	character	cFmt*40		/'(A9,3I2,F6.2,F5.2,I4,I5,I4,3I5,I4,F8.5)'/

	external	EARTH8
	include		'sun.h'

	save		cSrce,iY,iM,iD,UT,DP,iHLA,iHLO,iGLA,iGLO,iROT,iVEL,iDVEL,GVAL,iRec,cSRCVGe

	iU = iU					! Avoids compiler warning on VMS
	read (iU,cFmt,iostat=iReadNagoyan8) cSrce,iY,iM,iD,UT,DP,iHLA,iHLO,iGLA,iGLO,iROT,iVEL,iDVEL,GVAL

	write(cSRCVGe(1:66),cFmt)
     &		cSrce,iY,iM,iD,UT,DP,iHLA,iHLO,iGLA,iGLO,iROT,iVEL,iDVEL,GVAL

	if(iY.lt.50) iY = iY + 100
	iRec = (1-iFirst)*iRec+1		! iFirst=1 initializes record counter

	return

	entry iProcessNagoyan8(iR,nCar,JDCar,NCoff,XCend,Radius,
     &		cSRCVG,cSrc,nYr,Doy8,XCsc,XCobs,xLat,xVV,xGG,dRA,
     &		rLngSun,rDisSun,rLngP,rLatP,rEloP)

	if (iVEL .eq. 0) then			! No velocity available: reject
	    iProcessNagoyan8 = 0
	    return
	end if

	iR   = iRec

	cSRCVG = cSRCVGe
	cSrc = cSrce
	nYr  = 1900+iY
	cMon = '   '				! Make sure iM is used
	call DATE_DOY(0,nYr,cMon,iM,iD,iDoy)
	Doy8  = iDoy+dble(UT/24.)

 	XCsc = XMAP_SC_POS8(EARTH8,nYr,Doy8,-nCar,JDCar)

	if (XCsc .eq. BadR4()) then		! Time outside JDCar range
		iProcessNagoyan8 = 0
		return
	end if
						! Sub-Earth position
	if (XCsc .gt. XCend) then		! Sub-Earth past XCend
	    iProcessNagoyan8 = 2
	    return
	end if

	XCobs = iROT+(360.-iGLO)/360.		! Point-P at DP
	XCobs = XCobs+Sun__Spiral*(Radius-DP)/iVEL ! Point-P traced to Radius
	XCobs = XCobs-NCoff			! Subtract offset

	xLat = iGLA				! Heliographic lat(P)
	xVV  = iVEL				! IPS velocity
	xGG  = GVAL				! Set to 0 if GVAL not available

	rLng = iHLO				! Heliocentric ecliptic lng(P)-lng(Earth)
	rLat = iHLA				! Heliocentric ecliptic lat(P)
	rP   = DP				! Distance Sun-P (AU)
	call POINT_ON_LOS(rLng,rLat,rP,rEloSun,rEloP,i)
c	call PointOnLos(rLng,rLat,rP,rEloSun,rEloP,i)

	rLngP = rLng				! Geocentric ecliptic lng(P)-lng(Sun)
	rLatP = rLat				! Geocentric ecliptic lat(P)

	rEloP = -i*rEloP			! rEloP>0 : East of Sun
						!      <0 : West of Sun
	call SunNewcomb8(0,nYr,Doy8,dLngSun,dLatSun,dDisSun)
	rDisSun = dDisSun
	rLngSun = dLngSun			! Geocentric ecliptic lng(Sun)
	rSun    = rLngSun

	rLng = rSun+rLng			! Geocentric ecliptic lng(P)
	Doy  = sngl(Doy8)
	call ECLIPTIC_EQUATOR(0,nYr,Doy,rLng,rLat) ! rLng = RA of P

	rLat = dLatSun				! Geocentric ecliptic lat(Sun)
	call ECLIPTIC_EQUATOR(0,nYr,Doy,rSun,rLat) ! rSun = RA of Sun

	dRA = rLng-rSun				! RA(P)-RA(Sun)
	if (abs(dRA) .gt. 180.) dRA = dRA-sign(1.,dRA)*360.

	iProcessNagoyan8 = 1

	return
	end
C+
C NAME:
C	iReadOotyn8
C PROCEDURE:
C >	Structure of Ooty files:
C	iY,iM,iD	integer 	year, month and day of observation
C	IST		integer		time of observation in IST - hr and min
C					(the longitude difference between
C					Greenwich and the IST longitude is 82.5
C					deg corresponding to 5.5 hrs time difference. 
C	UT		real		universal time (should be IST-5.5)
C	cSrce		string		descriptor for radio source
C	iRA50(3)	integer		1950 Right Ascension (h,m,s)
C	iDEC50(3)	integer		1950 declination (deg,m,s)
C	iRA(3)		integer		Right Ascension at equinox of date
C	iDEC(3)		integer		declination at equinox of date
C	ELO		real		source elongation (deg)
C	SIZE		real		angular size of the compact component
C					of the source (seconds of arc)
C	iVEL		integer		IPS velocity (km/s)
C	AR		real		axial ratio of the solar wind density
C					irregularities
C	SN		real		signal to noise ratio of the spectrum
C	SCINDX		real		scintillation index
C	GVAL		real		disturbance factor
C
C From an email by Janardhan Padmanabhan:
C
C I have selected only those sources in this list with elongations
C greater than or equal to 10 deg (the turnover in the m vs epsilon 
C curve at 327 MHz is around 10 deg) and signal to noise greater than or 
C equal to 10db. However, please note that source sizes derived for 
C sources with S/N less than 15 db will have larger errors.  Therefore 
C if source sizes are an important input in your program, then use 
C only data with S/N greater than or equal to 15 db.
C
C NOTE: 
C for 8 observations in the file I was unable to extract the scintillation 
C indicies. 7 of these are the august observations. I've marked these in the
C list with a "N/A" (not available) in the last column. 
C-
	function iReadOotyn8(iU,iFirst)
	character	cMon*3,
     &			cSrce*8,cSrc*8,cSRCVGe*105,cSRCVG*105
	real*8		JDCar(nCar),
     &			dLngSun, dLatSun, dDisSun
     	real*8		Doy8,Doy88
	integer		iRA50 (3), iRA (3),
     &			iDEC50(3), iDEC(3)
	external	EARTH8
	save		iY,iMO,iD,UT,cSrce,iRA50,iDEC50,iRA,
     &			iDEC,ELO,SIZE,iVEL,SN,SCINDX,GVAL,iRec,cSRCVGe
	include		'sun.h'

	iU = iU
C	read (iU,'(2X,3I2,2X,2I2,F8.3,2X,A8,1X,2(1X,6I3),F7.2,F7.3,I6,F6.1,2F7.3)',
	read (iU,'(2X,3I2,F8.3,2X,A8,1X,2(1X,6I3),F7.2,F7.3,I6,F6.1,2F7.3)',
     &		iostat=iReadOotyn8) iY,iMO,iD,UT,cSrce,iRA50,iDEC50,iRA,
     &		iDEC,ELO,SIZE,iVEL,SN,SCINDX,GVAL

C	The following converts back to ooty time

C	cMon = '   '
C	nnYr  = 2000+iY
C	call DATE_DOY(0,nnYr,cMon,iMO,iD,iDoy)			! Convert to iDOY
C	Doy88 = dble(float(iDoy)) + dble(UT/24.) + dble(5.5/24.) 	! Convert to Ooty DOY 
C	iDoy = int(sngl(Doy88))
C	fract = sngl(Doy88 - dble(float(iDoy)))
C	iH  = int(fract*24.0)
C	iMI = nint((fract*24.0 - float(iH))*60.)
C	call DATE_DOY(1,nnYr,cMon,iMO,iD,iDoy)			! Convert back to Ooty iMO, and ID from iDoy
C	iY  = nnYr - 2000
	
C	write(*,'(2X,3I2,2X,2I2,F8.3,2X,A8,1X,2(1X,6I3),F7.2,F7.3,I6,F6.1,2F7.3)')
C     &		iY,iMO,iD,iH,iMI,UT,cSrce,iRA50,iDEC50,iRA,iDEC,ELO,SIZE,iVEL,SN,SCINDX,GVAL

C	End convert back to Ooty time

C	write(cSRCVGe(1:111),'(2X,3I2,2X,2I2,F8.3,2X,A8,1X,2(1X,6I3),F7.2,F7.3,I6,F6.1,2F7.3)')
C     &		iY,iMO,iD,iH,iMI,UT,cSrce,iRA50,iDEC50,iRA,iDEC,ELO,SIZE,iVEL,SN,SCINDX,GVAL

	write(cSRCVGe(1:105),'(2X,3I2,F8.3,2X,A8,1X,2(1X,6I3),F7.2,F7.3,I6,F6.1,2F7.3)')	! save file in UT as read
     &		iY,iMO,iD,UT,cSrce,iRA50,iDEC50,iRA,iDEC,ELO,SIZE,iVEL,SN,SCINDX,GVAL

C	print *, 'cSRCVGe',cSRCVGe
C	if (iY .lt. 50) iY = iY+100

	iRec = (1-iFirst)*iRec+1		! iFirst=1 initializes record counter

	return
C+
C NAME:
C	iProcessOotyn8
C CALLING SEQUENCE:
	entry iProcessOotyn8(iR,nCar,JDCar,NCoff,XCend,Radius,
     &		cSRCVG,cSrc,nYr,Doy8,XCsc,XCobs,xLat,xVV,xGG,dRA,
     &		rLngSun,rDisSun,rLngP,rLatP,rEloP)
C-
	if (iVEL .le. 0) then			! No velocity available: reject
	    iProcessOotyn8 = 0			! (never happens?)
	    return
	end if

	cSRCVG = cSRCVGe
	cSrc = cSrce
	nYr  = 2000+iY				! error was 1900 BVJ 08/24/08
	cMon = '   '				! Make sure iMO is used
	call DATE_DOY(0,nYr,cMon,iMO,iD,iDoy)
	Doy8 = dble(float(iDoy))+dble(UT/24.)		! to process the most recent Ooty data (08/22/08)
C	Doy8 = dble(float(iDoy)) + dble(float(iH)/24.0 + float(iMI)/1440. - 5.5/24.) ! Installed 08/04/2008 BVJ
	XCsc = XMAP_SC_POS8(EARTH8,nYr,Doy8,-nCar,JDCar)

	if (XCsc .eq. BadR4()) then		! Time outside JDCar range
		iProcessOotyn8 = 0
		return
	end if
						! Sub-Earth position
	if (XCsc .gt. XCend) then		! Sub-Earth past XCend
	    iProcessOotyn8 = 2
	    return
	end if

	iR   = iRec

	call SunNewcomb8(0,nYr,Doy8,dLngSun,dLatSun,dDisSun)
						! Geocentric eclip. coord(Sun)
	rDisSun = dDisSun
	rLngSun = dLngSun			! Geocentric eclip. lng(Sun)

	rLng = dLngSun
	rLat = dLatSun				! Geocentric ecliptic lat(Sun)
	call ECLIPTIC_EQUATOR8(0,nYr,Doy8,rLng,rLat)
	dRA = rLng				! RA (Sun)

	rLng  = (iRA(1)+(iRA(2)+iRA(3)/60.)/60.)*15.	! RA (P) (deg)
	rLat  = iDEC(1)+(iDEC(2)+iDEC(3)/60.)/60.	! Decl (P)

	dRA = rLng-dRA				! RA(P)-RA(Sun)
	if (abs(dRA) .gt. 180.) dRA = dRA-sign(1.,dRA)*360.

	call ECLIPTIC_EQUATOR8(1,nYr,Doy8,rLng,rLat)
						! Geocentric eclip. coord(P)
	rLng  = rLng-rLngSun
	rLngP = rLng				! Geocentric ecliptic lng(P)-lng(Sun)
	rLatP = rLat				! Geocentric ecliptic lat(P)

	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,nYr,Doy8,rLng,rLat)
						! Heliographic coord(P)
	xLat = rLat				! Heliographic lat(P)

	XCobs = XMAP_OBS_POS(XCsc,rLng)		! Point-P at rP
	XCobs = XCobs+(Radius-rP)/iVEL*SUN__SPIRAL ! Point-P traced to Radius

	xVV = iVEL				! IPS velocity
	xGG = GVAL

	iProcessOotyn8 = 1

	return
	end
C+
C >	Structure of UCSD files:
C	iY		integer		year-1900
C	IDOY		integer		day of year
C	IH,IM		integer		time of day
C	cSrce		string		descriptor for radio source
C	IQUAL		integer		quality indicator
C	VEL		real		IPS velocity (km/s)
C	DVEL		real		error in velocity (km/s)
C	VELpeak		real		IPS peak velocity (km/s)
C	DVELpeak	real		error in peak velocity (km/s)
C	SCALE		real		spatial scale of IPS pattern (km)
C	ISITE		integer		site parameter
C	ELO		real		Sun-Source angle (elongation; deg)
C	ILONG		integer		Carrington longitude point P (mapped to Sun)
C	XC		real		modified Carrington variable point P (x100) (mapped to Sun)
C	HLAT		real		heliographic latitude point P (deg)
C	ILONGe		real		Carrington longitude sub-Earth point
C	XCe		real		modified Carrington variable sub-Earth point (x100)
C
C NOTES:
C	ILONG = (1-(XC-int(XC)))*360.
C	XCe = XMAP_SC_POS8(EARTH8,nYr,Doy8,nCar,JDCar)
C	Both appear to be true to within 2 degree in longitude, or so.
C
C       For large elongation the point P was probably taken at a fixed distance
C       from the Earth. In principle, the distance can be reconstructed from
C       ILONG and HLAT by reversing the traceback to the solar surface.
C       I tried to do this for a number of elongations above 90 degrees and
C       found point-P distances close to 0.3 AU:
C	Elongation    Dist (AU)			Elongation    Dist (AU)
C	91.00000      0.3017747			115.0000      0.3055834
C	92.00000      0.3120534			116.0000      0.2949757
C	93.00000      0.3077490			117.0000      0.2932528
C	94.00000      0.3076252			121.0000      0.2888033
C	95.00000      0.2509942			122.0000      0.2779170
C	96.00000      0.2847953			123.0000      0.2820099
C	97.00000      0.2842951			124.0000      0.2752993
C	98.00000      0.2923858			125.0000      0.2946435
C	99.00000      0.3047100			126.0000      0.2908247
C	100.0000      0.2982190			127.0000      0.2950612
C	101.0000      0.2989129			128.0000      0.2925285
C	102.0000      0.3032739			129.0000      0.2909675
C	103.0000      0.2805566			130.0000      0.2826101
C	104.0000      0.2944881			131.0000      0.2786521
C	105.0000      0.2971804			132.0000      0.3057403
C	106.0000      0.2883710			133.0000      0.2940556
C	107.0000      0.3060628			134.0000      0.2966698
C	108.0000      0.2852099			135.0000      0.2718044
C	110.0000      0.2984203			136.0000      0.272730
C	111.0000      0.2984303			137.0000      0.276677
C	112.0000      0.2980161			140.0000      0.2886397
C	113.0000      0.2884684
C	114.0000      0.2940932
C-
	function iReadUCSDn8(iU,iFirst)
	character	cSrce*9,cSrc*9,cSRCVGe*115,cSRCVG*115
	real*8		JDCar(nCar),
     &			dLngSun, dLatSun, dDisSun
     	real*8		Doy8
	external	EARTH8
	save		iY,iD,iH,iM,cSrce,iQUAL,VEL,DVEL,R1,R2,R3,I1,
     &			ELO,ILON,XCobs_in,HLAT,RLONsc,XCsc_in,iRec,cSRCVGe
	include		'sun.h'

	iU = iU
	read (iU,'(2I4,I3,I2,A,I3,5F9.1,I2,F7.1,I4,F9.1,2F7.1,F9.1)',iostat=iReadUCSDn8)
     &		iY,iD,iH,iM,cSrce,iQUAL,VEL,DVEL,R1,R2,R3,I1,
     &		ELO,ILON,XCobs_in,HLAT,RLONsc,XCsc_in

	write(cSRCVGe(1:115),'(2I4,I3,I2,A,I3,5F9.1,I2,F7.1,I4,F9.1,2F7.1,F9.1)')
     &		iY,iD,iH,iM,cSrce,iQUAL,VEL,DVEL,R1,R2,R3,I1,
     &		ELO,ILON,XCobs_in,HLAT,RLONsc,XCsc_in

	XCobs_in = XCobs_in/100.

	iRec = (1-iFirst)*iRec+1		! iFirst=1 initializes record counter

	return

	entry iProcessUCSDn8(iR,nCar,JDCar,NCoff,XCend,Radius,
     &		cSRCVG,cSrc,nYr,Doy8,XCsc,XCobs,xLat,xVV,xGG,dRA,
     &		rLngSun,rDisSun,rLngP,rLatP,rEloP)

	if (iQUAL .lt. 21) then			! Data not good enough: reject
	    iProcessUCSDn8 = 0
	    return
	end if

	iR   = iRec

	cSRCVG = cSRCVGe
	cSrc = cSrce
	nYr  = 1900+iY
	Doy8  = iD+dble((iH+iM/60.)/24.)
	XCsc = XMAP_SC_POS8(EARTH8,nYr,Doy8,-nCar,JDCar) ! Sub-Earth position

	if (XCsc .eq. BadR4()) then		! Time outside JDCar range
		iProcessUCSDn8 = 0
		return
	end if

	if (XCsc .gt. XCend) then		! Sub-Earth past XCend
	    iProcessUCSDn8 = 2
	    return
	end if

	XCobs = XCobs_in+(Radius-Sun__RAu)/VEL*Sun__Spiral-NCoff
						! Point-P at Radius
	xLat = HLAT				! Heliographic lat(P)
	xVV  = VEL				! IPS velocity (km/s)
	xGG  = 0.				! No disturbance factor

	call SunNewcomb8(0,nYr,Doy8,dLngSun,dLatSun,dDisSun)
	rDisSun = dDisSun
	rLngSun = dLngSun			! Geocentric ecliptic lng(Sun)
	rSun    = rLngSun
C-------
C POINT_ON_LOS calculates the elongation based on the DP value assumed here.
C The cut-off at 0.30 appears to be the best choice and reproduces the
C elongations on the input file to within a few degrees for most sources.

	DP = max(0.30,cosd(ELO))		! Geocentric distance P (at least 0.3 AU)

	DP = 1.+DP*DP-2*DP*cosd(ELO)		! Cos-rule in Earth-Sun-P
	DP = sqrt(max(0.,DP))			! Distance Sun-P (AU)

	rLng = XCobs_in+(DP-Sun__RAu)/VEL*Sun__Spiral-NCoff
	rLng = rLng-int(rLng)			! Fractional Carrington rotation
	rLng = 360.*(1.-rLng)			! Heliographic lng(P)
	rLat = HLAT				! Heliographic lat(P)
	call ECLIPTIC_HELIOGRAPHIC8(1,nYr,Doy8,rLng,rLat)
						! rLng = Heliocentric ecliptic lng(P)
						! rLat = Heliocentric ecliptic lat(P)
	rLng = rLng-(180.+rSun)			! Heliocentric ecliptic lng(P)-lng(Earth)
	rP   = DP				! Distance Sun-P (AU)
	call POINT_ON_LOS(rLng,rLat,rP,rEloSun,rElo,i)

c	print *, ELO,rElo,ELO-rElo

	rLngP = rLng				! Geocentric ecliptic lng(P)-lng(Sun)
	rLatP = rLat				! Geocentric ecliptic lat(P)
	rEloP = -i*ELO				! rEloP>0 : East of Sun
						!      <0 : West of Sun
	rLng = rSun+rLng			! Geocentric ecliptic lng(P)

	Doy	= sngl(Doy8)
	call ECLIPTIC_EQUATOR(0,nYr,Doy,rLng,rLat) ! rLng = RA of P

	rLat = dLatSun				! Geocentric ecliptic lat(Sun)
	call ECLIPTIC_EQUATOR(0,nYr,Doy,rSun,rLat) ! rSun = RA of Sun

	dRA = rLng-rSun				! RA(P)-RA(Sun)
	if (abs(dRA) .gt. 180.) dRA = dRA-sign(1.,dRA)*360.

	iProcessUCSDn8 = 1

	return
	end