;+
; NAME:
;	vu_read
; PURPOSE:
;	Read output file from tomography
; CATEGORY:
;	I/O
; CALLING SEQUENCE:
	FUNCTION vu_read, hdr,	$
		nobad	= nobad 	, $
		get_roi = get_roi	, $
		get_full= get_full	, $
		fss 	= fss		, $
		shfts	= shfts 	, $
		losxs	= losxs 	, $
		silent	= silent	, $
		count	= count 	, $
		roi_offset=roi_offset
; INPUTS:
;	hdr			array; type: structure
;					information extracted from file headers
; OPTIONAL INPUTS:
;	/nobad		by default 'bad' values are replaced by !values.f_nan.
;					The 'bad value' indicator is read from the file header, or,
;					if not found in the header is assumed to be -9999.9990
;					Set /nobad to suppress this 'bad value' replacement
;	/get_full	extract full array (not just region of interest)
;	/get_roi	extract region of interest only (this is the default)
;
;	roi_offset=roi_offset
;				scalar; type: float: default: 0.0
;					offset to be added to the region of interest specified in the
;					header (not used if /get_full is set)
; OUTPUTS:
;	Result		array; type: float
;					3D matrix (invalid values marked by large negative value)
;					Either velocity/density or radial/azimuthal magnetic field
; OPTIONAL OUTPUT PARAMETERS:
;	count		scalar; type: integer
;					# files read; this will be either 0 of the number of valid
;					nv3d* or wson* files. If there is a read error on the first of
;					these files then count=0 is returned and none of the output
;					arrays will exist. If it happens for the second or later file
;					then its slot in the output arrays is filled with bad values.
;	silent=silent
;				if set, informational messages are suppressed
;
;	The following arrays may not be present in the file:
;
;	fss=fss		array[*,*,2]; type: float
;					source surface maps
;	shfts=shfts
;				array[*,*,*]; type: float
;					3D shift matrix (shft[*,*,0]=0)
;	losxs=losxs
;				array[*,*,2]; type: float
;					# line of sight crossings at source surface
; INCLUDE:
	@compile_opt.pro				; On error, return to caller
; CALLS:
;	flt_read, BadValue, IsType, SuperArray, destroyvar, TimeGet, CarringtonT
;	vu_get, vu_set, InitVar
; PROCEDURE:
; >	The files are read using procedure flt_read. If there is a read error
;	then the value zero is returned.
; >	The shift array is assumed to be stored at the end of the file
; MODIFICATION HISTORY:
;	AUG-1999, Paul Hick (UCSD/CASS)
;	MAR-2001, Paul Hick (UCSD/CASS)
;		rewrite to deal with new headers of tomography files
;	NOV-2001, Paul Hick (UCSD/CASS)
;		added keyword roi_offset to adjust the region of interest
;		specified in the header.
;	AUG-2002, Paul Hick (UCSD/CASS)
;		modification to allow for reading of files with different
;		number of radial distances. The lowest number of distances
;		is read from each file.
;	JUN-2003, Paul Hick (UCSD/CASS)
;		Fixed bug which prevented reading of source surface map fss correctly.
;	NOV-2003, Paul Hick (UCSD/CASS)
;		Fixed bug in determining indices for ROI. If the ROI edges fell exactly
;		halfway two gridpoints, one edge would sometimes be rounded up, the
;		other rounded down.
;	NOV-2003, Paul Hick (UCSD/CASS)
;		Modified handling of read errors. If the error occurs on the first file
;		then nothing is returned; if it happens for a later file its slot
;		in the output arrays is filled with bad values.
;	OCT-2006, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;		For magnetic data the third component is added as zeroes if it
;		is absent. The radial (1st) and tangential (2nd) component are always
;		present. The normal (3rd) component is not (yet) present in the
;		magnetic files.
;-

count = n_elements(hdr)

InitVar, nobad	 , /key
InitVar, get_roi , /key
InitVar, get_full, /key
InitVar, silent	 , 0

get_roi = get_roi OR (1-get_full)

InitVar, roi_offset, 0.0

; Pick up the dimensions. nLng and nLat must be the same for all headers
; (this was checked by vu_select). The lowest nRad value is used.

nLng  = vu_get(hdr[0], /nlng)
nLat  = vu_get(hdr[0], /nlat)

nRad  = vu_get(hdr, /nrad)
nRad0 = min(nRad)
OneRad = ( where(nRad NE nRad0) )[0] EQ -1

; Read all the files. At this point we know that all files are valid
; nv3d* or wson* files, so there should be no read error.

nx = 0
ny = 0

FOR ifile=0,count-1 DO BEGIN

	hdri = hdr[ifile]

	sts = 0
	status = 0
	;CASE flt_read( vu_get(hdri,/file), data, comments=comments, nx=nx, ny=ny, silent=silent+1) OF
	;IF vu_get(hdri, /prefix) EQ 'nv3h' THEN BEGIN
	;	command = "sts = nv3h_read_gz( vu_get(hdri,/file), data, comments=comments, nx=nx, ny=ny, silent=silent+1)"
	;	status = execute(command,1,0)
	;ENDIF

	; status = 0 if nv3h_read failed or if the prefix is not nv3h
	IF status EQ 0 THEN sts = flt_read( vu_get(hdri,/file), data, comments=comments, nx=nx, ny=ny, silent=silent+1)

	CASE sts OF
	0: BEGIN

		; Read errors should be darn near impossible at this point (the file exists and
		; is a valid nv3d* or wson* file). If an error occurs anyway on the first file,
		; then we return zero (no data read); if an error occurs for ifile > 0 (after
		; succesfully reading the first file then we fill the slot 'ifile' with bad values.

		tmp = 'unexpected read error on '+hide_env(vu_get(hdri,/file))

		CASE ifile OF
		0: BEGIN		
			say, threshold=2, silent=silent, [tmp+', no data read','... ABORTING']
			destroyvar, ff, fss, shfts, losxs
			count = 0
			RETURN, -1						; Exit on read error
		END
		ELSE: BEGIN
			say, threshold=2, silent=silent, tmp+', assuming bad data'

			ff[*,*,*,*,ifile] = BadValue(ff)

			IF ss_present   THEN fss  [*,*,*,ifile] = BadValue(fss  )
			IF losx_present THEN losxs[*,*,*,ifile] = BadValue(losxs)
			IF shft_present THEN shfts[*,*,*,ifile] = BadValue(shfts)
		END
		ENDCASE

	END
	1: BEGIN

		;say, threshold=0, silent=silent, TimeGet(vu_get(hdri,/uttime), /ymd )

		; Check for 3D density and velocity matrix (must be present),
		; or magnetic field data. For velocity/normalized density there may be
		; two hits: one for the 3D data and one for the source surface.

		CASE vu_get(hdri, /plasma) OF
		0: target = '; Radial (r^2*B_r)/tangential (r*B_t) magnetic field'
		1: target = '; Velocity (km s^-1)/normalized density (cm^-3)'
		ENDCASE

		i = where( strmid(comments,0,strlen(target)) EQ target )
		ss_present = n_elements(i) EQ 2

		; Check for los crossings maps (may not be present)

		target = '; Velocity/normalized density line of sight crossings at source surface'

		i = (where( strmid(comments,0,strlen(target)) EQ target ))[0]
		losx_present = i NE -1

		; The shift array is stored as a single block of nLat*nRad records at the end of the file
		; (may not be present)

		target = '; Shift in Carrington variable (add to go back to source surface)'
		i = (where( strmid(comments,0,strlen(target)) EQ target ))[0]
		shft_present = i NE -1

		; The 3D data are at the start of the matrix

		nRad = vu_get(hdri, /nrad)					; Pick up nRad

		nrows_present = (size(data,/dim))[1]
		nrows_expect2 = nLat*(2*nRad+2*ss_present+2*losx_present+(2-1)*shft_present*nRad)
		nrows_expect3 = nLat*(3*nRad+3*ss_present+3*losx_present+(3-1)*shft_present*nRad)

		; Expected data:
		;	Always		: nLat*nRad*2 rows for main data array
		;	ss_present	: nLat*     2 rows for source surface
		;	losx_present: nLat*     2 rows for los crossings
		;	shft_present: nLat*nRad   rows for shift matrix

		IF nrows_present  LT nrows_expect2 THEN BEGIN
			print, 'Found   :', nrows_present,' rows'
			print, 'Expected:', nrows_expect2,' rows'
			tmp = 'corrupt file '+hide_env(vu_get(hdri,/file))
			message, tmp+', no data read'
		ENDIF

		nlayer = nrows_present EQ nrows_expect2 ? 2 : nrows_present EQ nrows_expect3 ? 3: 0

		m = 0
		n = nLat*nRad*2-1

		f = reform( data[*,m:n], nLng,nLat,nRad,2)

		; Two source surface maps may follow the 3D data

		IF ss_present THEN BEGIN
			m = n+1
			n = m+2*nLat-1
			fs = reform( data[*,m:n], nLng,nLat,2)
		ENDIF

		; Line of sight crossings and shift matrix may be present in nv3d* files.

		IF losx_present THEN BEGIN
			m = n+1
			n = m+2*nLat-1
			losx = reform( data[*,m:n], nLng,nLat,2)
		ENDIF

		IF shft_present THEN BEGIN
			m = n+1
			n = m+nLat*nRad-1
			shft = reform( data[*,m:n], nLng,nLat,nRad)
		ENDIF

		; If the file contains more than nRad0 distances, discard the upper
		; portions. Make sure to update the header !

		IF nRad NE nRad0 THEN BEGIN
			f = f[*,*,0:nRad0-1,*]						; Discard everything above nRad0
			IF shft_present THEN shft = shft[*,*,0:nRad0-1]

			hdr[ifile] = vu_set(hdri, nrad = nRad0)		; Update header
			say, threshold=1, silent=silent, vu_get(hdri,/file)+ $
				', retain '+strcompress(nRad0, /rem)+'/'+strcompress(nRad, /rem)+' distances'
		ENDIF

		xcrange = vu_get(hdri, /array_range)
		xcroi   = vu_get(hdri, /roi_range)+get_roi*roi_offset	; Offset applied only if get_roi is set

		CASE get_roi OF

		0: say, threshold=0, silent=silent, vu_get(hdri,/prefix)+' '+$
				TimeGet(vu_get(hdri,/uttime), /ymd )+'; full:'+strjoin(strcompress(xcrange))

		1: BEGIN

			i = (xcroi-xcrange[1])/(xcrange[0]-xcrange[1])*(nLng-1)
			i = round(i[0])+[0, round(i[1]-i[0])]		; Round both edges in same direction
			xcroi = xcrange[1]+float(i)/(nLng-1)*(xcrange[0]-xcrange[1])

			j = i[0]
			i = i[1]

			say, threshold=0, silent=silent, vu_get(hdri,/prefix)+' '+$
				TimeGet(vu_get(hdri,/uttime), /ymd )+'; roi:'+strjoin( strcompress(xcroi) )+	$
				' ['+strcompress(i,/rem)+','+strcompress(j,/rem)+']/[0'+','+strcompress(nLng-1,/rem)+']'+ $
				' ('+strcompress(vu_get(hdri,/marker),/rem)+')'

			;IF i LT 0 OR j GE (size(f,/dim))[0] THEN BEGIN
			;	say, threshold=0, silent=silent, $
			;		[ 'longitude index range is [0,'+strcompress((size(f,/dim))[0],/rem)+']', $
			;		  'calculated ROI range is ['+strcompress(i,/rem)+','+strcompress(j,/rem)+']' ]
			;	message, vu_get(hdri,/file)+', error extracting region of interest'
			;ENDIF

			; Extract the center of the three maps

			f = f[i:j,*,*,*]

			IF ss_present   THEN fs   = fs  [i:j,*,*]
			IF losx_present THEN losx = losx[i:j,*,*]	; 3rd dim is 2
			IF shft_present THEN shft = shft[i:j,*,*]	; 3rd dim is nRad

			hdr[ifile].array = xcroi
			hdr[ifile].roi   = xcroi
			hdr[ifile] = vu_set(hdr[ifile],nLng = j-i+1)

		END

		ENDCASE

		; Check for bad values

		Bad = vu_get(hdri, /bad)

		IF NOT NoBad AND finite(Bad) NE 0 THEN BEGIN
			IF IsType(f, /float) THEN Bad = float(Bad)
			i = where(f EQ Bad)
			IF i[0] NE -1 THEN f[i] = BadValue(f)
		ENDIF


		; After reading the first file (ifile = 0) create the output arrays ff.

		CASE ifile OF

		0   : BEGIN
			base_size = size(f)
			ff = SuperArray(f, count, /trail)
			IF ss_present   THEN fss   = SuperArray(fs  , count, /trail)
			IF losx_present THEN losxs = SuperArray(losx, count, /trail)
			IF shft_present THEN shfts = SuperArray(shft, count, /trail)
		END

		ELSE: BEGIN

			IF (where(base_size ne size(f)))[0] NE -1 THEN BEGIN
				say, threshold=0, silent=silent, [vu_get(hdri,/file), 'probably an error processing header']
				whatis, ff,f
			ENDIF

			ff[*,*,*,*,ifile] = f

			IF ss_present   THEN fss  [*,*,*,ifile] = fs
			IF losx_present THEN losxs[*,*,*,ifile] = losx
			IF shft_present THEN shfts[*,*,*,ifile] = shft
		END

		ENDCASE

		; If all files have the same number of levels, and only the data
		; array is present, then setting nx and ny will speed up reading
		; of files with flt_read (skips the preread).

		IF OneRad AND NOT ss_present AND NOT losx_present AND NOT shft_present THEN BEGIN
			nx = nLng
			ny = nLat*nRad*2
		ENDIF ELSE BEGIN
			nx = 0
			ny = 0
		ENDELSE

	END
	ENDCASE

ENDFOR

ff = reform(ff, /overwrite)

say, threshold=1, silent=silent, 'got '+strcompress(count, /rem)+' files'

IF hdr[0].magnetic THEN BEGIN

	; For magnetic data check for the presence of the 3rd (normal) component
	; If not present add a zero component.

	IF size(ff,/n_dim) GE 4 THEN			$
		IF (size(ff,/dim))[3] EQ 2 THEN 	$
			boost, ff, pushdim=3, plus=1
ENDIF

RETURN, ff  &  END