; Copyright (c) 1998-2011 A.P. Hitchcock  All rights reserved
;+
;NAME:
;  READ_XRF_IMG
;
;LAST CHANGED: -----------------------------------  13-Mar-11 (aph)
;
;PURPOSE:
;	This function reads in XRF images from Twinmic STXM; Inca (Oxford SDD) on Tescan SEM
;
;CATEGORY:
;	STAND ALONE: read in utility
;
;CALLING SEQUENCE:
;	Result = READ_XRF_IMG([file=file, NOF=nof, group=group,TYPE=TYPE, $
;                  VERBOSE=VERBOSE, SILENT=SILENT, PARAM=PARAM, _extra=e])
;
;CALLED FROM AXIS:
;	Read~Images~XRF
;
;INPUTS: none. All input parameters are passed as keywords.
;
;KEYWORDS:
;	FILE 	filename to read
;	NOF 	no filter (default filter = '*.axb')
;	GROUP	group leader (Axis_ID if called from aXis2000 and axis)
;	SILENT  if set, do not display feedback on progress
;	PARAM	if set, ask user for parameters of Twinmic image
;	TYPE	type of image (sum, X_mom, Y_mom) selected automatically
;	SCOPE	instrument generating data (twinmic, inca)
;	VERBOSE print log information about processing
;	EXTRA 	passed on parameters
;
;OUTPUTS:
;	RESULT = AXIS 1d or 2d structure
;
;COMMON BLOCKS:
;	@AXIS_COM	standard set of common blocks
;	@BSIF_COM   for NSLSL *.nc
;
;MODIFICATION HISTORY:
; (14-Dec-09 aph) adapted from AXB_LOAD
; (14-Feb-10 aph) extend to SEM Oxford INCA XRF images (*.tsv)
; (27-Mar-10 aph) add flip to inca images to make consistent with Tescan SEM images
;        add trim and multi-file *.tsv read
; (25-May-10 aph) inca images read in - adapt to FIELD0001 versus FIELD001 names
; (10-mar-11 aph) auto find sub-images from INCA (*.tsv)
; (13-mar-11 aph) get multi-file read in going
;-

Function READ_XRF_IMG, file=file, NOF=nof, group=group, silent=silent, $
       type=type, scope=scope, PARAM=PARAM, verbose=verbose, _extra=e
@axis_com
@bsif_com
on_error,2

IF NOT keyword_set(file) then begin
	fltr='*.img'
	if keyword_set(scope) then begin
		if scope EQ 'inca' then fltr='*.tsv'
	endif
	if keyword_set(nof) then fltr='*'
	file = pickfile2(/READ, FILTER=fltr, /LPATH, DEFPATH=defpath)
ENDIF
s = ''
if file EQ '' THEN RETURN, s  ; bail-out if no filename
check = findfile(file)		 ; bail-out if non-existent file
if strlen(check(0)) EQ 0 then begin
	axis_log, 'Cannot find ' + file
	return, s
endif
t = ax_name(file)					; file format: yymmdd_hhmmss_(type).ssv)
fileshort = t(1)
if NOT keyword_Set(scope) then begin
	if t(2) EQ 'img' then scope ='twinmic' else scope = 'inca'
endif

CASE scope OF
 'twinmic' : BEGIN
	openr, lun, file, /get_lun		; open as text
	ssvt = ''
	readf,lun, ssvt
	; print, ssvt
	ssvz=parsestring(ssvt)		; convert to array of numbers
	t = size(ssvz)
	if t(0) EQ 1 then begin
		nx = t(1)						; assumes image file varies x-dimension most rapidly
	endif else begin					; use ROTATE afterward to obtain correct sense of image
		axis_log, 'error on read of Twinmic file ' + fileshort
		return, s
	endelse
	ny = 1
	dt = ssvz
	while not  eof(lun) do begin
		readf,lun, ssvz
	;	print, ssvz
		dt = [dt,ssvz]
		ny = ny + 1
	endwhile
	d = rotate(float(reform(dt,nx,ny)),0)
	close, lun & free_lun, lun

; extract parameters from an associated parameter file
	fltr='*.prm'
	filep = pickfile2(/READ, FILTER=fltr, /LPATH, DEFPATH=defpath)
	if filep EQ '' THEN goto, skip_param
	check = findfile(filep)		 ; bail-out if non-existent file
	if strlen(check(0)) EQ 0 then goto, skip_param
	openr, lun, filep, /get_lun		; open as text
	; ------- Read header file into a single string (max = 32600 characters)  --------
	line = ' '
	all = ''
	count = 0
	while not  eof(lun) do begin
		readf,lun, line
	;	if keyword_set(verbose) then print, line
		all = all + line
		count = count + 1			; the number of lines is used in spectral readin
	endwhile
	close, lun & free_lun, lun
	on_ioerror,null

	x_mid = float(ax_strmid(all, 'X center of image (micron) =', 'Y'))
	y_mid = float(ax_strmid(all, 'Y center of image (micron) =', 'X'))
	x_range = float(ax_strmid(all, 'X scan range (micron) =', 'Y'))
	y_range = float(ax_strmid(all, 'Y scan range (micron) =', 'No'))
	x_step = x_range/nx
	y_step = y_range/ny
	x = findgen(nx)
	y = findgen(ny)
	x= x_mid - 0.5*x_range  + x*x_step
	y= y_mid - 0.5*y_range  + y*y_step
	dwell = float(ax_strmid(all, 'Dwell time (ms) =', 'Total'))
	if keyword_set(group) then begin
		energy = get_num(prompt='Photon energy (eV)', val=energy, group=group)
	endif else begin
		energy = get_num(prompt='Photon energy (eV)', val=energy)
	endelse

	skip_param:				; no parameter file - ask user if param is set
	if keyword_set(param) then begin
		if keyword_set(group) then begin
			energy = get_num(prompt='Photon energy (eV)', val=energy, group=group)
	;		dwell = get_num(prompt='Dwell (ms)', val=dwell, group=group)
			dl = get_text(prompt='Header', val=dl, group=group)
		endif else begin
			energy = get_num(prompt='Photon energy (eV)', val=energy)
	;		dwell = get_num(prompt='Dwell (ms)', val=dwell)
			dl = get_text(prompt='Header', val=dl)
		endelse
	 endif
  END

; -------- read in INCA XRF map --------------------
'inca' : BEGIN
	tmp = read_ascii(file)
;	print, tag_names(tmp)
;	return, tmp
	d = tmp.(0)
	d = reverse(d,2)
	t = size(d)
	nx = t(1)  & ny = t(2)
	x = findgen(nx)
	y = findgen(ny)
; ---- set default values in case user has not set up row, col, limits
	nx_start = 0
	nx_stop = nx-1
	nxt = nx_stop - nx_start+1
	ny_start = 0
	ny_stop = ny-1
	nyt = ny_stop - ny_start+1
	map_width = nxt
; ---- determine row and column limits based on excluding completely non-zero rows, columns
	if not keyword_set(silent)then begin		;  ASSUME silent means multiple and thus all subsequent are same
 ; ------------ determine the indices of non-zero rows
		test = 0  & i = 0
		while test EQ 0 do begin
			for j = 0, nx-1 do test = test + d(j,i)
	;		print, i, test
			i = i + 1
		endwhile
		ny_start = i
		if not keyword_set(silent) then axis_log, 'first non-zero row is ' + strtrim(string(ny_start),2)
		i = ny_start
		while test GT 0 AND i LT ny-1 do begin
			i = i + 1
			test = 0
			for j = 0, nx-1 do test = test + d(j,i)
;			print, i, test
		endwhile
		ny_stop = i-1
		if not keyword_set(silent) then axis_log, 'last non-zero row is ' + strtrim(string(ny_stop),2)

; ------------ determine the indices of non-zero columns
		nx_start = 0
		nx_stop = nx-1
		line_sum = fltarr(nx)
		print, 'Find index of first non-zero column'
		for i=0, ny-1 do line_sum = line_sum + d(i,*)
		test = 0  & j = 0
		while test EQ 0 do begin
			test = test + line_sum(j)
;			print, j, test
			j = j + 1
		endwhile
		nx_start = j-1
;		print, 'index of first non-zero', nx_start, ' preceding and that point ', line_sum(nx_start-1), line_sum(nx_start)
		if not keyword_set(silent) then axis_log, 'first non-zero column is ' + strtrim(string(nx_start),2)

		index = where(line_sum GT 0, count)
		if count NE 0 then nx_stop = max(index)
;		print, 'index of last non-zero', nx_stop, ' that point and next ', line_sum(nx_stop), line_sum(nx_stop+1)
		if not keyword_set(silent) then axis_log, 'last non-zero column is ' + strtrim(string(nx_stop),2)

		nxt = nx_stop - nx_start+1
		if n_elements(xrf_energy)  EQ 0 then  energy = 20000.
		if n_elements(xrf_dwell)   EQ 0 then  dwell  = 100.

		if keyword_set(group) then begin
;			 nx_start = get_num(prompt='start column', val=nx_start, group=group)
;			 nx_stop  = get_num(prompt='stop column', val=nx_stop, group=group)
;			 ny_start = get_num(prompt='start row', val=ny_start, group=group)
;			 ny_stop  = get_num(prompt='stop row', val=ny_stop, group=group)
			 map_width = get_num(prompt='map width (um)', val=map_width, group=group)
			 energy =  get_num(prompt='beam energy (eV)', val=energy, group=group)
			 dwell =  get_num(prompt='dwell (ms)', val=dwell, group=group)
		endif else begin
;			 nx_start = get_num(prompt='start column', val=nx_start)
;			 nx_stop  = get_num(prompt='stop column', val=nx_stop)
;			 ny_start = get_num(prompt='start row', val=ny_start)
;			 ny_stop  = get_num(prompt='stop row', val=ny_stop)
			 map_width = get_num(prompt='map width (um)', val=map_width)
			 energy =  get_num(prompt='beam energy (eV)', val=energy)
			 dwell =  get_num(prompt='dwell (ms)', val=dwell)
		endelse
    endif else begin		; OK silent is set so we assume same as first file read
	; ---- if row limit is defined and keyword is silent use prior saved value
			if n_elements(xrf_cols_start) GT 0 then nx_start = xrf_cols_start
			if n_elements(xrf_cols_stop)  GT 0 then nx_stop  = xrf_cols_stop
			if n_elements(xrf_rows_start) GT 0 then ny_start = xrf_rows_start
			if n_elements(xrf_rows_stop)  GT 0 then ny_stop  = xrf_rows_stop
			if n_elements(xrf_map_width)  GT 0 then map_width  = xrf_map_width
			if n_elements(xrf_energy)  GT 0 then energy  = xrf_energy
			if n_elements(xrf_dwell)  GT 0 then dwell  = xrf_dwell
	endelse

	if ny_stop GT ny-1 then ny_stop = ny-1
	nyt = ny_stop-ny_start+1
	if not keyword_set(silent) then axis_log, '# of non-zero rows is ' $
	    + strtrim(string(fix(nyt)),2) + ' of ' + strtrim(string(ny),2)
	if nx_stop GT nx-1 then nx_stop = nx-1
	nxt = nx_stop-nx_start+1
	if not keyword_set(silent) then axis_log, '# of non-zero columns is ' $
	    + strtrim(string(fix(nxt)),2) + ' of ' + strtrim(string(nx),2)

; ------ trim to non-zero rows
	if nyt LT ny then begin		; remove rows below nyt
		dt = d(*,ny_start:ny_stop)
		d = dt
		ny = nyt
		y = findgen(ny)
;		print, 'y', d(0,*)
	endif
; ------ trim to non-zero columns
	if nxt LT nx then begin		; remove columns above nxt
; ===== KLUGE ===== KLUGE ===== KLUGE ===== KLUGE ===== KLUGE ===== KLUGE ===== KLUGE ===== KLUGE
		kluge = 5	; for some reason the read-in starts 5 lines BEFORE IT SHOULD
; ===== KLUGE ===== KLUGE ===== KLUGE ===== KLUGE ===== KLUGE ===== KLUGE ===== KLUGE ===== KLUGE
		dt = d(nx_start+kluge+1:nx_stop,*)
		d = dt
		nx = nxt-1-kluge
		x = findgen(nx)
	endif

;----- define pixel size
	pix_siz = map_width/nx		; NB  this is stored in axis_com
	x = findgen(nx)*pix_siz
	y = findgen(ny)*pix_siz

;	update stored value of column and row limits and width of map
	xrf_cols_start = nx_start
	xrf_cols_stop  = nx_stop
	xrf_rows_start = ny_start
	xrf_rows_stop  = ny_stop
	xrf_map_width = map_width
	xrf_energy = energy
	xrf_dwell = dwell

  END

  ELSE: axis_log, 'READ_XRF_IMG: no code for scope type ' + scope

ENDCASE

;  ------------ convert to axis 2d structure
dl = fileshort		; default header

if dwell GT 1000 then begin
	xl = string(FORMAT='("x (um)     E = ",f6.0," eV     dwell = ",f7.0," sec")', $
      energy, dwell/1000.)
endif else begin
    xl = string(FORMAT='("x (um)     E = ",f6.0," eV     dwell = ",f7.0," ms")', $
      energy, dwell)
endelse
yl = scope + '  XRF map'

text = string(format= $
  '(i4," x",i4," XRF image read from ",A,".   E =",f8.2," eV. Dwell =",g13.6," ms.")', $
   nx, ny, fileshort, energy, dwell)
if not keyword_set(silent) then axis_log, text

s = {t:'2d', x:x, y:y, d:d, e: energy, xl:xl, yl:yl, dl: dl}	; add E

return, s

END
