PRO fan_read, file=FILE, clear=CLEAR

;+
; NAME:
;	fan_read
; PURPOSE:
;	Read ASCII file containing data for time series
; CATEGORY:
;	I/O
; CALLING SEQUENCE:
;	fan_read, file=FILE, clear=CLEAR
; CALLS:
;	flt_string
; OPTIONAL INPUT PARAMETERS:
;	FILE	Name of ASCII file with time series data (if omitted user is prompted)
;	CLEAR	If set, all previously read data are dumped and only the new
;		data file remains in memory)
; OUTPUTS:
;	Through COMMON blocks
; INCLUDE:
	@compile_opt.pro		; On error, return to caller
; COMMON BLOCKS:
	common	CFANCY,		IPMAX,IP,BTYP,BLOW,BUP
	common	TFANCY,		YEAR0,DOY0,TUNIT,T0,F0,T1,F1,T2,F2,T3,F3,T4,F4,T5,F5
; RESTRICTIONS:
;	Times with the value -1 will be ignored
; PROCEDURE:
;	The input file should have the following layout:
;	- Header containing time origin (Year,doy)
;	- Subheader containing nr of rows NROW and columns NCOL of data
;	- NROW records, each containing NCOL numbers
;	The same groupings may be repeated on the same file.
;	The odd columns should contains the times; the even columns the
;	function values (e.g. intensities).
; MODIFICATION HISTORY:
;	MAY-1991, Paul Hick (ARC)
;-

IF n_elements(IP) EQ 0 OR keyword_set(CLEAR) THEN IP = -1

IF IP EQ 5 THEN BEGIN echo, 'ALL DATA ARRAYS ALREADY FULL'  &  RETURN  & ENDIF

IF NOT keyword_set(FILE) THEN BEGIN
	FILE = 'FANCY.TXT'
	echo, 'File name', FILE, /file
ENDIF

openr, LUN, /get_lun, FILE

WHILE NOT EOF(LUN) DO BEGIN
	ARR = ''								; Read first non-empty record of block of data
	REPEAT readf, LUN, ARR UNTIL strlen(strtrim(ARR)) GT 0

; Analyze first non-empty record

	IF strmid(ARR,0,4) EQ 'HXIS' THEN BEGIN
		ARR = strmid(ARR,5,strlen(ARR)-5)	; Drop 5-char start of record: 'HXIS '
		NEWYEAR = flt_string(ARR,strvec=A)		; Should have 2+3*NHX entries
		A = fix(A)
		NEWYEAR = A[0]						; 1st entry: year
		NEWDOY  = A[1]						; 2nd entry: doy
		OFFDOY  = 0.
		NHX = (n_elements(A)-2)/3			; # time series in file
		B1  = A[0*NHX+2:1*NHX+1]			; NHX entries
		B2  = A[1*NHX+2:2*NHX+1]			; NHX entries
		B3  = A[2*NHX+2:3*NHX+1]			; NHX entries
		NEWTUNIT = 'HRS'
	ENDIF ELSE BEGIN
		NHX = 0								; Not a HXIS file
		A = flt_string(ARR,strvec=A)
		NEWYEAR = fix(A[0])
		NEWDOY  = fix(A[1])
		IF n_elements(A) EQ 2 THEN	$		; A[1] is doy with fraction for time of day
			OFFDOY = A[1]-NEWDOY	$
		ELSE						$		; A[2:4] is hour, minutes and seconds
			OFFDOY = (A[2]+(A[3]+A[4]/60.)/60.)/24.
		ARR = strupcase(ARR)				; Find unit of time
		IF strpos(ARR,'DAY') NE -1 THEN NEWTUNIT = 'DAY' ELSE	$
		IF strpos(ARR,'HRS') NE -1 THEN NEWTUNIT = 'HRS' ELSE	$
		IF strpos(ARR,'MIN') NE -1 THEN NEWTUNIT = 'MIN' ELSE	$
		IF strpos(ARR,'SEC') NE -1 THEN NEWTUNIT = 'SEC' ELSE	$
										NEWTUNIT = 'HRS'
	ENDELSE

	; The time origin is set only once. If a new time origin is set in a second file than
	; the difference with the time origin already in memory is stored in DJD.

	IF n_elements(YEAR0) EQ 0 THEN BEGIN	; If not set yet, set time origin and unit of time
		YEAR0 = NEWYEAR
		DOY0  = NEWDOY
		TUNIT = NEWTUNIT
		message, /info, 'Time origin :'+string(YEAR0)+' DOY'+string(DOY0)+' , units : '+TUNIT
	ENDIF
	CvTime, /indoy, yr=YEAR0	, doy=DOY0			, jd=JD0	; Old time origin in Julian days
	CvTime, /indoy, yr=NEWYEAR	, doy=NEWDOY+OFFDOY	, jd=NEWJD	; New time origin in Julian days
	DJD = NEWJD-JD0							; Offset of new series start time with time origin

	CASE NEWTUNIT OF
	'DAY': DJD = DJD
	'HRS': DJD = 24*DJD
	'MIN': DJD = 1440*DJD
	'SEC': DJD = 86400*DJD
	ENDCASE

	A = ''
	readf, LUN, A				; Read record with row and column information
	A = fix( flt_string(A) )
	NROW = A[0]
	NCOL = A[1]
	IF NHX EQ 0 THEN BEGIN		; Not a HXIS file
		B1 = replicate(-1,NCOL/2)
		B2 = B1 				; NCOL/2 is # time series
		B3 = B1
	ENDIF

	ARR = fltarr(NCOL,NROW)		; Create array to contain time series
	A = fltarr(NCOL)
	FOR I=0,NROW-1 DO BEGIN		; Read next block of data (NROW records)
		readf, LUN, A
		ARR[*,I] = A
	ENDFOR

	IPOLD = IP
	FOR I=0,NCOL-1,2 DO BEGIN	; Pick up the times from odd columns
		I2 = I/2				; # time series
		T  = ARR[I,*]
		A  = where(T NE -1)		; Column of times
		IF A[0] EQ -1 THEN BEGIN; Skip empty time series
			echo, 'EMPTY TIME SERIES IGNORED'
			goto, SKIP
		ENDIF

		T = T[A]
		F = (ARR[I+1,*])[A]		; Column of function values
		T = DJD+T				; Add offset to current time origin

		IF IP EQ -1 THEN BEGIN	; No series loaded yet
			IP   = 0
			SIP  = strcompress(IP,/remove_all)
			BTYP = B1[I2]
			BLOW = B2[I2]
			BUP  = B3[I2]
			STAT = 'T'+SIP+' = T  &  F'+SIP+' = F'

		ENDIF ELSE BEGIN
			B = -1
			IF NHX NE 0 THEN BEGIN		; For HXIS data check data type
				B = where(BTYP EQ B1[I2] AND BLOW EQ B2[I2] AND BUP EQ B3[I2])
			ENDIF

			IF B[0] EQ -1 THEN BEGIN	; No HXIS data or new data type
				IP   = IP+1
				SIP  = strcompress(IP,/remove_all)
				BTYP = [BTYP,B1[I2]]
				BLOW = [BLOW,B2[I2]]
				BUP  = [BUP ,B3[I2]]
				STAT = 'T'+SIP+' = T & F'+SIP+' = F'

			ENDIF ELSE BEGIN
				A = 0B
				echo, 'Merge time series', A, /logical
				IF A THEN BEGIN
					SIP  = strcompress(B[0],/remove_all)
					TS   = 'T'+SIP
					FS   = 'F'+SIP
					STAT = TS+' = ['+TS+',T] & '+FS+' = ['+FS+',F] & A = sort('+TS+') & '+TS+' = '+TS+'(A) & '+FS+' = '+FS+'(A)'
				ENDIF ELSE BEGIN
					IP   = IP+1
					SIP  = strcompress(IP,/remove_all)
					BTYP = [BTYP,B1[I2]]
					BLOW = [BLOW,B2[I2]]
					BUP  = [BUP ,B3[I2]]
					STAT = 'T'+SIP+' = T & F'+SIP+' = F'
				ENDELSE

			ENDELSE

		ENDELSE
		IF NOT execute(STAT) THEN message, 'Problem reading '+FILE

		IF IP EQ 5 THEN BEGIN
			free_lun, LUN
			echo, 'ALL DATA ARRAYS ARE NOW FULL
			A = NCOL/2+IPOLD-IP
			IF A NE 0 THEN echo, '('+strcompress(A)+' TIME SERIES HAVE NOT BEEN READ)'
			RETURN
		ENDIF

SKIP:
	ENDFOR

ENDWHILE

free_lun, LUN

RETURN  &  END