; Copyright (c) 1998-2016 A.P. Hitchcock  All rights reserved
;+
;NAME:
;		READ_SDF
;
;LAST CHANGED: ----------------------------------- 11-Apr-16
;
; PURPOSE:
;	This function reads self-defining format ascii files written by "ALS 532"-class instruments
; using STXM_Control as their operating system (as of 30-Jul-09: ALS (3), CLS, SLS, Bessy, Stanford)
; Spectra, images, stacks and linescans have *.hdr and a set of *.xim files
; Images and stacks with non-square pixels are extended to have square pixels
; Image modes have separate file for each {energies.channels.regions}
; Spectra are a *#.xsp file for each region (#) containing multi-column spectra.
; Images, stacks are normalized to 400 mA ring current (N400)
;
; CATEGORY:
;	Input / output utilty; stand alone or from ax_sdf / axis2000
;
; CALLING SEQUENCE:
;	Result = READ_SDF(file,[ filter=filter, channel = channel, region = region, verbose=verbose, $
;             group=group, header_only=header_only, no_display=no_display, no_Save=no_save, map=map, $
;             no_interpolate=no_interpolate, one_image=one_image, _extra=e])
;
; INPUTS: none required
;
; KEYWORDS:
; 	FILE		optional file name
;	FILTER	   user defined filter
;	CHANNEL		preselected data channel
;	REGION		preselected region
;	ALL_REGIONS read in all regions
;	VERBOSE		print out details of parsing structure (testing)
;	GROUP		group leader (Axis_ID if called from aXis2000 and axis)
;	HEADER_ONLY extract and return header contents only
;	IMAGE_NUM	2 number aray defining background and signal for maps (first image is 1)
;	MAP			treat NEXAFS Image Scan as a set of 2 images to derive OD difference as a map
;	NO_DISPLAY	if set, do not store or display in curbuf
;	NO_INTERPOLATE  if set, do not convert rectangular to square pixels
;   NO_SAVE  	if set, do not save output
;   SILENT		do not print any tracking information
;   ONE_IMAGE	if defined, index of image to be read (first image is number 1)
;	_EXTRA		other passed through parameters
;
; OUTPUTS:
;	structure variable; either 1d (spectra) or 2d (image, linescan) or 3d (stack, as *.ncb)
;
; COMMON BLOCKS:
;	STACK_PROCESS_COM	common for stack_process
;	AXIS_COM	standard set of common blocks
;	BSIF_COM	stack common
;	VOLUME_DATA	stack data set
;	SDF_COM		common for SDF read in
;	SDFORMAT	reference number for naming muti-element arrays of structures
;
; MODIFICATION HISTORY:
; (11-feb-01 aph) first version
; (23-feb-01 aph) arrays of structures implemented
; (25-feb-01 aph) option to return structure only implemented
; (27-jun-01 aph) adapt to read focus and OSA scans
; (12-jul-01 aph) enable read-in of spectral modified format
; (29-jul-01 aph) retired read_bl5 spectral read-in; new format with *.hdr is in place
;               regions implemented for spectra
; (14-aug-01 aph) fixed errors found by Ivo in spectra & linescan read-in
; (17-aug-01 aph) adapt to REGIONS (kluge needed !)
; (22-aug-01 aph) stack readin (directly convert to *.ncb file)
; (09-nov-01 aph) activate all region read in for spectra
; (17-feb-02 aph) motor scan processing
; (28-feb-02 aph) ensure works for partial stacks
; (12-may-02 aph) correct header name format
; (28-may-02 aph) delete last column - set of zero's which may cause processing problems
; (11-dec-02 aph) adapt for VLM Image readin; (NB group_ID problem if run stand alone)
; (02-apr-03 aph) 'OSA focus scan' read-in
; (29-apr-03 aph) add tag  in log and axis windows to identify bad images
; (19-jun-04 aph) auto sort energies on readin of stack - NB PROBLEM with IDL sort function !!
; (22-aug-04 aph) replace 'GoTo' with  'ToGo' to avoid problem with IDL 6 throwing error on 'GoTo'
;                 kluge solution to out-of-order sort
; (29-aug-04 aph) restricted sort correction to IDL 6 using !version.release check
;        NB - array_equal not in pre-IDL6 ; sort still does not work correctly !
; (21-nov-04 aph) read Detector scans
; (20-feb-05 aph) add MAP function to read in a 2-image stack and take OD difference
;                 add one_image finction to read one image from a stack
; (09-mar-05 aph) fix IDL 5.2 partial stack read-in problem
; (24-mar-05 aph) simplify label for 1-image read-in (filename indicates it)
; (17-jun-05 aph) channels and regions working for map and 1-image read; default region = 1
; (24-aug-05 aph) extend to 26 channels of data
; (25-aug-05 aph) add Pattern and Chart read-in
; (15-dec-05 aph) adapt to non-axis and non-interactive use (for tomography)
;				  add definition of images for mapping; add keywords SILENT, NO_INTERPOLATE
;				  correct error on interpolation to square pixels
; (22-Mar-06 aph) Images, stacks normalized to 400 mA ring current (N400); add chart scan read in
; (04-May-06 aph) normalize stacks image-by-image to the ring current
; (11-mar-07 aph) switch on panel to set flag to control ring current normalization
; (18-May-07 aph) correct read-in of 1-image for multi-channel data (B. Watts)
; (14-Jun-07 aph) systematically close_lun, free_lun (replaced 2 close,/all commands) to allow stack read and release
; (14-Jul-07 aph) allow transfer of stack to stack_analyze even if do not write raw stack
; (13-Oct-08 aph) convert from stack-analyze to stack_process
; (01-Feb-09 aph) add interpolation of images to recorded (x,y) values (TRIGRID and 2-pass (1d line-by-line)
; (31-Jul-09 aph) clean-up transfer from stack-analyze to stack_process
;                 change sdf_normalize so it adapts to changes in norm_value
; (22-Jan-10 aph) remove print of group, axisID except when verbose is on
; (19-Feb-10 aph) add region name to label
; (16-Apr-10 aph) fix auto-read of channel & one-image for XRF stack read
; (17-Apr-10 aph) add no_display option
; (20-Apr-10 aph) fix no_display so it works
; (02-feb-11 aph) stacks with non-square pixels are converted (congrid) to square pixels
; (19-feb-11 aph) fix bad i-ring normalization in stxm_11 read-ins
; (01-Jun-11 aph) clean-up tracking
; (27-Jul-11 aph) remove old (commented out) use of read_bl5 (intermediate data format for STXM5322 in 2001)
; (31-jul-11 aph) converted map to use 1,2,3 (people) numbering instead of 0,1,2,3 (IDL)
;                 so that the one_image and map image numbering protocols are the same
; (08-May-13 aph) adapt read-in to changes Russ Berg made in file format for XMCD scans
; (28-Aug-13 aph) provide default file name for stacks when read in
; (09-Nov-13 aph) sdf_lastfile, sdf_path defined in main routine; add keyword no_save
; (11-Nov-13 aph) ignore request to normalize to ring current if beam current not defined or below a threshold
; (21-Dec-13 aph) reversed change in what sdf_lastfile is - was main name; had changed to full (path, main * extension)
;                which ended up having the read_sdf widget displaying full version in the name box
; (16-jun-14 aph) corrected error in Dec-13 version where sdf_lastfile was defined with the '.hdr' extension
; (21-Jul-14 aph) ensure stand alone operation works; variabilize the ring current normalization value
; (24-Oct-14 aph) modify map to deal with rectangular pixels
; (29-Nov-14 aph) introduce scan.type ='Ptycho Image Scan' images; add read_image_scan() function
; (02-Jul-15 aph) clean up feedback of ring current read-in - silent keyword replaced with
;                 SDF_silent common variable
;                 fix problem of slow read-in (pixel by pixel scaling !)
; (12-Jul-15 aph) corrected Pixel_Diff ==> Pix_diff (in sdf-com)
; (20-mar-16 aph) fixed error in setting norm_flag (found reading aXis_pass_tes_data\read\sdf\linescan)
; (11-Apr-16 aph) move make_Square_pixels (for images using rectangular ones) to ax_interp_square
;-

Function sdf_normalize,d, scandef, current=current
; --------------- normalize to  ring current to give response as per a pre-defined value(400 mA)
on_error,2
@stack_process_com
@axis_com
@bsif_com
@sdf_com

if norm_flag EQ 0 then begin		; only normalize if norm_flag is set
	norm_string = '   not ring current normalized'
	return,d
endif
norm = norm_value		; set in aXis2000 (31-jul-09 to 500 mA)
;                         - can change in ax_sdf widget (20-Jul-14)

if not keyword_set(current) then begin
	test = strupcase(tag_names(scandef))
	q = where(test EQ 'STORAGERINGCURRENT', count)
	if count GT 0 then begin
		Iring= ScanDef.STORAGERINGCURRENT
;		axis_log, 'Ring current (mA) ' + string(Iring, format='(F6.1)')
		scale = float(norm/Iring)
		norm_string = '  Normalized to ' + strtrim(fix(norm_value),2) + ' mA'
		if Iring LT 0.001 then begin
			axis_log, 'ring current < 1 uA; not ring current normalized'
			norm_string = '   not ring current normalized'
			scale = 1.
		endif
	endif else begin
		axis_log, '   not ring current normalized'
		scale=1.
		norm_string = '   not ring current normalized'
	endelse
endif else begin
	Iring = current
;	axis_log, 'Ring current (mA) ' + string(Iring, format='(F6.1)')
	norm_string = '  Normalized to ' + strtrim(fix(norm_value),2) + ' mA'
	scale = float(norm)/float(Iring)
	if current LE 0.001 then begin
		axis_log, 'ring current < 1 uA; not ring current normalized'
		norm_string = '   not ring current normalized'
		scale = 1.
	 endif
endelse
; normalize - OK, we now have the right scale factor
; NB Jul 20, 2014 - if start with norm_flag undefined
;           d*scale  ==> scalar. IDL should multiply array*scalar and end up with an array
; 02-jul-2015  - problem is that scale is a 1-element array and you need to identify element 0
; help, d
t = size(d)
if t(0) eq 2 then begin		; for 2D data
	x_ind = t(1)  & y_ind = t(2)
	d = d*scale(0)
endif else begin		; for 1D data
	x_ind = t(1)
	d = d*scale(0)
endelse

RETURN, d
END



; ==========================================================================

Function read_sdf_ax, ScanDef, region, y_axis=y_axis
on_error,2
@stack_process_com
@axis_com
@bsif_com
@sdf_com
;
; ----------------- extract appropriate axis from the header
if keyword_set(y_axis) then image = 1 else image = 0
CASE REGION OF
  0 : BEGIN
	x = float(ScanDef.PAxis_0.points)
	if image then y = float(ScanDef.QAxis_0.points)
	END
  1 : BEGIN
  	x = float(ScanDef.PAxis_1.points)
	if image then y = float(ScanDef.QAxis_1.points)
	END
  2 : BEGIN
  	x = float(ScanDef.PAxis_2.points)
	if image then y = float(ScanDef.QAxis_2.points)
	END
  3 : BEGIN
  	x = float(ScanDef.PAxis_3.points)
	if image then y = float(ScanDef.QAxis_3.points)
	END
  4 : BEGIN
  	x = float(ScanDef.PAxis_4.points)
	if image then y = float(ScanDef.QAxis_4.points)
	END
  5 : BEGIN
  	x = float(ScanDef.PAxis_5.points)
	if image then y = float(ScanDef.QAxis_5.points)
	END
  6 : BEGIN
  	x = float(ScanDef.PAxis_6.points)
	if image then y = float(ScanDef.QAxis_6.points)
	END
  7 : BEGIN
  	x = float(ScanDef.PAxis_7.points)
	if image then y = float(ScanDef.QAxis_7.points)
	END
  8 : BEGIN
  	x = float(ScanDef.PAxis_8.points)
	if image then y = float(ScanDef.QAxis_8.points)
	END
  9 : BEGIN
  	x = float(ScanDef.PAxis_9.points)
	if image then y = float(ScanDef.QAxis_9.points)
	END
  ELSE : BEGIN
  	axis_log, 'Region ' + strtrim(string(fix(region)),2) + $
  	          ' exceeds maximum software set up to read (9)'
  ENDELSE
ENDCASE
if image then return, y else return, x
END

; ---------------------------------------------------------------------------
; ---------  used for reading both Image Scan and Ptycho Image scan --------
; --- first generated 29-Nov-14 (aph).  NOT ALL OPTIONS TESTED !!
; ---------------------------------------------------------------------------

function read_image_scan, file, filter=filter, channel=channel, $
		 n_regions = n_regions, let = let, scandef = scandef, fileshort = fileshort, $		; line added to get it to work 29-nov-14 (aph)
		 region = region, all_regions = all_regions, $
         verbose=verbose, group=group, header_only=header_only, map=map, $
         image_num = image_num, one_image = one_image, no_display=no_display, $
         no_interpolate=no_interpolate, no_save=no_save, xy_correct=xy_correct,_extra=e
on_error,2
@stack_process_com
@axis_com
@bsif_com
@sdf_com

COMMON volume_data, image_stack
common sdformat, starnum

if keyword_set(all_regions) then begin
	n_read = n_regions & region = -1
endif else begin
	n_read = 1 & region = region - 1
endelse
CurBuf = CurBuf - 1
start_file = file
for j_region = 0, n_read-1 DO BEGIN
	CurBuf = CurBuf + 1 & if CurBuf EQ 10 then CurBuf = 0
	region = region + 1
;	--- determine file name -----------------
	t = ax_name(start_file)
	if n_regions gt 1 then begin
;				print, ' identifier for CLS-STXM XMCD 1-E image scans (Flags = : Multi-Region Image Stack', ScanDef.Flags
		if ScanDef.Flags EQ 'Multi-Region Image Stack' then $
			img_file = t(0) + t(1) + '_' + let + '000' + strtrim(string(fix(region)), 2) + '.xim' else $
			img_file = t(0) + t(1) + '_' + let + strtrim(string(fix(region)), 2) + '.xim'
	endif else begin
		img_file = t(0) + t(1) + '_' + let +'.xim'
	endelse
	d = text_read1(img_file)
;	------ pack into AXIS 2d structure -----------------
	x = read_sdf_ax(ScanDef, region)
	y = read_sdf_ax(ScanDef, region, /y_axis)

	nx = n_elements(x) & ny = n_elements(y)
;			print, nx, ny
	energy = ScanDef.StackAxis.min

; -------- xy_correct using channels (b,c) for (xpos, ypos) data ; 01-feb-09 aph)
;			print, 'xy_correct_Flag = ',xy_correct_flag
	if xy_correct_flag EQ 1 then begin
		x_file = t(0) + t(1) + '_b.xim'
		xa = text_read1(x_file)
		y_file = t(0) + t(1) + '_c.xim'
		ya = text_read1(y_file)
; select method
  		xy_method = dialog_message('2d Trigrid (No), or x, then y (Yes)',/question, $
  			    Title='Method to correct (xy) values')
  		If xy_method EQ 'Yes' then begin
  				check_spline=dialog_message('Linear (Yes) or Spline (No) ',/question)
  				if check_spline EQ 'No' then  xy_spline = 1 else xy_spline = 0
  			for i=1, ny-1 do begin
  				d(*,i) = interpol(d(*,i), xa(*,i), x, spline=xy_spline)
  			endfor
  			check_y=dialog_message('Also correct y-values ? ',/question)
  			If check_y EQ 'Yes' then begin
	  			for j=0, nx-1 do begin
	  				d(j,*) = interpol(d(j,*), ya(j,*), y, spline=xy_spline)
	  			endfor
  			endif
  		endif else begin	 ; Trigrid - both (x,y) used simultaneously
			mxa=max(abs(xa)) & mya=max(abs(ya))
			tol = 1.0e-9*max(mxa,mya)
			if keyword_set(group) then $
	  			tol = get_num(prompt='triangulate tolerance', val=tol, group=group) $
	   			else tol = get_num(prompt='triangulate tolerance', val=tol)
			triangulate, xa, ya, tr, tol=tol
			d= trigrid(xa,ya,d,tr,xout=x,yout=y)
			if not SDF_silent then axis_log, 'xy_corrected with tol = '+ $
			   strtrim(string(tol,format='(G13.2)'),2)
		endelse
	endif

; -------- normalize to ring current --------------------
	if norm_flag then d = sdf_normalize(d, scandef)
	dl = fileshort + '   ' + SDF_Channel_List(fix(n_data)-1)
	dwell = ScanDef.ImageScan.EnergyRegions[0].DwellTime
	xl = string(FORMAT='("x (um)     E = ",f8.3," eV     dwell = ",f8.2," ms")', $
 			 energy, dwell)
	t = ax_name(img_file)
	fileshort = t(1)
	if n_regions GT 1 then dl = dl + ' region ' + strtrim(string(region+1),2)
    s = {t:'2d', x:x, y:y, d:d, e: energy, xl:xl, yl: '  um     ch '+let + norm_string(0), dl:dl}

; --- check if rectangular pixels - if so interpolate to square
; ----- the make_square_pixel code was moved to ax_interp_square  (11-Apr-16 aph)
	if NOT keyword_set(no_interpolate) then s = ax_interp_square(s)

; ----- document read in of an image; store in axis2000 CurBuf if active and requested
	text = string(format= $
       '(i4," x",i4," image read from ",A,". E =",f7.2," eV. Dwell =",f8.2," ms.")', $
        nx, ny, fileshort, energy, dwell)
	if not SDF_silent then axis_log, text
	if  n_elements(axis_ID) EQ 1 and axis_ID NE 0 then begin
		if not keyword_set(no_display) then begin
			HANDLE_VALUE, Data(Curbuf), s, /SET
	 		Plotbuf, CurBuf
	 	endif
 	endif
endfor
return, s
END

function read_sdf, file, filter=filter, channel=channel, $
		 region = region, all_regions = all_regions, $
         verbose=verbose, group=group, header_only=header_only, map=map, $
         image_num = image_num, one_image = one_image, no_display=no_display, $
         no_interpolate=no_interpolate, no_save=no_save, xy_correct=xy_correct,_extra=e
on_error,2
@stack_process_com
@axis_com
@bsif_com
@sdf_com

COMMON volume_data, image_stack
common sdformat, starnum

;------------- establish print feedback level for all READ_SDF routines  -----------------
;
; if keyword_set(silent) then SDF_silent=1 else SDF_silent=0
; NOW set in ax-sdf
   print, 'SDF_silent = ',SDF_silent

; ------ minimum X- Y-pixel size difference (fraction of the pixel size)
;        to trigger remeshing to square pixels
   pix_diff = 0.005

; ----------------- need a default value for region --------- (17-jun-05 aph)
if n_elements(region) EQ 0 then region = 1
; print, 'Region is ', region

if NOT keyword_set(group) then group = 0 else begin
	if n_elements(axis_ID) EQ 1 and axis_ID NE 0 then group = axis_ID
 endelse
if keyword_set(verbose) then print, 'group, axis_ID', group, axis_ID

; the preceding is a kluge !!
; want to make the group be flexible - if just called from ax_sdf widget, group is sdf_ID, otherwise
; should be axis_ID.

; if run outside of aXis2000 then need to set xy_correct_flag
 if  n_elements(axis_ID) EQ 0 then begin
	xy_correct_flag=0
	if keyword_set(xy_correct) then xy_correct_flag=1	; just use the flag, so can use common or keyword
endif

t = size(starnum)
if t(1) EQ 0 then starnum = 0   ; start of unique names - SHOULD USE A CHECK for existence

if n_elements(file) eq 0 then begin  ;popup file dialog box
   if not keyword_set(filter) then fltr='*.hdr' else fltr = filter
   file=PICKFILE2(/Read, FILTER=fltr, /LPATH, DEFPATH=defpath)
endif
s = ''
if strlen(file) LE 0 THEN return, s  ; bail-out if no filename
; ------- Open the header to be read -------
; first force to be a header
t = ax_name(file)
fileshort = t(1)
sdf_lastfile = t(1)
file_hdr = t(0) + t(1) + '.hdr'
t = findfile(file_hdr)
IF strlen(t(0)) EQ 0 then begin
     axis_log, file_hdr  + ' - header file not found'
	   return, s
ENDIF

; -----------
on_ioerror,nofile
get_lun,lun
openr,lun,file_hdr
; ------- Read header file into a single string (max = 32600 characters)  --------
line = ' '
all = ''
count = 0
while not  eof(lun) do begin
	readf,lun, line
;		print, line				; debug to ensure full file read-in
; replace reserved keyword 'GoTo'  with 'ToGo'
	t = strpos(line,'GoTo')
	if t GE 0 then 	StrPut, line, 'ToGo', t
	all = all + line
	count = count + 1			; the number of lines is used in spectral readin
endwhile
close, lun & free_lun, lun
on_ioerror,null


; format for recursive structure parsing
all = '{'+all+'};'
; return, all

; --------- parse the structure --------------
if keyword_set(verbose) then verbose = 1 else verbose = 0

reg = strpos(all,'Regions')
if reg GT 0 then begin
	b = strmid(all,0, reg-1) + '}'
	ScanDef = parse_sdf(b,verbose=verbose)
	ScanDef = ScanDef.ScanDefinition
;	return, ScanDef
	q = strpos(ScanDef.Flags,'Image')
	if q GE 0 then image = 1 else image = 0	; FLAG to differentiate spectra and images
	if image then stax = strpos(all,'StackAxis') $
	   else stax = strpos(all,'Channels')
; -- isolate regions information
	regions = strmid(all, reg, stax-reg)
	t = strpos(regions,'(') & tc = strpos(regions,',')
	n_regions = fix(strmid(regions,t+1,tc-t-1))
	for i = 0, n_regions-1 do begin
		pax = strpos(regions,'PAxis') & paxe = strpos(regions,'};',pax)
		paxis = strmid(regions,pax,paxe-pax)
		p = parse_sdf('{'+paxis,verbose=verbose)
		pn = 'PAxis_'+strtrim(string(i),2)
		ScanDef = create_struct(ScanDef,pn,p.paxis)
		regions = strmid(regions,paxe)
		if image then begin	; point scans do not have QAxis
			qax = strpos(regions,'QAxis') & qaxe = strpos(regions,'};',qax)
			qaxis = strmid(regions,qax,qaxe-qax)
			q = parse_sdf('{'+qaxis,verbose=verbose)
			qn = 'QAxis_'+strtrim(string(i),2)
			ScanDef = create_struct(ScanDef,qn,q.qaxis)
			regions = strmid(regions,qaxe)
		endif
	endfor
	ScanDef = create_Struct(ScanDef, 'N_regions', n_regions)
	end_part = '{' + strmid(all,stax-1)
	ScanDef_end1 = parse_sdf(end_part,verbose=verbose)
	end1 = tag_names(ScanDef_end1)
	q = where(end1 EQ 'CHANNELS', count)
	if count GT 0 then begin
;		print, ' Channels replaced by Channel_labels'
		end1(q(0)) = 'CHANNEL_LABELS'
		for i = 0, n_elements(end1)-1 do begin
			tl = end1(i)
			ScanDef = create_struct(ScanDef,tl, scandef_end1.(i))
		endfor
	endif else begin
		for i = 0, n_tags(ScanDef_end1)-1 do $
		ScanDef = create_Struct(ScanDef, end1(i), ScanDef_end1.(i))
	endelse
	SD_names = tag_names(ScanDef)
	if verbose then for i = 0, n_tags(ScanDef)-1 do print, SD_names(i)
	end_part = '{END= {' + strmid(all,stax-1)
	ScanDef_end2 = parse_sdf(end_part,verbose=verbose)
	end2 = tag_names(ScanDef_end2)
	q = where(end2 EQ 'BEAMFEEDBACK', count)
	if count GT 0 then begin
		if verbose then for i = q(0), n_elements(end2)-1 do print, end2(i)
		for i = q(0)+1, n_tags(ScanDef_end2)-1 do $
		ScanDef = create_Struct(ScanDef, end2(i), ScanDef_end2.(i))
	endif
; return, ScanDef
; ------ handle old type files
endif else  ScanDef = parse_sdf(all)
; help, /struct, ScanDef

;------------ kluge to adapt for reading in VLM-image format -----
	if n_tags(ScanDef) EQ 1 then ScanDef = ScanDef.ScanDefinition

if keyword_set(header_only) then begin
; ------- print header information
;	XDisplayFile,file_hdr
	return, ScanDef
endif

; ------------- read in VLM Image as a special case - *.hdr files is not structured right
if ScanDef.Type EQ 'VLM Image' then begin
	start_file = file
;	--- determine file name -----------------
	t = ax_name(start_file)
	img_file = t(0) + t(1) + '_a.xim'
	d = text_read1(img_file)
;	------ pack into AXIS 2d structure -----------------
	xmin = ScanDef.PAxis.min
	xmax = ScanDef.PAxis.max
	nx   = ScanDef.PAxis.nPoints
	xstep = (xmax - xmin)/nx
	x = findgen(nx)*xstep + xmin
	ymin = ScanDef.QAxis.min
	ymax = ScanDef.QAxis.max
	ny   = ScanDef.QAxis.nPoints
	ystep = (ymax - ymin)/ny
	y = findgen(ny)*ystep + ymin
; --- check if rectangular pixels - if so interpolate to square
	if NOT keyword_set(no_interpolate) then begin
		dx = x(1) - x(0)
		dy = y(1) - y(0)
		pix_avg =moment([dx,dy])
		if abs(dx - dy)/pix_avg(0) GT pix_diff then begin
			if dx GT dy then begin
				axis_log, ' Rectangular pixels: x interpolated'
				nx = fix((x(n_elements(x)-1) - x(0))/dy)
				x = interpol(x,nx)
			endif else begin
				axis_log, ' Rectangular pixels: y interpolated'
				ny = fix((y(n_elements(y)-1) - y(0))/dx)
				y = interpol(y,ny)
			endelse
			nx = n_elements(x) & ny = n_elements(y)
			td = congrid(d,nx,ny, /interp, cubic=-0.5)
			d = td
		endif
	endif
	nx = n_elements(x) & ny = n_elements(y)
	energy = 2.
	dwell = ScanDef.Dwell
	xl = string(FORMAT='("x (um)     E = ",f8.3," eV     dwell = ",f6.2," ms")', $
 			 energy, dwell)
	t = ax_name(img_file)
	fileshort = t(1)
	dl = fileshort
    s = {t:'2d', x:x, y:y, d:d, e: energy, xl:xl, yl:'um', dl:dl}
	text = string(format= $
       '(i4," x",i4," image read from ",A,". E =",f7.2," eV. Dwell =",f6.2," ms.")', $
        nx, ny, img_file, energy, dwell)
	if not SDF_silent then axis_log, text
	if  n_elements(axis_ID) EQ 1 and axis_ID NE 0 then begin
		HANDLE_VALUE, Data(Curbuf), s, /SET
 		Plotbuf, CurBuf
 	endif
	return, s
endif

; ------- set up list of channels
n_channels = n_elements(ScanDef.Channel_labels)
if n_channels EQ 0 then n_channels = n_elements(ScanDef.ScanDefinition.channels)
if keyword_set(verbose) then axis_log, 'STXM data from '+ fileshort
if keyword_set(verbose) then axis_log,  strtrim(string(n_channels),2)+ ' channels'
SDF_Channel_List = strarr(n_channels)
for i = 0, n_channels-1 do SDF_Channel_List(i) = ScanDef.channel_labels[i].name

; identify channel which is to be read
if NOT keyword_set(channel) then begin
	for i = 0, n_channels-1 do axis_log, string(i+1) + '  ' + SDF_Channel_List(i)
	n_data =  0
	while n_data LT 1 OR n_data GT n_channels do begin
		if keyword_set(group) then $
		  n_data = get_num(prompt='Choose data channel', val=n_data, group=group) $
		   else n_data = get_num(prompt='Choose data channel', val=n_data)
	endwhile
endif else 	n_data = channel

IF NOT SDF_silent then print, 'channel is ', n_data
case n_data of
	1 : let = 'a'
	2 : let = 'b'
	3 : let = 'c'
	4 : let = 'd'
	5 : let = 'e'
	6 : let = 'f'
	7 : let = 'g'
	8 : let = 'h'
	9 : let = 'i'
	10 : let = 'j'
	11 : let = 'k'
	12 : let = 'l'
	13 : let = 'm'
	14 : let = 'n'
	15 : let = 'o'
	16 : let = 'p'
	17 : let = 'q'
	18 : let = 'r'
	19 : let = 's'
	20 : let = 't'
	21 : let = 'u'
	22 : let = 'v'
	23 : let = 'w'
	24 : let = 'x'
	25 : let = 'y'
	26 : let = 'z'    ; pronounced 'zee'
else : axis_log, 'Only set up for 26 channels max'
endcase

;print, 'Reading channel ' + strtrim(string(fix(n_data)),2) + ' letter ' + let
IF NOT SDF_silent then axis_log, 'Reading channel ' + SDF_Channel_List(fix(n_data)-1)
; -------  identify number of regions and allow user to select a region
;t = where(tag_names(ScanDef.ScanDefinition) EQ 'REGIONS', count)
; if count GT 0 then begin		; only execute region dialog if scandefinition contains regions
if  NOT keyword_set(region) then begin
	n_regions = fix(ScanDef.n_regions)
	if n_regions GT 1 and not keyword_set(all_regions) then begin
			axis_log, strtrim(string(n_regions),2) + ' regions'
		if keyword_set(group) then $
		  region = get_num(prompt='Choose region', val=region, group=group) $
			   else region = get_num(prompt='Choose region', val=region)
	endif else region = 1
endif
region = region - 1

if keyword_set(all_regions) then begin
	region = 0
;	axis_log, ' Reading all regions'
endif else begin
	text = strtrim(string(fix(region)),2)
;	axis_log, ' Reading region ' + text
endelse

; ------- get axes for desired region
x = read_sdf_ax(ScanDef, region)
if image then y = read_sdf_ax(ScanDef, region, /y_axis)
nx = n_elements(x) & ny = n_elements(y)

; print, 'reading in region ', string(region+1, format='(I2)'),'.    ', nx, '  x  ', ny, '  points '
;return, ScanDef

; --------- Read in data according to the type of scan -----

TYPE = ScanDef.type
; print, 'reading self_defining_data file of type ', type

CASE TYPE of
   'Image Scan': BEGIN
	   s = read_image_scan(file, filter=filter, channel=channel, $
		 n_regions = n_regions, let = let, scandef = scandef, fileshort = fileshort, $
		 region = region, all_regions = all_regions,$
         verbose=verbose, group=group, header_only=header_only, map=map, $
         image_num = image_num, one_image = one_image, no_display=no_display, $
         no_interpolate=no_interpolate, no_save=no_save, xy_correct=xy_correct,_extra=e)
	   return, s
	END

   'Ptycho Image Scan' : BEGIN
	   s = read_image_scan(file, filter=filter, channel=channel, $
		 n_regions = n_regions, let = let, scandef = scandef, fileshort = fileshort, $
		 region = region, all_regions = all_regions, $
         verbose=verbose, group=group, header_only=header_only, map=map, $
         image_num = image_num, one_image = one_image, no_display=no_display, $
         no_interpolate=no_interpolate, no_save=no_save, xy_correct=xy_correct,_extra=e)
       s.dl = 'Ptycho STXM image  ' + s.dl
	   return, s
	END


   'Focus Scan' : BEGIN
		t = ax_name(file)
		img_file = t(0) + t(1) + '_' + let +'.xim'
		d = text_read1(img_file)
;	------ pack into AXIS 2d structure -----------------
		energy = ScanDef.StackAxis.min
		dwell = ScanDef.FocusScan.DwellTime
		xl = string(FORMAT='("(um)     E = ",f8.3," eV     dwell = ",f6.2," ms")', $
     			 energy, dwell)
		t = ax_name(file)
		fileshort = t(1)
	s = {t:'2d', x:x, y:y, d:d, e: energy, xl:xl, yl:'Zone Plate Z (um)  ch '+let, dl: fileshort}
		text = string(format= $
	  '(i4," x",i4," focus scan read from ",A,". E =",f7.2," eV. Dwell =",f6.2," ms.")', $
	   nx, ny, file, energy, dwell)
	if not SDF_silent then axis_log, text
	return, s
	END

	'Detector Scan' : BEGIN
		t = ax_name(file)
		img_file = t(0) + t(1) + '_' + let +'.xim'
		d = text_read1(img_file)
;	------ pack into AXIS 2d structure -----------------
		energy = ScanDef.StackAxis.min
		dwell = ScanDef.DetScan.DwellTime
		xl = string(FORMAT='("X (um)     E = ",f8.3," eV     dwell = ",f6.2," ms")', $
     			 energy, dwell)
		t = ax_name(file)
		fileshort = t(1)
	s = {t:'2d', x:x, y:y, d:d, e: energy, xl:xl, yl:'detector Y (um)  ch '+let, dl: fileshort}
		text = string(format= $
	  '(i4," x",i4," detector scan read from ",A,". E =",f7.2," eV. Dwell =",f6.2," ms.")', $
	   nx, ny, file, energy, dwell)
	if not SDF_silent then axis_log, text
	return, s
	END

	'OSA Scan' : BEGIN
		t = ax_name(file)
		img_file = t(0) + t(1) + '_' + let +'.xim'
		d = text_read1(img_file)
;	------ pack into AXIS 2d structure -----------------
		energy = ScanDef.StackAxis.min
		dwell = ScanDef.OSAScan.DwellTime
		xl = string(FORMAT='("OSA x (um)     E = ",f8.3," eV     dwell = ",f6.2," ms")', $
     			 energy, dwell)
		t = ax_name(file)
		fileshort = t(1)
	s = {t:'2d', x:x, y:y, d:d, e: energy, xl:xl, yl:'OSA y (um)  ch '+let, dl: fileshort}
		text = string(format= $
	  '(i4," x",i4," OSA scan read from ",A,". E =",f7.2," eV. Dwell =",f6.2," ms.")', $
	   nx, ny, file, energy, dwell)
	if not SDF_silent then axis_log, text
	return, s
	END

	'OSA Focus Scan' : BEGIN
		t = ax_name(file)
		img_file = t(0) + t(1) + '_' + let +'.xim'
		d = text_read1(img_file)
;	------ pack into AXIS 2d structure -----------------
		energy = ScanDef.StackAxis.min
		dwell = ScanDef.OSAFocusScan.DwellTime
		xl = string(FORMAT='("OSA-x (um)     E = ",f8.3," eV     dwell = ",f6.2," ms")', $
     			 energy, dwell)
		t = ax_name(file)
		fileshort = t(1)
	s = {t:'2d', x:x, y:y, d:d, e: energy, xl:xl, yl:'Zone Plate Z (um)  ch '+let, dl: fileshort}
		text = string(format= $
	  '(i4," x",i4," OSA focus scan read from ",A,". E =",f7.2," eV. Dwell =",f6.2," ms.")', $
	   nx, ny, file, energy, dwell)
	if not SDF_silent then axis_log, text
	return, s
	END


	'XANES': BEGIN
		test = ''
		energy = float(x)
		npts = n_elements(energy)
;		print, 'E-axis', npts, ' points from ', energy(0), ' to ', energy (n_elements(energy)-1), ' eV.'
		d = fltarr(npts)  & i = 0
		ldata = fltarr(n_channels+1)
;		Print, file, '       # points: ', strtrim(string(npts),2),'      # data channels: ', strtrim(string(n_channels),2)
		dwell_time = ScanDef.ImageScan.EnergyRegions[0].DwellTime
		if dwell_time EQ 0 then dwell_time = 1.

		if keyword_set(all_regions) then begin
			n_read = n_regions & region = -1
		endif else begin
			n_read = 1 & region = region - 1
		endelse
		CurBuf = CurBuf - 1
		start_file = file
		for j_region = 0, n_read-1 DO BEGIN
			CurBuf = CurBuf + 1 & if CurBuf EQ 9 then CurBuf = 0
			region = region + 1
;	--- determine file name -----------------
			t = ax_name(start_file)
			if n_regions gt 1 then begin
				file = t(0) + t(1) + '_' + strtrim((string(fix(region), format='(A)')),2) + '.xsp'
			endif else begin
				file = t(0) + t(1) + '.xsp'
				test = findfile(file) & test = test(0)
			    if test EQ '' then file = t(0) + t(1) + '_0.xsp'
			endelse
			if not SDF_silent then print, 'Testing file name: ', file
			test = findfile(file) & test = test(0)
			if not SDF_silent then print, 'Result = ', test
			if test EQ '' then $
			  file = pickfile2(title = 'Select spectral file name', FILTER='*.xsp', /LPATH, DEFPATH=defpath)
			if file EQ '' then return, ''
			if not SDF_silent then print, 'opening data file: ', file
;			Print, file, '       # points: ', strtrim(string(npts),2),'      # data channels: ', strtrim(string(n_channels),2)
			openr,lun,file, /get_lun
			for i = 0, npts-1 do begin
				readf,lun, line
	;					print, 'line ', i,': ', line
				reads, line, ldata
				energy(i) = ldata(0)
				d(i)  = ldata(n_data)/dwell_time
			endfor
			close, lun & free_lun, lun
	;		plot, energy,d

	; --------------- generate axis 1d structure -----
			t = ax_name(file)
			fileshort = t(1)
			xl = string(FORMAT='("Energy(eV)    dwell = ",f7.2," ms")', dwell_time)
			yl = 'STXM SDF: CHANNEL ' + let + '.  kHz   )'
;		---------- normalize to ring current if in the header ------------------
			if norm_flag then d = sdf_normalize(d, scandef)
			dl = fileshort + '   ' + SDF_Channel_List(fix(n_data)-1)
			if n_regions GT 1 then dl = dl + ' region ' + strtrim(string(region+1),2)
			tmp = {t:'1d', x:energy, d:d, dn:d, xl:xl, yl:'kHz   ' +  norm_string, dl:dl}
			if not SDF_silent then axis_log, 'read SDF NEXAFS spectrum: ' + file
			npts = n_elements(tmp.x)
			text = string(format='(i4," pts. ",f9.4," to ",f9.4," eV. Dwell= ",f8.2," ms.")', $
			    npts, tmp.x(0), tmp.x(npts-1), dwell_time)
			if not SDF_silent then axis_log, text
			if  n_elements(axis_ID) EQ 1 and axis_ID NE 0 then begin
				HANDLE_VALUE, Data(Curbuf), tmp, /SET
		 		Plotbuf, CurBuf
		 	endif
		endfor
		return, tmp
	END


	'NEXAFS Point Scan': BEGIN
	; ------- open file for selected region - each region is in separate file
		test = ''
		energy = float(x)
		npts = n_elements(energy)
;		print, 'E-axis', npts, ' points from ', energy(0), ' to ', energy (n_elements(energy)-1), ' eV.'
		d = fltarr(npts)  & i = 0
		ldata = fltarr(n_channels+1)
		dwell_time = ScanDef.ImageScan.EnergyRegions[0].DwellTime
		if dwell_time EQ 0 then dwell_time = 1.
		if keyword_set(all_regions) then begin
			n_read = n_regions & region = -1
		endif else begin
			n_read = 1 & region = region - 1
		endelse
		CurBuf = CurBuf - 1
		start_file = file
		for j_region = 0, n_read-1 DO BEGIN
			CurBuf = CurBuf + 1 & if CurBuf EQ 10 then CurBuf = 0
			region = region + 1
;	--- determine file name -----------------
			t = ax_name(start_file)
			if n_regions gt 1 then begin
				file = t(0) + t(1) + '_' + strtrim((string(fix(region), format='(A)')),2) + '.xsp'
			endif else begin
				file = t(0) + t(1) + '.xsp'
				test = findfile(file) & test = test(0)
			    if test EQ '' then file = t(0) + t(1) + '_0.xsp'
			endelse
			test = findfile(file) & test = test(0)
			if test EQ '' then $
			  file = pickfile2(title = 'Select spectral file name', FILTER='*.xsp', /LPATH, DEFPATH=defpath)
			if file EQ '' then return, ''
;			print, 'opening data file: ', file
;			Print, file, '       # points: ', strtrim(string(npts),2),'      # data channels: ', strtrim(string(n_channels),2)
			openr,lun,file, /get_lun
			for i = 0, npts-1 do begin
				readf,lun, line
				reads, line, ldata
				energy(i) = ldata(0)
				d(i)  = ldata(n_data)/dwell_time
			endfor
			close, lun & free_lun, lun

	; --------------- generate axis 1d structure -----
			t = ax_name(file)
			fileshort = t(1)
			xl = string(FORMAT='("point", I2, "  Energy(eV)    dwell = ",f7.2," ms")', n_data+1, dwell_time)
			yl = 'STXM: CHANNEL ' + let + '.  kHz   '
;		---------- normalize to ring current if in the header ------------------
			if norm_flag then d = sdf_normalize(d, scandef)
			dl = fileshort + '   ' + SDF_Channel_List(fix(n_data)-1)
			if n_regions GT 1 then dl = dl + ' region ' + strtrim(string(region+1),2)
			tmp = {t:'1d', x:energy, d:d, dn:d, xl:xl, yl:yl + norm_string, dl:dl}
			if not SDF_silent then axis_log, 'read SDF NEXAFS spectrum: ' + fileshort
			npts = n_elements(tmp.x)
			reg_txt=region + 1
			text = string(format='("Region ",i2," ",i4," pts. ",f9.4," to ",f9.4," eV. Dwell= ",f7.2," ms.")', $
			    reg_txt, npts, tmp.x(0), tmp.x(npts-1), dwell_time)
			if not SDF_silent then axis_log, text

			if  n_elements(axis_ID) EQ 1 and axis_ID NE 0 then begin
				HANDLE_VALUE, Data(Curbuf), tmp, /SET
		 		Plotbuf, CurBuf
		 	endif
		endfor
		return, tmp
	END


; ------------------------------ STACK readin ------------------------------------
	'NEXAFS Image Scan': BEGIN
		if keyword_Set(map) then goto, stack_map
		if keyword_Set(one_image) then goto, stack_image
		if not SDF_silent then axis_log, 'Reading in STXM Stack ' + fileshort
		dwell_time = ScanDef.ImageScan.EnergyRegions[0].DwellTime
		energy = float(ScanDef.StackAxis.points)
		npts = n_elements(energy)
		nx = n_elements(x) & ny = n_elements(y)
		image_stack = fltarr(nx,ny,npts)
		filename_ev_msec_list = strarr(npts)
		t = ax_name(file) & sourcepath = t(0) & filebase = t(1)
        widget_control,hourglass=1
		for i = 0, npts-1 do begin
			zero ='000'
			t_num = strtrim(strcompress(string(i)),2)
			strput, zero, t_num, 3-strlen(t_num)
			if n_regions gt 1 then begin
				file_short = filebase +'_' + let + zero + strtrim(string(fix(region)),2) + '.xim'
			endif else begin
				file_short = filebase +'_' + let + zero + '.xim'
			endelse
			filenm = sourcepath + file_short
	; ----- check for file existence (recovery of part stacks)
			test = findfile(filenm)
			if strlen(test(0)) EQ 0 then goto, part_stack
;			text = string(i) + '  ' + file_short			; give feedback on progress
;			if not SDF_silent then axis_log, text
			d = text_read1(filenm)
; ------------ normalize to ring current  - skip if norm_flag = 0  ------
			if norm_flag then begin
				iring = sdf_stack_cur(i)
				if iring LE 1 then iring = 1.				; (11-Nov-13 do not normalize if iring LE 1
				d = sdf_normalize(d, scandef, current=iring)
			endif
; --------------
			image_stack(*,*,i) = d
			filename_ev_msec_list(i) = file_short + ' ' + $
			    string(energy(i), format = '(F7.2)') + '  ' +  string(dwell_time, format='(F5.2)')

; --------- feedback on the image that was read
			text = string(i, format='(I3)') + '   E (eV) = ' + string(energy(i), format = '(F7.2)')
			if norm_flag then text =  text + '  I(mA) = ' + strtrim(string(Iring, format = '(F6.1)'),2)
			axis_log, text
		endfor

		goto, full_stack
		part_stack:
			t = image_stack
			text = string(format='("Incomplete NEXAFS Image scan. Only ", i3, " of ", i3," images exist.")', i, npts)
			axis_log, text
			npts = i
			image_stack = fltarr(nx,ny,npts)
			image_stack = t(*,*,0:npts-1)
			t = filename_ev_msec_list
			filename_ev_msec_list = strarr(npts)
			filename_ev_msec_list = t(0:npts-1)
			if not SDF_silent then axis_log, 'info on images ' + filename_ev_msec_list
		full_stack:
; define parameters needed for stack write out
		ev = energy
		svec = size(image_stack)
		n_cols = nx  & n_rows = ny
		x_start = x(0)
		x_stop = x(n_elements(x)-1)
		y_start = y(0)
		y_stop = y(n_elements(y)-1)
		xstep = (x_stop - x_start)/n_cols
		ystep = (y_stop - y_start)/n_rows

; --- check if rectangular pixels - if so interpolate to square
		if NOT keyword_set(no_interpolate) then begin
			dx = x(1) - x(0)
			dy = y(1) - y(0)
			pix_avg =moment([dx,dy])
			if abs(dx - dy)/pix_avg(0) GT pix_diff then begin
				if dx GT dy then begin
					 axis_log, ' Rectangular pixels: x interpolated'
					nx = fix((x(n_elements(x)-1) - x(0))/dy)
					x = interpol(x,nx)
				endif else begin
					axis_log, ' Rectangular pixels: y interpolated'
					ny = fix((y(n_elements(y)-1) - y(0))/dx)
					y = interpol(y,ny)
				endelse
				nx = n_elements(x) & ny = n_elements(y)
				image_stack2 = fltarr(nx,ny,npts)
				for i = 0, npts-1 do begin
					d = reform(image_stack(*,*,i))
					td = congrid(d, nx, ny, /interp, cubic=-0.5)
					image_stack2(*,*,i) = td
				endfor
				image_stack = image_stack2
			endif
		endif
; ------------ sort based on ENERGY if needed ----------  (19jun04 aph) - AUTO
	goto, skipE		; for now lets ignore !!
		sort_ev = sort(ev)             ;(29-aug-04 aph) adapt to run version)
	; ----------- ERROR in sort - puts first as last element ---------
		last = n_elements(ev)-1
		te = sort_ev(last)								; KLUGE !!
		for i=last-1, 0,-1 do sort_ev(i+1) = sort_ev(i)	; KLUGE !!
		sort_ev(0) = te									; KLUGE !!
		test = indgen(n_elements(ev))
;		print, 'sort_ev ', sort_ev
;		print, 'test ', test
; --------- array_equal only available in IDL 6 -------
		if !version.release GE 6.0 then same_order = array_equal(sort_ev, test) $
		      else same_order = 1
		if same_order NE 1 then begin
;			print, 'original order', ev
			axis_log, 'Energies in stack are out-of-order. Sorting .....'
			timg = image_stack
			tlbl = filename_ev_msec_list
			for i = 0,n_elements(ev)-1 do begin
				image_stack(*,*,i) = timg(*,*,sort_ev(i))
				filename_ev_msec_list(i) = tlbl(sort_ev(i))
				ev(i) = ev(sort_ev(i))
			endfor
	;		ev = ev(sort(ev))
	; KLUGE due to error [in IDL 6.0 sort ??]  - last element not correct !!!!!
			if !version.release GE 6.0 then begin
				print, 'KLUGE!:  transferring last element to first'
				ti = image_stack(*,*,last)
				tl = filename_ev_msec_list(last)
				te = ev(last)
				for i = last-1,0,-1 do begin
					image_stack(*,*,i+1) = image_stack(*,*,i)
					filename_ev_msec_list(i+1) = filename_ev_msec_list(i)
					ev(i+1) = ev(i)
				endfor
				image_stack(*,*,0) = ti
				filename_ev_msec_list(0) = tl
				ev(0) = te
			endif
;			print, 'sorted order', ev
		endif
		skipE:
; ------- set default name to the file number ---------
		filename = string(fix(strmid(fileshort, strlen(fileshort)-3)))
		if not keyword_set(no_save) then begin
		    filename = pickfile2(filter = '*.ncb', file = filename, title = 'binary stack file ')
	        if strlen(filename) GT 0 then begin
				t = ax_name(filename)
	        	filename = t(0) + t(1) +'.ncb'
	        	stack_wb, filename
	        endif else return, filename      ; filename='stack.ncb' PRIOR to 28-Aug-13, if canceled, gave problems
	        if not keyword_set(no_display) then Stack_process, filename, /no_align, /binary
        endif
        return, filename

; ***************************  read-in 1 2-image stack as a MAP *********
stack_map:
		text = 'Reading in STXM 2-image map' + file
		if not SDF_silent then axis_log, text
		energy = ScanDef.StackAxis.points
		npts = n_elements(energy)
		nx = n_elements(x) & ny = n_elements(y)
		if npts GT 2 then begin
			if not SDF_silent then axis_log,  'warning: stack map has more than 2 energies'
		endif
		t = ax_name(file) & sourcepath = t(0) & filebase = t(1)

; ------ determine image number0 from pre-defined input, or assume it is 1, 2 (** changed 29-Jul-11)
		if not keyword_set(image_num) then image_num = [1,2]
; ------ read in reference ----------
		zero ='000'
		t = string(fix(image_num(0)-1))
		t_num0 = strtrim(strcompress(t),2)
		strput, zero, t_num0, 3-strlen(t_num0)
		if n_regions gt 1 then begin
			file_short = filebase +'_' + let + zero + strtrim(string(fix(region)),2) + '.xim'
		endif else begin
			file_short = filebase +'_' + let + zero + '.xim'
		endelse
		filenm = sourcepath + file_short
		test = findfile(filenm)				; ----- check for file existence
		if strlen(test(0)) EQ 0 then return,0
		if not SDF_silent then print, 'reading in background file: ', filenm
		ref = text_read1(filenm)


; ------ read in signal ----------
		zero ='000'
		t = string(fix(image_num(1)-1))
		t_num1 = strtrim(strcompress(t),2)
		strput, zero, t_num1, 3-strlen(t_num1)
		if n_regions gt 1 then begin
			file_short = filebase +'_' + let + zero + strtrim(string(fix(region)),2) + '.xim'
		endif else begin
			file_short = filebase +'_' + let + zero + '.xim'
		endelse
		filenm = sourcepath + file_short
		test = findfile(filenm)				; ----- check for file existence
		if strlen(test(0)) EQ 0 then return,0
		if not SDF_silent then print, 'reading in signal file ', filenm
		signal = text_read1(filenm)

; ------- assume first image is reference, second image is signal ---------
		Io_ref = max(ref)
;		Io_ref = get_num(prompt = 'Io level for reference',val = Io_ref, group=axis_ID)
		ref = -1.*alog(float(ref)/Io_ref)
		Io_sig = max(signal)
;		Io_sig = get_num(prompt = 'Io level for signal',val = Io_sig, group=axis_ID)
		signal = -1.*alog(float(signal)/Io_sig)
		signal = signal - ref

; -------- return the OD- difference as an axis image -------
;	------ pack into AXIS 2d structure -----------------
		dwell = ScanDef.ImageScan.EnergyRegions[0].DwellTime
		E_val=energy(t_num1)	; assumes second E is more interesting one
		xl = string(FORMAT='("x (um)     E = ",f8.3," eV     dwell = ",f6.2," ms")', $
     			 E_val, dwell)
		dl = fileshort + ' Stack map '
		if n_regions GT 1 then dl = dl + ' region ' + strtrim(string(region+1),2)
; --- check if rectangular pixels - if so interpolate to square  aph added 23-oct-14
		if NOT keyword_set(no_interpolate) then begin
			dx = x(1) - x(0)
			dy = y(1) - y(0)
			pix_avg =moment([dx,dy])
			if abs(dx - dy)/pix_avg(0) GT pix_diff then begin
				if dx GT dy then begin
					axis_log, ' Rectangular pixels: x interpolated'
					nx = fix((x(n_elements(x)-1) - x(0))/dy)+1
					x = interpol(x,nx)
				endif else begin
					axis_log, ' Rectangular pixels: y interpolated'
					ny = fix((y(n_elements(y)-1) - y(0))/dx)+1
					y = interpol(y,ny)
				endelse
				nx = n_elements(x) & ny = n_elements(y)
				td = congrid(signal,nx,ny, /interp, cubic=-0.5)
				signal = td
			endif
		endif
		nx = n_elements(x) & ny = n_elements(y)
		dwell = ScanDef.ImageScan.EnergyRegions[0].DwellTime
		xl = string(FORMAT='("x (um)     E = ",f8.3," eV     dwell = ",f6.2," ms")', $
     			 E_val, dwell)
; -------- pack into image structure
	    s = {t:'2d', x:x, y:y, d:signal, e: E_val, xl:xl, yl:'um  ch '+let, dl:dl}
		if  n_elements(axis_ID) EQ 1 and axis_ID NE 0 then begin
			HANDLE_VALUE, Data(Curbuf), s, /SET
	 		Plotbuf, CurBuf
	 	endif
	 	return, s

; ***************************  read-in 1 image  from a stack  *********
stack_image:
		energy = ScanDef.StackAxis.points
		text = 'STXM: reading in a single image from a stack, E = ' $
		     + string(strtrim(string(energy(one_image-1),format='(f7.2)'),2))
		if NOT SDF_silent then axis_log, text

		if keyword_set(all_regions) then begin
			n_read = n_regions & region = -1
		endif else begin
			n_read = 1 & region = region - 1
		endelse
		CurBuf = CurBuf - 1
; -------- define image file name
		t = ax_name(file) & sourcepath = t(0) & filebase = t(1)
		zero ='000'
		t_num = strtrim(strcompress(string(fix(one_image-1))),2)	; subtract 1 as incremented to allow index=0
		strput, zero, t_num, 3-strlen(t_num)
; -------- read in multiple regions if they exist
		for j_region = 0, n_read-1 DO BEGIN
			CurBuf = CurBuf + 1 & if CurBuf EQ 10 then CurBuf = 0
			region = region + 1
;	--- determine file name -----------------
		if n_regions gt 1 then begin
			file_short = filebase +'_' + let + zero + strtrim(string(fix(region)),2) + '.xim'
		endif else begin
			file_short = filebase +'_' + let + zero + '.xim'
		endelse
		filenm = sourcepath + file_short
		test = findfile(filenm)				; ----- check for file existence
		if strlen(test(0)) EQ 0 then return,0
		if not SDF_silent then print, 'reading in signal file ', filenm
		d = text_read1(filenm)

;	------ pack into AXIS 2d structure -----------------
		x = read_sdf_ax(ScanDef, region)
		y = read_sdf_ax(ScanDef, region, /y_axis)
; --- check if rectangular pixels - if so interpolate to square
		if NOT keyword_set(no_interpolate) then begin
			dx = x(1) - x(0)
			dy = y(1) - y(0)
			pix_avg =moment([dx,dy])
			if abs(dx - dy)/pix_avg(0) GT pix_diff then begin
				if dx GT dy then begin
					axis_log, ' Rectangular pixels: x interpolated'
					nx = fix((x(n_elements(x)-1) - x(0))/dy)+1
					x = interpol(x,nx)
				endif else begin
					axis_log, ' Rectangular pixels: y interpolated'
					ny = fix((y(n_elements(y)-1) - y(0))/dx)+1
					y = interpol(y,ny)
				endelse
				nx = n_elements(x) & ny = n_elements(y)
				td = congrid(d,nx,ny, /interp, cubic=-0.5)
				d = td
			endif
		endif
		nx = n_elements(x) & ny = n_elements(y)
		E_val = energy(one_image-1)
		dwell = ScanDef.ImageScan.EnergyRegions[0].DwellTime
		xl = string(FORMAT='("x (um)     E = ",f8.3," eV     dwell = ",f6.2," ms")', $
     			 E_val, dwell)
;		---------- normalize to ring current if in the header ------------------
		if norm_flag then d = sdf_normalize(d, scandef)
		dl = fileshort
		if n_regions GT 1 then dl = dl + ' region ' + strtrim(string(region+1),2)
	    s = {t:'2d', x:x, y:y, d:d, e: E_val, xl:xl, yl:norm_string + '   um     ch '+let, dl:dl}
		text = string(format= $
	       '(i4," x",i4," image read from ",A,". E =",f7.2," eV. Dwell =",f6.2," ms.")', $
	        nx, ny, fileshort, E_val, dwell)
		if NOT SDF_silent then axis_log, text
		if n_elements(axis_ID) EQ 1 and axis_ID NE 0 then begin
			if not keyword_set(no_display) then begin
				HANDLE_VALUE, Data(Curbuf), s, /SET
		 		Plotbuf, CurBuf
		 	endif
	 	endif
	endfor
	return, s
	END

	'NEXAFS Line Scan': BEGIN
		dwell = ScanDef.ImageScan.EnergyRegions[0].DwellTime
		t = ax_name(file)
		img_file = t(0) + t(1) + '_' + let +'.xim'
		d = float(text_read1(img_file))
		xl = string(FORMAT='("x (um)   dwell = ",f6.2," ms")', $
     			 dwell)
; ------ remove column at highest energy ----
		nx = n_elements(x)-1
     	x = x(0:nx-1)
     	dt = d
     	d = fltarr(nx,n_elements(y))
     	d = dt(0:nx-1,*)
;	---------- normalize to ring current if in the header ------------------
		if norm_flag then d = sdf_normalize(d, scandef)
		dl = fileshort
		s = {t:'2d', x:x, y:y, d:d, e: 0, xl:xl, yl:norm_string + '   um     ch '+let, dl: dl}
		text = string(format= $
	  	    '(" Linescan read from ",A,". Emin =",f7.2," eV. Emax =",f7.2," eV.  Dwell =",f6.2," ms.")', $
	   		   file, x(0), x(n_elements(x)-1), dwell)
		if NOT SDF_silent then axis_log, text
		return,s
	END


	'Pattern Gen': BEGIN        ; added 25-Aug-05
		nx = n_elements(x)
		ny = n_elements(y)   ; as of 25-aug-05 Qaxis is just 1 point (undefined meaning)
		energy = ScanDef.StackAxis.min
		dwell = ScanDef.Dwell
		xl = string(FORMAT='("X (um)     E = ",f8.3," eV     dwell = ",f6.2," ms")', $
     			 energy, dwell)
		t = ax_name(file)
		img_file = t(0) + t(1) + '_' + let +'.xim'
;		y = fltarr(2)
;		y(0)=0.
;		y(1)= 10.
;		ny = n_elements(y)
		d = float(text_read1(img_file))
		s = {t:'1d', x:x, y:y, d:d, e: energy, xl:xl, yl:'Y (um)  ch '+let, dl: 'PATTERN ' + fileshort}
		text = string(format= $
		  '(i4," x",i4," pattern_gen image read from ",A,". E =",f7.2," eV. Dwell =",f6.2," ms.")', $
		   nx, ny, file, energy, dwell)
		if NOT SDF_silent then axis_log, text
		return, s
	END

	'Chart': BEGIN		; added 22-Mar-06
		dwell_time = ScanDef.Dwell
; ------ read in data from *.xsp file - potentially multi-channel
		npts = n_elements(x)
		d = fltarr(npts)  & i = 0
		ldata = fltarr(n_channels+1)
;		Print, file, '       # points: ', strtrim(string(npts),2),'      # data channels: ', strtrim(string(n_channels),2)
		if dwell_time EQ 0 then dwell_time = 1.
		n_read = 1 & region = region - 1
		start_file = file
;	--- determine file name -----------------
		t = ax_name(start_file)
		if n_regions gt 1 then begin
			file = t(0) + t(1) + '_' + strtrim((string(fix(region), format='(A)')),2) + '.xsp'
		endif else begin
			file = t(0) + t(1) + '.xsp'
			test = findfile(file) & test = test(0)
		    if test EQ '' then file = t(0) + t(1) + '_0.xsp'
		endelse
		test = findfile(file) & test = test(0)
		if test EQ '' then $
		  file = pickfile2(title = 'Select chart file name', FILTER='*.xsp', /LPATH, DEFPATH=defpath)
		if file EQ '' then return, ''
;		print, 'opening data file: ', file
;			Print, file, '       # points: ', strtrim(string(npts),2),'      # data channels: ', strtrim(string(n_channels),2)
		openr,lun,file, /get_lun
		for i = 0, npts-1 do begin
			readf,lun, line
			reads, line, ldata
			x(i) = ldata(0)
			d(i)  = ldata(n_data)/dwell_time
		endfor
		close, lun & free_lun, lun
; --------- prepare structure for returning data
		t = ax_name(file)
		fileshort = t(1)
		xl = string(FORMAT='("    dwell = ",f7.2," ms")', dwell_time)
		yl = 'STXM: CHANNEL ' + let + '.  counts per msec (rate in kHz)'
		dl = fileshort + 'Chart '
		if n_regions GT 1 then dl = dl + ' region ' + strtrim(string(region+1),2)
		tmp = {t:'1d', x:x, d:d, dn:d, xl:xl, yl:yl, dl:dl}
		if not SDF_silent then axis_log, 'read STXM Chart Scan: ' + fileshort
		npts = n_elements(tmp.x)
		text = string(format='(i4," pts. ",g13.6," to ",g13.6," microns. Dwell = ",f7.2," ms.")', $
		    npts, tmp.x(0), tmp.x(npts-1), dwell_time)
		if NOT SDF_silent then axis_log, text
		if  n_elements(axis_ID) EQ 1 and axis_ID NE 0 then begin
			HANDLE_VALUE, Data(Curbuf), tmp, /SET
		 	Plotbuf, CurBuf
		ENDIF
		return, tmp
	END


	; -------- Motor Scan or bailout -------------------------
	ELSE: BEGIN
	; assume this is a motor scan - check for existence of ScanDef.MotorScan tag
	; NB (29-Nov-14) this is an efficient way to handle Motor Scans since the ScsnDef.type
	; variable does not say 'Motor Scan' but rather indicates the name of the motor scanned
		test = where(tag_names(ScanDef)	EQ 'MotorScan', count)
		if count EQ -1 then begin
			axis_log, 'Unknown scan type'
			return, scandef
		endif
		if not SDF_silent then axis_log, 'Motor Scan: ' + ScanDef.MotorScan.Motor
		dwell_time = ScanDef.Dwell
; ------ read in data from *.xsp file - potentially multi-channel
		npts = n_elements(x)
		d = fltarr(npts)  & i = 0
		ldata = fltarr(n_channels+1)
;		axis_log, file, '       # points: ', strtrim(string(npts),2),'      # data channels: ', strtrim(string(n_channels),2)
		if dwell_time EQ 0 then dwell_time = 1.
		n_read = 1 & region = region - 1
		start_file = file
;	--- determine file name -----------------
		t = ax_name(start_file)
		if n_regions gt 1 then begin
			file = t(0) + t(1) + '_' + strtrim((string(fix(region), format='(A)')),2) + '.xsp'
		endif else begin
			file = t(0) + t(1) + '.xsp'
			test = findfile(file) & test = test(0)
		    if test EQ '' then file = t(0) + t(1) + '_0.xsp'
		endelse
		test = findfile(file) & test = test(0)
		if test EQ '' then $
		  file = pickfile2(title = 'Select spectral file name', FILTER='*.xsp', /LPATH, DEFPATH=defpath)
		if file EQ '' then return, ''
;		axis_log, 'opening data file: ', file
;			axis_log, file, '       # points: ', strtrim(string(npts),2),'      # data channels: ', strtrim(string(n_channels),2)
		openr,lun,file, /get_lun
		for i = 0, npts-1 do begin
			readf,lun, line
			reads, line, ldata
			x(i) = ldata(0)
			d(i)  = ldata(n_data)/dwell_time
		endfor
		close, lun & free_lun, lun
; --------- prepare structure for returning data
		t = ax_name(file)
		fileshort = t(1)
		xl = ScanDef.MotorScan.Motor + string(FORMAT='("    dwell = ",f7.2," ms")', dwell_time)
		yl = 'STXM: CHANNEL ' + let + '.  counts per msec (rate in kHz)'
		dl = fileshort
		if n_regions GT 1 then dl = dl + ' region ' + strtrim(string(region+1),2)
		tmp = {t:'1d', x:x, d:d, dn:d, xl:xl, yl:yl, dl:dl}
		if not SDF_silent then axis_log, 'read STXM Motor Scan: ' + fileshort
		npts = n_elements(tmp.x)
		text = string(format='(i4," pts. ",g13.6," to ",g13.6," microns. Dwell = ",f7.2," ms.")', $
		    npts, tmp.x(0), tmp.x(npts-1), dwell_time)
		if NOT SDF_silent then axis_log, text
		if n_elements(axis_ID) EQ 1 and axis_ID NE 0 then begin
			HANDLE_VALUE, Data(Curbuf), tmp, /SET
		 	Plotbuf, CurBuf
		ENDIF
		return, tmp
	ENDELSE
ENDCASE
return, scandef

nofile:
axis_log, 'Open or read error on ' + file
close,/all
return, s
end

