; Copyright (c) 1998-2015 A.P. Hitchcock  All rights reserved
;+
;NAME:
;		READ_BESSY_STACK.PRO
;
;LAST CHANGED: ----------------------------------- 	13-Jan-15
;
; PURPOSE:
;	This function reads a sequence of images measured with Bessy STXM (Thieme)
; and converts the set to an  aXxis2000 *.ncb format
;
; CATEGORY:
;	data read routine. can execute outside of aXis2000
;
; CALLING SEQUENCE:
;	Result = READ_BESSY_STACK(directory=directory, filter=filter, normalize=normalize, silent=silent, txm=txm)
;
; INPUTS:
;
; KEYWORDS:
;	AXIS:	if set to 1, indicates called from aXis
;	DIRECTORY:  name of directory
;	FILTER: extension (default for images is '*.pgm')
;	NORMALIZE: if set, modify intensities to standard ring current (value of normalize parameter
;	SILENT: if set, do not print tracking comments
;	TXM:  if set, data is from the Bessy TXM
;
; OUTPUTS:
;	No explicit outputs.
;
; COMMON BLOCKS:
; AXIS_COM	standard set of common blocks
; STACK_PROCESS_COM
; BSIF_COM
; volume_data, image_stack
;
; SIDE EFFECTS: none
;
; RESTRICTIONS: none
;
; PROCEDURE:
;
; MODIFICATION HISTORY:
; (13-Oct-08 aph) first version
; (22-Sep-11 aph) extend to read in Bessy TXM stacks
; (29-Oct-11 aph) complete development from 22-sep-11
; (13-Jan-15 aph) adapt to take larger data sets; fix normalization
;-

function read_bessy_stack, directory=directory, filter=filter, normalize=normalize, $
       axis=axis, silent=silent, txm=txm
@axis_com
@bsif_com
@STACK_PROCESS_COM
COMMON volume_data, image_stack

on_error,2

; determine if AXIS is running (therefore may have called ax_svd )
; either when AXIS keyword is supplied or if any widget active
if  keyword_set(axis) then axis_on = 1 else axis_on = widget_info(/active)

check=0
if not keyword_set(filter) then begin
	if keyword_set(txm) then filter = '*.spe' else filter='*.pgm'
endif

if keyword_set(normalize) then ring_norm=normalize else ring_norm=300.
if keyword_set(txm) then begin
; ----- select files & set DefPath to this path --------------------------
	file=DIALOG_PICKFILE(/Read, /multiple, FILTER=filter)
	if strlen(file(1)) GT 0 THEN BEGIN
		name_bits = ax_name(file(0))	; sort the name into its parts
		path = name_bits(0)
		lpath= path
		defpath=path
		print, 'files: ', file
		file2 = file(1:n_elements(file)-1)
		file = [file2, file(0)]
		print, 'files: ', file
	 endif else begin
		axis_log, 'Directory ' + path +  ' is not defined, or does not exist'
		return, 0
	endelse
	bin = 2
	if axis_on then begin
		bin = get_num(prompt='bin factor ', val=bin, group = axis_ID)
	endif else begin
		bin = get_num(prompt='bin factor ', val=bin)
	endelse

; -------- read the first image to get dimensions --------
	n_img = n_elements(file)
	ev=fltarr(n_img)
	filename_ev_msec_list = strarr(n_img)
	tmp = read_bessy_txm(file=file(0), /silent, bin=bin, normalize = ring_norm)
	text_mask = '00'
	if n_img GE 100 then text_mask='000'
	if n_img GE 1000 then text_mask='0000'
	text_n_img = ' of '+ strtrim(string(n_img),2)
	if not keyword_set(silent) then begin
		text_i = strtrim(string(1),2)
		text = text_mask
		strput, text, text_i,strlen(text_i)
		axis_log, 'image '+ text + text_n_img
	endif

	t = ax_name(file(0))
	filename_ev_msec_list(0) = t(1)+ string(tmp.e, format='(G13.6)')
	ev(0)=tmp.e
	n_cols = n_elements(tmp.x)
	n_rows = n_elements(tmp.y)
	x_start = tmp.x(0)
	x_stop  = tmp.x(n_cols-1)
	y_start = tmp.y(0)
	y_stop  = tmp.y(n_rows-1)
	xstep = (x_stop - x_start)/n_cols
	ystep = (y_stop - y_start)/n_rows

;	help, tmp, /struct
	text = string(n_cols, n_rows, n_img, format='(3I5)')
	axis_log, 'Stack size is (col, row, e) ' + text
	mem = string((double(n_cols)*double(n_rows)*double(n_img))/1.e6, format='(F8.2)')
	axis_log, 'Estmated memory needed (32-bit words) Mb ' + mem
	image_stack=fltarr(n_cols, n_rows, n_img)


	image_stack(*,*,0) = tmp.d

; --------  read in the rest of the stack ------
	ax_wait
	for i=1, n_img-1 do begin
		tmp = read_bessy_txm(file=file(i), /silent,bin=bin, normalize = 300.)
		if not keyword_set(silent) then begin
			text_i = strtrim(string(i+1),2)
			text = text_mask
			strput, text, text_i,strlen(text_i)
			axis_log, 'image '+ text_i + text_n_img
		endif
		image_stack(*,*,i) = tmp.d
		ev(i) = tmp.e
		t = ax_name(file(i))
		filename_ev_msec_list(i) = t(1)+ string(tmp.e, format='(G13.6)')
	endfor
	text = string(n_cols, n_rows, n_img, format='(3I5)')
	axis_log, 'Stack dimension is (col, row, e) ' + text
	mem = string((double(n_cols)*double(n_rows)*double(n_img))/1.e6, format='(F8.2)')
	axis_log, 'Memory (32-bit words) Mb ' + mem

endif else begin
; ------------------- READ in old Bessy STXM stacks (Oct 2008) ---------
	if n_elements(directory) eq 0 then begin  ;popup file dialog box
	   file=PICKFILE2(/Multiple, FILTER=filter, /LPATH, DEFPATH=defpath)
	endif

	if strlen(file(0)) GT 0 THEN BEGIN
		name_bits = ax_name(file(0))	; sort the name into its parts
		path = name_bits(0)
		fileshort = name_bits(1)
	    ext = name_bits(2)
	 endif else begin
		axis_log, 'Directory ' + path +  ' is not defined, or does not exist'
		return, 0
	endelse

	if not keyword_set(silent) then axis_log, 'Reading Bessy STXM stack: ' + path

	; =================
	; construct a list of the files to read (all *.pgm in the directory)
	cd, path
	fnames=file_search(filter, count=a)
	if a LE 0 then begin
		axis_log, 'no *.pgm files found'
		return, check
	endif
	n_imgs = n_elements(fnames)
	filename_ev_msec_list = strarr(n_imgs)

	if not keyword_set(normalize) then begin
		norm=300.  ; normalize to standard start-of-fill value (or top-up)
		t=dialog_message('Normalize to ring current ?',/question)
		if t EQ 'Yes' then begin
			if n_elements(axis_ID) EQ 0 then norm = get_num(prompt='standard ring current', val=norm) $
			       else norm = get_num(prompt='standard ring current', val=norm, group = axis_ID)
		endif
	endif else norm = normalize
	if not keyword_set(silent) then axis_log, 'Normalize to ' + string(norm) + ' mA'
	; --------- read first image to get image size and other parameters
	ev=fltarr(n_imgs)
	image=read_bessy_stxm(file=fnames(0), /silent, normalize=norm)
	i=0
	t=ax_name(fnames(i))
	filename_ev_msec_list(i) = t(1)+ string(image.e, format='(G13.6)')
	ev(i)=image.e
	n_cols = n_elements(image.x)
	n_rows = n_elements(image.y)
	x_start = image.x(0)
	x_stop  = image.x(n_cols-1)
	y_start = image.y(0)
	y_stop  = image.y(n_rows-1)
	xstep = (x_stop - x_start)/n_cols
	ystep = (y_stop - y_start)/n_rows
	if not keyword_set(silent) then axis_log, string(i) + '   ' + filename_ev_msec_list(i)
	image_stack=fltarr(n_elements(image.x), n_elements(image.y), n_imgs)
	image_stack(*,*,i) = image.d

	; -------- now process the other images
	for i=1, n_imgs-1 do begin
		image=read_bessy_stxm(file=fnames(i), /silent, normalize=norm)
		ev(i)=image.e
		image_stack(*,*,i) = image.d
		t=ax_name(fnames(i))
		filename_ev_msec_list(i) = t(1)+ string(image.e, format='(G13.6)')
		if not keyword_set(silent) then axis_log, string(i) + '   ' + filename_ev_msec_list(i)
	endfor

ENDELSE


 ; return, image_stack
 ; end

; --------  write image_stack -----------------------
stackname = pickfile2(/read, prompt='Stack name', filter='*.ncb')
if strlen(stackname) NE 0 then begin
	tn=ax_name(stackname)
	stackname=tn(0)+tn(1)+'.ncb'
	stack_wb, stackname
	t=ax_name(stackname)
	if not keyword_set(silent) then axis_log, 'generated Bessy stxm stack ' + stackname
endif
check=1
return, check
end