; Copyright (c) 1998-2015 A.P. Hitchcock  All rights reserved
;+
;NAME:
;		READ_TXM_ALBA
;
;LAST CHANGED: ----------------------------------- 	29 Mar 2015
;
; PURPOSE:
;	This function reads an image or stack from an hdf format file
;
; CATEGORY:
;	Data read-in. Stand alone or used in aXis2000
;
; CALLING SEQUENCE:
;	Result = READ_TXM_ALBA(file=file, one=one)
;
; KEYWORDS:
;	file	filename
;	norm	ring current to which to normalize   (100 mA is initial ring current)
;	one		if set, read only 1 image
;   silent	if set, do not print tracing comments
;	three	if set, use 3rd format established by Mistral (Mark Rosanes)
;
; OUTPUTS:
;	axis format image if one image; *.(ncb,dat) if stack
;
; LIMITATIONS
;	for large files, this code will not work on 32-bit architecture systems (XP)
;
; COMMON BLOCKS:
;	axis
;	stack_process_com
;	volume_data
;	bsif_com
;
; MODIFICATION HISTORY:
; (29-Mar-15 aph) first verion
;-


FUNCTION read_txm_alba, file=file, one=one, silent=silent, norm=norm, three = three
@axis_com
@stack_process_com
COMMON volume_data, image_stack
@bsif_com
on_error,2

; ----- determine if aXis2000 is running to set a local switch
if aXis_ID GT 0 then begin
	print, 'aXis_ID ', aXis_ID
	axis_on = 1
endif else axis_on=0

; -------- establish ring current to normalize to		(Alba starts at 100 mA then decays)
if not keyword_set(norm) then norm = 100.
norm = float(norm)

IF NOT keyword_set(file) then begin
	fltr = '*.hdf5'
   	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_path = t(0)
fileshort = t(1)

; ------ determine type of structure - there are 3 as of March 2015
a = H5_browser(file, /dialog_read)
test = tag_names(a)
type = 'null'
for i = 0, n_elements(test)-1 do begin
;	print, test(i)
	if STRUPCASE(test(i)) EQ 'NXTOMO' then type = 'one'			; first one sent
	if STRUPCASE(test(i)) EQ 'FASTALIGNED' then type = 'three'	; Mark Rosanes says this will be default
endfor
if type EQ 'null' then begin
	axis_log,'ALBA.Mistral: unknown hdf5 data structure'
	return, a
endif

if keyword_set(one) then begin	; single image

CASE type OF

'three': BEGIN
; XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
; -------- read in one image from type 3 hdf5 file
	nx = a.fastaligned.spec_aligned.pixel_columns._data
	ny = a.fastaligned.spec_aligned.pixel_rows._data

; --------- get number of images and ask user which one they want if more than 1
	npts = n_elements(a.fastaligned.ExpTimes._data)
	n_img=0
  	if npts GT 1 then begin
		text = strtrim(string(fix(npts)),2) + ' images. Select the image to read'
  		axis_log, fileshort +' contains ' + text
;  		print, 'axis_on = ', axis_on
  		if axis_on EQ 1 then $
  			 n_img = get_num(Prompt = text, val = n_img, group = axis_ID) $
		else n_img = get_num(Prompt = text, val = n_img)
		if n_img LT 0 or n_img GT npts then begin
			axis_log, 'non-existent image selected. Aborting.'
			return, s
		endif
	endif
; ------- get pixel size
	xstep = 0.06  &  ystep = 0.06	; default values
	test = tag_names(a.fastaligned.spec_aligned)
	for i = 0, n_elements(test)-1 do begin
		if test(i) EQ 'X_PIXEL_SIZE' then $
		     xstep = a.fastaligned.spec_aligned.x_pixel_size._data(n_img)
		if test(i) EQ 'Y_PIXEL_SIZE' then $
			ystep = a.fastaligned.spec_aligned.y_pixel_size._data(n_img)
	endfor
;	help, xstep, ystep

; ------ get dwell time
	dwell_time =  a.fastaligned.ExpTimes._data(n_img)

; ------- get beam current
	ring_current=a.fastaligned.currents._data(n_img)

; ----- get incident photon energy (eV) -----
	ev = a.fastaligned.energy._data(n_img)

; ---- get data
  	d = a.fastaligned.spec_aligned._data(*,*,n_img)*norm/ring_current
  	END

; XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
; -------- read in one image from type 1 hdf5 file
	'one' :	BEGIN ; first format sent; derived from their tomo format

  	npts = a.nxtomo.instrument.sample.data.number_of_frames._data
	nx = a.nxtomo.instrument.sample.data.image_width._data
	ny = a.nxtomo.instrument.sample.data.image_height._data

; ---- select image to read
	n_img=0
  	if npts GT 1 then begin
		text = strtrim(string(fix(npts)),2) + ' images. Select the image to read'
  		axis_log, fileshort +' contains ' + text
;  		print, 'axis_on = ', axis_on
  		if axis_on EQ 1 then $
  			 n_img = get_num(Prompt = text, val = n_img, group = axis_ID) $
		else n_img = get_num(Prompt = text, val = n_img)
		if n_img LT 0 or n_img GT npts then begin
			axis_log, 'non-existent image selected. Aborting.'
			return, s
		endif
	endif

	; ------- get pixel size
		xstep = a.NXTOMO.instrument.sample.x_pixel_size._data(0)
		ystep = a.NXTOMO.instrument.sample.y_pixel_size._data(0)

; ------ get dwell time
		dwell_time =  A.NXTOMO.instrument.sample.ExpTimes._data(n_img)

; -------- ring current
	ring_current =  A.NXTOMO.instrument.sample.current._data(n_img)

; ----- get incident photon energy (eV) -----
  		ev = A.NXTOMO.instrument.source.energy._data(n_img)
; ---- get data
	  	d = (a.nxtomo.instrument.sample.data._data(*,*,n_img) + 32767.)*norm/ring_current
	END
	ENDCASE

	x_start = 0.  & x_stop = nx*xstep
	y_start = 0.  & y_stop = ny*ystep
	x = findgen(nx) & x = x_start +x*xstep
	y = findgen(ny) & y = y_start +y*ystep
	help, x, y

	t = size(d)
	xl = 'Position (um)     ' + string(ev, format='(F7.2)') + ' eV   ' $
	      + string(dwell_time, format='(F6.2)')+ ' sec'
	yl = 'Alba-Mistral TXM HDF5'
	if npts GT 1 then dl = fileshort + '  Image ' + strtrim(string(fix(n_img)),2) $
	            else  dl = fileshort
	s = {t:'2d', x:x, y:y, d:d, e:ev, xl:xl, yl:yl, dl:dl}
;	help, s, /struct

; -------- print tracing for trouble shooting
	IF NOT keyword_set(silent) then  begin
		print, 'Read type ', type,' Alba - Mistral image from ', file
		print, 'nx, ny, #images ', nx, ny, npts
		print, 'Ring current (mA) = ', ring_current
		print, 'Dwell time (s) = ', dwell_time
		print, 'Energy (eV) = ', ev
	ENDIF
RETURN, s

endif else begin
; XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
;  ------------  READ in full  HDF5 format Alba TXM stack
; XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX

CASE type OF

'three': BEGIN
; XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
; -------- read in stack from type 3 hdf5 file
	nx = a.fastaligned.spec_aligned.pixel_columns._data
	ny = a.fastaligned.spec_aligned.pixel_rows._data

; --------- get number of images and ask user which one they want if more than 1
	npts = n_elements(a.fastaligned.ExpTimes._data)

; ------- get pixel size
	xstep = 0.06  &  ystep = 0.06	; default values
	test = tag_names(a.fastaligned.spec_aligned)
	for i = 0, n_elements(test)-1 do begin
		if test(i) EQ 'X_PIXEL_SIZE' then $
		     xstep = a.fastaligned.spec_aligned.x_pixel_size._data(0)
		if test(i) EQ 'Y_PIXEL_SIZE' then $
			ystep = a.fastaligned.spec_aligned.y_pixel_size._data(0)
	endfor

; ------ get dwell time
	dwell_time =  a.fastaligned.ExpTimes._data

; ------- get beam current
	ring_current=a.fastaligned.currents._data

; ----- get incident photon energy (eV) -----
	ev = a.fastaligned.energy._data

; ---- get data and normalize
	image_stack = a.fastaligned.spec_aligned._data
	for i = 0, npts-1 do begin
	  	image_stack(*,*,i) = image_stack(*,*,i)*norm/ring_current(i)
	ENDFOR
	 help, image_stack

  END

; XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
; -------- read in stack from type 1 hdf5 file

'one':	BEGIN 				; first format sent; derived from their tomo format
	npts = a.nxtomo.instrument.sample.data.number_of_frames._data
	nx = a.nxtomo.instrument.sample.data.image_width._data
	ny = a.nxtomo.instrument.sample.data.image_height._data

; ------- get pixel size
	xstep = a.NXTOMO.instrument.sample.x_pixel_size._data(0)
	ystep = a.NXTOMO.instrument.sample.y_pixel_size._data(0)

; -------- ring current
	ring_current =  A.NXTOMO.instrument.sample.current._data

; --------- read data; NB  in Alba images it is a bipolar integer so must add 2^32-1
	image_stack = a.nxtomo.instrument.sample.data._data  + 32767.
	for i = 0, npts-1 do begin
	  	image_stack(*,*,i) = image_stack(*,*,i)*norm/ring_current(i)
	endfor

; ------ get dwell for each image
	dwell_time =  A.NXTOMO.instrument.sample.ExpTimes._data

; ----- get incident photon energy (eV) -----
  	ev = A.NXTOMO.instrument.source.energy._data

END
ENDCASE

; ------------- prepare necessary parameters
n_cols=nx & n_rows=ny			; needed to pass elsewhere via bsif_com
n_data = npts
x_start = 0  & x_stop = nx*xstep
y_start = 0  & y_stop = ny*ystep
x = findgen(nx) & x = x_start +x*xstep
y = findgen(ny) & y = y_start +y*ystep

; ---------------------------------------
; START OF WRITE OUT OF binary stack file
;------------------------------------------

; -------- Create aXis2000 stack

filename_ev_msec_list = strarr(npts)
for i = 0, npts-1 do filename_ev_msec_list(i) = string(ev(i), format='(F8.2)') + ' eV      ' $
         + string(dwell_time(i), format='(F6.2)') + ' s'

; --------- ask user what file name they want - default name is same as files
fileout = file_path + fileshort + '.ncb'
sname=pickfile2(/write, file= fileout, filter='*.ncb', title = 'name of binary stack file')
stack_wb, sname, /silent

; -------- print tracing for trouble shooting
	IF NOT keyword_set(silent) then  begin
		axis_log, 'Abla.Mistral TXM stack read from ' + fileshort
		print, 'Read type ', type,' Alba - Mistral image from ', file
		print, 'Ring current(mA) range = ', ring_current(0), '  to  ', ring_current(npts-1)
		print, 'nx, ny, #images ', nx, ny, npts
		print, 'Dwell time (s) = ', dwell_time(0)
		print, 'Energy range (eV) = ', ev(0),'  to  ', ev(npts-1)
	ENDIF

; -------- launch stack_process ----------------
; stack_process, sname, /no_align, fpath=WritePath, /binary, /no_read

s=fileshort
return, s
endelse

END