; Copyright (c) 1998-2016 A.P. Hitchcock, B. Watts  All rights reserved
;+
;NAME:
; READ_NEXUS
;
;LAST CHANGED: ----------------------------------- 12-Apr-16 (aph)
;
;PURPOSE:
; This set of procedures is a widget to read in files
; in NXstxm, the HDF5 NeXus file format developed for STXM.
; The widget uses read_nexus to perform the file input.
; For more info on the NeXus file format, see www.nexusformat.org
;
; CATEGORY:
; Input / output utilty; stand alone or from ax_nexus / axis2000
;
; CALLING SEQUENCE:
; Result = READ_NEXUS(file,[ filter=filter, channel = channel, region = region, verbose=verbose, $
;             group=group, header_only=header_only, no_display=no_display, no_Save=no_save, $
;             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
; 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
; NO_DISPLAY  if set, do not store or display in curbuf
;   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
; NEXUS_COM   common for NeXus read in
;
;
;MODIFICATION HISTORY:
; (27-Apr-15 bw )  First version written using read_sdf.pro as a template
; (13-Jul-15 aph) first integration into axis2000 (02-July-15 version)
; (22-Jul-15 aph) add restriction to post-8.0 IDL versions
; (31-Jul-15 bw ) modified to use IDL6.3 routines; stack read-in corrected
; (14-Mar-16 bw ) improved with new NEXUS routines
; (27-Mar-16 aph) convert DOUBLES to FLOATs, to allow *.axb to be read
;                 get multi-region (multi-channel) working
; (04-Apr-16 aph) send back first image of stack if no_save is selected
; (05-Apr-16 aph) provide access to parameters ('Attributes')
; (12-Apr-16 aph) add ax_interp_square to properly display data with non-square pixels; veto this for OSA_focus and line scans
;-

PRO PUSHLIST, array, append
@nexus_com
on_error,2

     IF n_elements(array) EQ 0 THEN array=[append] ELSE array=[array,append]
END

PRO PUSHSTRUCT, array, append
@nexus_com
on_error,2

     IF n_elements(array) EQ 0 THEN array=CREATE_STRUCT(append) ELSE array=CREATE_STRUCT(array, append)
END

FUNCTION hdf5_object_type, Loc_id, name
@nexus_com
on_error,2

;Find the type of the referenced HDF5 object.
;Possible types are: 'GROUP', 'DATASET', 'TYPE' and 'UNKNOWN'
IF NOT keyword_set(name) THEN name='.'
test_result = H5G_GET_OBJINFO(Loc_id, name, /FOLLOW_LINK)
RETURN, test_result.TYPE
END

FUNCTION list_nexus_groups, Loc_id, name
@nexus_com
on_error,2

;List valid NeXus groups contained within the referenced group.
;Any HDF5 group having an NX_class attribute is considered valid.
IF NOT keyword_set(name) THEN name='.'
n_members = H5G_GET_NMEMBERS(Loc_id, name)
;list_of_nexus_groups = 0
FOR I = 0, n_members-1 DO BEGIN
 member_name = H5G_GET_MEMBER_NAME(Loc_id, name, I)
; print, "list", name+'/'+member_name
 IF hdf5_object_type(Loc_id, name+'/'+member_name) EQ 'GROUP' THEN BEGIN;NeXus groups
   M_id = H5G_OPEN(Loc_id, name+'/'+member_name)
   num_attr = H5A_GET_NUM_ATTRS(M_id)
   ClassAttr_id = H5A_OPEN_NAME(M_id, 'NX_class')
   member_class = H5A_READ(ClassAttr_id)
   H5A_CLOSE, ClassAttr_id
     new_group = create_struct('Name',member_name,'Class',member_class)
     PUSHLIST, list_of_nexus_groups,new_group
   H5G_CLOSE, M_id
 ENDIF
ENDFOR
RETURN, list_of_nexus_groups
END

FUNCTION get_nexus_attributes, Loc_id, name
@nexus_com
on_error,2

  IF hdf5_object_type(Loc_id, name) EQ 'DATASET' THEN BEGIN;NeXus fields and arrays
    D_id = H5D_OPEN(Loc_id, name)
    D_type = H5T_GET_CLASS(H5D_GET_TYPE(D_id))
    num_attr = H5A_GET_NUM_ATTRS(D_id)
    FOR J=0, num_attr-1 DO BEGIN
      A_id = H5A_OPEN_IDX(D_id, J)
      PUSHSTRUCT, A_struct, CREATE_STRUCT(H5A_GET_NAME(A_id), H5A_READ(A_id))
      H5A_CLOSE, A_id
    ENDFOR
    IF D_type EQ 'H5T_STRING' THEN String_Value = H5D_READ(D_id) ELSE String_Value = ''
    H5D_CLOSE, D_id
  ENDIF ELSE BEGIN ;NeXus groups
    D_id = H5G_OPEN(Loc_id, name)
    num_attr = H5A_GET_NUM_ATTRS(D_id)
    FOR J=0, num_attr-1 DO BEGIN
      A_id = H5A_OPEN_IDX(D_id, J)
      PUSHSTRUCT, A_struct, CREATE_STRUCT(H5A_GET_NAME(A_id), H5A_READ(A_id))
      H5A_CLOSE, A_id
    ENDFOR
    H5G_CLOSE, D_id
  ENDELSE
  IF num_attr EQ 0 THEN A_struct={dummy:0}
RETURN, A_struct
END

FUNCTION list_nexus_datasets, Loc_id, name
@nexus_com
on_error,2

;List HDF5 datasets contained within the referenced group and read strings along the way.
IF NOT keyword_set(name) THEN name='.'
n_members = H5G_GET_NMEMBERS(Loc_id, name)
;list_of_nexus_datasets = \NULL
FOR I = 0, n_members-1 DO BEGIN
  member_name = H5G_GET_MEMBER_NAME(Loc_id, name, I)
  IF hdf5_object_type(Loc_id, name+'/'+member_name) EQ 'DATASET' THEN BEGIN;NeXus groups
    D_id = H5D_OPEN(Loc_id, name+'/'+member_name)
    D_type = H5T_GET_CLASS(H5D_GET_TYPE(D_id))
    IF D_type EQ 'H5T_STRING' THEN String_Value = H5D_READ(D_id) ELSE String_Value = ''
    PUSHLIST, list_of_nexus_datasets, create_struct('Name',member_name,'Type',D_type,'String',String_Value)
    H5D_CLOSE, D_id
  ENDIF
ENDFOR
RETURN, list_of_nexus_datasets
END

FUNCTION get_nexus_regions, file
@nexus_com
on_error,2

file_id = H5F_OPEN(file)
group_list = list_nexus_groups(file_id, '/')
;region_list = []
FOR I=0, N_ELEMENTS(group_list)-1 DO BEGIN
  IF group_list[I].Class EQ 'NXentry' THEN BEGIN
    FLAG_definition_found = 0
    dataset_list = list_nexus_datasets(file_id, group_list[I].Name)
    FOR J=0, N_ELEMENTS(dataset_list)-1 DO BEGIN
;      PRINT, J, N_ELEMENTS(dataset_list), dataset_list[J]
      IF dataset_list[J].Name EQ 'definition' AND dataset_list[J].String EQ 'NXstxm' THEN BEGIN
        PUSHLIST, region_list,  '/'+group_list[I].Name
        FLAG_definition_found = 1
      ENDIF
    ENDFOR
   IF NOT FLAG_definition_found THEN BEGIN
     PRINT, 'Check for subgroups not yet implemented.'
   ENDIF
  ENDIF
ENDFOR
H5F_CLOSE, file_id
RETURN, region_list
END

FUNCTION get_nexus_region_info, file, region_path
@nexus_com
file_id = H5F_OPEN(file)
group_list = list_nexus_groups(file_id, region_path)
;channel_list = []
;energy_list = []
;scan_type = []
FOR I=0, N_ELEMENTS(group_list)-1 DO BEGIN
  IF group_list[I].Class EQ 'NXdata' THEN BEGIN
    PUSHLIST, channel_list, group_list[I].Name
  ENDIF
ENDFOR
Region_info = CREATE_STRUCT('Channels',channel_list)
dataset_list = list_nexus_datasets(file_id, region_path+'/'+channel_list(0))
FOR I=0, N_ELEMENTS(dataset_list)-1 DO BEGIN
  IF TOTAL(STRLOWCASE(dataset_list[I].Name) EQ ['e','energy','photon_energy']) THEN BEGIN
    D_id = H5D_OPEN(file_id, region_path+'/'+channel_list(0)+'/'+dataset_list[I].Name)
    energy_list = H5D_READ(D_id)
    H5D_CLOSE, D_id
    Region_info = CREATE_STRUCT(Region_info,'Energy',energy_list)
  ENDIF
  IF TOTAL(STRLOWCASE(dataset_list[I].Name) EQ ['stxm_scan_type','scan_type']) THEN BEGIN
    D_id = H5D_OPEN(file_id, region_path+'/'+channel_list(0)+'/'+dataset_list[I].Name)
    scan_type = H5D_READ(D_id)
    H5D_CLOSE, D_id
    Region_info = CREATE_STRUCT(Region_info,'Type',scan_type)
  ENDIF
ENDFOR
H5F_CLOSE, file_id
IF N_ELEMENTS(energy_list) GT 0 THEN BEGIN
  RETURN, Region_info
ENDIF ELSE BEGIN
  print, "Error: 'energy' not provided!"
ENDELSE
END

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
FUNCTION read_NXdata, file, region_path, channel_name, xy_correct=xy_correct, one_image=one_image
@axis_com
@bsif_com
@stack_process_com
@nexus_com
COMMON volume_data, image_stack

on_error,2

DATA_FLAG = 0

file_id = H5F_OPEN(file)
dataset_list = list_nexus_datasets(file_id, region_path+'/'+channel_name)
file_name_split = ax_name(file)
Data_File_Name = file_name_split(1)+'.'+file_name_split(2)
Attributes = get_nexus_attributes(file_id, region_path+'/'+channel_name)
Region_info = get_nexus_region_info(file, region_path)

IF TOTAL(TAG_NAMES(Region_info) EQ 'ENERGY') THEN BEGIN
  energy = Region_info.Energy
ENDIF ELSE BEGIN
  axis_log, "Error: 'energy' not provided!"
ENDELSE
IF TOTAL(STRLOWCASE(dataset_list.name) EQ 'sample_y') THEN BEGIN
  data_id = H5D_OPEN(file_id, region_path+'/'+channel_name+'/sample_y')
  sample_y = H5D_READ(data_id)
  H5D_CLOSE, data_id
ENDIF ELSE BEGIN
  axis_log, "Error: 'sample_y' not provided!"
ENDELSE
IF TOTAL(STRLOWCASE(dataset_list.name) EQ 'sample_x') THEN BEGIN
  data_id = H5D_OPEN(file_id, region_path+'/'+channel_name+'/sample_x')
  sample_x = H5D_READ(data_id)
  H5D_CLOSE, data_id
ENDIF ELSE BEGIN
  axis_log, "Error: 'sample_x' not provided!"
ENDELSE
IF TOTAL(STRLOWCASE(dataset_list.name) EQ 'count_time') THEN BEGIN
  data_id = H5D_OPEN(file_id, region_path+'/'+channel_name+'/count_time')
  count_time = H5D_READ(data_id)*1000; convert from seconds to milliseconds
  H5D_CLOSE, data_id
ENDIF ELSE BEGIN
  axis_log, "Error: 'count_time' not provided!"
ENDELSE
IF TOTAL(TAG_NAMES(Region_info) EQ 'TYPE') THEN BEGIN
  stxm_scan_type = Region_info.Type
ENDIF ELSE BEGIN
  axis_log, "Error: 'stxm_scan_type' not provided!"
ENDELSE


;print, "energy:", energy
;print, "count_time:", count_time
;print, "sample_y:", sample_y
;print, "sample_x:", sample_x

;Attributes attached to the NXdata group should describe where the data array is ('signal') and how its dimensions are labelled ('axes').
;Initial version were written another way, with attributes ('axis') attached to the fields within the NXdata group.
;First we try correct way, then fall back to methods that work for old, non-spec versions.
IF N_ELEMENTS(Attributes) NE 0 THEN BEGIN
  IF TOTAL(TAG_NAMES(Attributes) EQ 'AXES') THEN BEGIN
    i = 0
    FOR J=0, N_ELEMENTS(Attributes.axes)-1 DO BEGIN
      IF TOTAL(dataset_list.name EQ Attributes.axes[J]) THEN BEGIN
        data_id = H5D_OPEN(file_id, region_path+'/'+channel_name+'/'+Attributes.axes[J])
        Dim_Array = H5D_READ(data_id)
        ; The following two lines allow us to have a list of arrays of different sizes (it is actually a list of pointers)
        PUSHLIST, Axes_list, Ptr_New([0], /No_Copy) ; add pointer to list
        *Axes_list[i] = Dim_Array ; set pointer to data array
        H5D_CLOSE, data_id
        i = i+1
      ENDIF ELSE BEGIN
        print, "ERROR: dimension label not found: "+Attributes.axes[J]
      ENDELSE
    ENDFOR
  ENDIF ELSE BEGIN; Fallback to looking for "axis" attributes on each child field
    print, "'axes' attribute not found, will look for individual 'axis' attributes"
    i=0
    FOR J=0, N_ELEMENTS(dataset_list.name)-1 DO BEGIN; Go through set of fields and make a list of 'axis' attributes
      field_attr = get_nexus_attributes(file_id, region_path+'/'+channel_name+'/'+dataset_list[J].name)
      IF N_TAGS(field_attr) GT 0 THEN BEGIN
        IF TOTAL(TAG_NAMES(field_attr) EQ 'AXIS') THEN BEGIN
          IF i EQ 0 THEN BEGIN
            Axis_struct=CREATE_STRUCT('Axis', field_attr.axis, 'Target', dataset_list[J].name)
          ENDIF ELSE BEGIN
            Axis_struct=[Axis_struct,CREATE_STRUCT('Axis', field_attr.axis, 'Target', dataset_list[J].name)]
          ENDELSE
          i = i+1
        ENDIF
      ENDIF
    ENDFOR
;    print, Axis_struct
;    IF TOTAL(TAG_NAMES(Axis_struct) EQ 'Axis') THEN BEGIN
    IF N_ELEMENTS(Axis_struct) GT 0 THEN BEGIN
      FOR i=0, N_ELEMENTS(Axis_struct.Axis)-1 DO BEGIN; Go through list of 'axis' attributes and collect dimension labels
        FOR j=0, N_ELEMENTS(Axis_struct.Axis)-1 DO BEGIN
  ;        print, 'Compile Axes_list', i, Axis_struct[j].Axis, Axis_struct[j].Target
          IF Axis_struct[j].Axis EQ i+1 THEN BEGIN
            PUSHLIST, target_list, Axis_struct[j].Target
            data_id = H5D_OPEN(file_id, region_path+'/'+channel_name+'/'+Axis_struct[j].Target)
            Dim_Array = H5D_READ(data_id)
            ; The following two lines allow us to have a list of arrays of different sizes (it is actually a list of pointers)
            PUSHLIST, Axes_list, Ptr_New([0], /No_Copy) ; add pointer to list
            *Axes_list[i] = Dim_Array ; set pointer to data array
            H5D_CLOSE, data_id
            BREAK
          ENDIF
        ENDFOR
      ENDFOR
      Attributes = CREATE_STRUCT(Attributes, 'axes', target_list)
    ENDIF ELSE BEGIN
      print, "'axis' attributes not found either. This file is not written correctly."
    ENDELSE
  ENDELSE
  IF TOTAL(TAG_NAMES(Attributes) EQ 'SIGNAL') THEN BEGIN
    IF TOTAL(dataset_list.name EQ Attributes.signal) THEN BEGIN
;      print, "Signal points to real dataset "+Attributes.signal
      DATA_FLAG = 1
      data_id = H5D_OPEN(file_id, region_path+'/'+channel_name+'/'+Attributes.signal)
      Data_Array = H5D_READ(data_id)
      dataspace_id = H5D_GET_SPACE(data_id)
      H5S_CLOSE, dataspace_id
      H5D_CLOSE, data_id
    ENDIF
  ENDIF ELSE BEGIN; Fallback to assuming dataset is named "data"
    IF TOTAL(dataset_list.name EQ 'data') THEN BEGIN
      DATA_FLAG = 1
      data_id = H5D_OPEN(file_id, region_path+'/'+channel_name+'/data')
      Data_Array = H5D_READ(data_id)
      dataspace_id = H5D_GET_SPACE(data_id)
      H5S_CLOSE, dataspace_id
      H5D_CLOSE, data_id
    ENDIF ELSE BEGIN
      print, "ERROR: Dataset not found"
    ENDELSE
  ENDELSE
ENDIF
norm_string = ''
IF norm_flag EQ 1 THEN BEGIN
  norm_flag = 0
  group_list = list_nexus_groups(file_id, region_path)
  FOR G=0, N_ELEMENTS(group_list)-1 DO BEGIN
    IF group_list[G].Class EQ 'NXmonitor' AND group_list[G].Name EQ 'control' THEN BEGIN
      norm_flag = 1
      norm_string = 'N - ' + strtrim(fix(norm_value[0]),2) + ' mA'
      data_id = H5D_OPEN(file_id, region_path+'/control/data')
      Norm_Array = H5D_READ(data_id)
      mean_value = MEAN(Norm_Array, /NAN)
      index_NAN = where(NOT FINITE(Norm_Array), count_NAN)
      IF count_NAN NE 0 THEN Norm_Array[index_NAN] = mean_value
      index_0 = where(Norm_Array EQ 0, count_0)
      IF count_0 NE 0 THEN Norm_Array[index_0] = mean_value
      H5D_CLOSE, data_id
      IF DATA_FLAG EQ 1 THEN Data_Array = Data_Array*norm_value[0]/Norm_Array
      BREAK
    ENDIF
  ENDFOR
  IF NOT norm_flag THEN BEGIN
    print, "normalisation data 'control' not found."
  ENDIF
ENDIF




;IF keyword_set(xy_correct) THEN BEGIN
;  print, 'read_NXdata: XY correct =', xy_correct
;  group_list = list_nexus_groups(file_id, region_path)
;  print, group_list
;  FOREACH G, group_list DO BEGIN
;    IF G.Class EQ 'NXinstrument' THEN BEGIN
;      axis_log, 'Applying XY correction... (very slow)'
;      nx = n_elements(sample_x)
;      ny = n_elements(sample_y)
;      data_id = H5D_OPEN(file_id, region_path+'/'+G.Name+'/sample_x/data')
;      fine_X_points = REFORM(H5D_READ(data_id),nx,ny)
;      H5D_CLOSE, data_id
;      data_id = H5D_OPEN(file_id, region_path+'/'+G.Name+'/sample_y/data')
;      fine_Y_points = REFORM(H5D_READ(data_id),nx,ny)
;      H5D_CLOSE, data_id
;      data_id = H5D_OPEN(file_id, region_path+'/'+G.Name+'/'+channel_name+'/data')
;      Data_list = REFORM(H5D_READ(data_id),nx,ny)
;      H5D_CLOSE, data_id
;      Krig_switch=0
;      tmp = file_Search('krig2d_fast.pro',COUNT=Krig_switch)
;      IF N_ELEMENTS(Data_list) GT 500 THEN BEGIN
;       step = 20
;       overlap = 5
;       Data_Array[*,*]=0
;        FOR i=-overlap,nx-step/2-1,step/2 DO BEGIN
;          FOR j=-overlap,ny-step/2-1,step/2 DO BEGIN
;;            print, i,j, ' of ' ,nx,ny
;            boundary_values = [sample_x[MAX([0,i])],sample_y[MAX([0,j])],sample_x[MIN([(i+step),nx-1])],sample_y[MIN([(j+step),ny-1])]]
;            tx = MIN([(i+step),nx-1])-MAX([0,i])
;            ty = MIN([(j+step),ny-1])-MAX([0,j])
;            IF Krig_switch THEN BEGIN
;              Data_temp = KRIG2D_FAST( Data_list(MAX([0,i]):MIN([(i+step),nx-1]),MAX([0,j]):MIN([(j+step),ny-1])), fine_X_points(MAX([0,i]):MIN([(i+step),nx-1]),MAX([0,j]):MIN([(j+step),ny-1])), fine_Y_points(MAX([0,i]):MIN([(i+step),nx-1]),MAX([0,j]):
;            ENDIF ELSE BEGIN
;              Data_temp = KRIG2D( Data_list(MAX([0,i]):MIN([(i+step),nx-1]),MAX([0,j]):MIN([(j+step),ny-1])), fine_X_points(MAX([0,i]):MIN([(i+step),nx-1]),MAX([0,j]):MIN([(j+step),ny-1])), fine_Y_points(MAX([0,i]):MIN([(i+step),nx-1]),MAX([0,j]):MIN([
;            ENDELSE
;            Data_Array[MAX([0,i+overlap]):MAX([0,i+overlap])+tx-overlap,MAX([0,j+overlap]):MAX([0,j+overlap])+ty-overlap] = Data_temp[overlap-1:tx-1,overlap-1:ty-1]
;          ENDFOR
;        ENDFOR
;      ENDIF ELSE BEGIN
;        IF Krig_switch THEN BEGIN
;          Data_Array = KRIG2D_FAST( Data_list, fine_X_points, fine_Y_points, EXPONENTIAL=[2*(sample_x[nx-1]-sample_x[0])/nx,0.2,1], BOUNDS=boundary_values, NX=nx, NY=ny)
;        ENDIF ELSE BEGIN
;          Data_Array = KRIG2D( Data_list, fine_X_points, fine_Y_points, EXPONENTIAL=[2*(sample_x[nx-1]-sample_x[0])/nx,0.2,1], BOUNDS=boundary_values, NX=nx, NY=ny)
;        ENDELSE
;      ENDELSE
;      BREAK
;    ENDIF
;  ENDFOREACH
;ENDIF

IF DATA_FLAG THEN BEGIN
;  print, "Reading ", STRING(stxm_scan_type)+" with "+STRING(format='(i1)',SIZE(Data_Array, /N_DIMENSIONS))+ " dimensions."
  IF keyword_set(one_image) THEN BEGIN
    print, ' (single image #',strtrim(fix(one_image),2),'/',strtrim(fix(n_elements(energy),2)),')'
    energy = energy[one_image-1]
    count_time = count_time[one_image-1]
    Data_Array = Data_Array[*,*,one_image-1]
    Axes_list = Axes_list[1:*]
    stxm_scan_type = 'image slice'
  ENDIF
  CASE SIZE(Data_Array, /N_DIMENSIONS) OF
    1:BEGIN
;        print, '(1D scan).'
        n_e = n_elements(energy)
        log_text = string(format='(A," with ",i5," points read from ",A,". Dwell =",f6.2," ms.")',STRUPCASE(STRMID(stxm_scan_type,0,1))+STRMID(stxm_scan_type,1), n_e, Data_File_Name, count_time[0])
        if NOT keyword_set(silent) then axis_log, log_text
        X_string = string(FORMAT='("Point", i2, "  Energy (eV)    Dwell = ",f7.2," ms")', n_e, count_time[0])
        Y_string = norm_string+'    (counts)'
        RETURN_STRUCTURE = {t:'1d', x:*Axes_list[0], d:Data_Array, xl:X_string, yl:Y_string, dl:Data_File_Name+'   '+channel_name}
      END
    2:BEGIN
;        print, '(2D scan)'
        nx = n_elements(*Axes_list[1]) & ny = n_elements(*Axes_list[0])
        log_text = string(format='(A," (",i4," x ",i4," ) read from ",A,". Dwell =",f6.2," ms.")',STRUPCASE(STRMID(stxm_scan_type,0,1))+STRMID(stxm_scan_type,1), ny, nx, Data_File_Name, count_time[0])
        IF NOT keyword_set(silent) THEN axis_log, log_text
        IF n_elements(energy) GT 1 THEN BEGIN ; if linescan
          	E_string = string(FORMAT='("Energy (eV)    Dwell = ",f7.2," ms")', count_time[0])
          	L_string = norm_string+'   '+Attributes.axes[1]+' (um)'
          	tmp = {t:'2d', x:float(*Axes_list[0]), y:float(*Axes_list[1]), $
               	d:TRANSPOSE(float(Data_Array)), e:float(energy), xl:E_string, yl:L_string, dl:Data_File_Name+'   '+channel_name}
        ENDIF ELSE BEGIN ; if image
          	X_string = Attributes.axes[1]+string(FORMAT='(" (um)     E = ",f8.3," eV     dwell = ",f6.2," ms")', energy, count_time)
          	Y_string = norm_string+'   '+Attributes.axes[0]+' (um)'
          	tmp =  {t:'2d', x:float(*Axes_list[1]), y:float(*Axes_list[0]), d:float(Data_Array), e:float(energy(0)), xl:X_string, yl:Y_string, dl:Data_File_Name+'   '+channel_name}
        ENDELSE
  ;      print, 'nexus_head.Type = ', nexus_head.Type
       if nexus_head.Type NE 'osa focus' AND  nexus_head.Type NE 'sample line spectrum' then RETURN_STRUCTURE = ax_interp_square(tmp) else RETURN_STRUCTURE = tmp
      END
    3:BEGIN
;        print, '(3D scan)'
        nx = n_elements(*Axes_list[2]) & ny = n_elements(*Axes_list[1]) & n_e = n_elements(*Axes_list[0])
        log_text = string(format='(i3," x",i3," stack with ",i3," energies read from ",A,".")',nx, ny, n_e, Data_File_Name)
        IF NOT keyword_set(silent) THEN axis_log, log_text
        IF Nexus_no_save EQ 0 THEN BEGIN
            filename = pickfile2(filter='*.ncb', lfile=file_name_split(1), title='binary stack file ', WRITE=write)
            IF strlen(filename) GT 0 THEN BEGIN
              t = ax_name(filename)
              filename = t(0) + t(1) +'.ncb'
              image_stack = Data_Array
              x_start = (*Axes_list[2])[0]
              x_stop =  (*Axes_list[2])[N_ELEMENTS(*Axes_list[2])-1]
              y_start = (*Axes_list[1])[0]
              y_stop =  (*Axes_list[1])[N_ELEMENTS(*Axes_list[1])-1]
              ev = *Axes_list[0]
              filename_ev_msec_list = strarr(n_elements(ev))
              FOR i = 0, n_elements(ev)-1 DO BEGIN
                filename_ev_msec_list(i) = file_name_split(1) + ' ' + string(ev(i), format = '(F7.2)') + '  ' +  string(count_time(i), format='(F5.2)')
              ENDFOR
              stack_wb, filename
            ENDIF ELSE RETURN, filename
            IF NOT keyword_set(no_display) THEN Stack_process, filename, /no_align, /binary
        ENDIF ELSE BEGIN	; If  KeyWordSet no_save then send back currently selected image, of region, in channel
        	if n_elements(one_image) GT 0 then img_i = one_image + 1 else  img_i = 0
        	tmp = {t:'2d', x:findgen(n_elements(*Axes_list[2])), y: findgen(n_elements(*Axes_list[1])), d: Data_Array(*,*,img_i), $
        	       xl: ' ', yl: ' ', e: *Axes_list[0], dl:'stack'}
        	RETURN, tmp
        ENDELSE
        RETURN, filename
      END
    ELSE: BEGIN
      PRINT, "Don't know how to deal with ", N_ELEMENTS(data_dimensions), ' dimensions of data.'
      ;RETURN_STRUCTURE = {}
    END
  ENDCASE

ENDIF ELSE BEGIN
  PRINT, 'No data found in ', file, ' | ', region_path, '/', channel_name
  ;RETURN_STRUCTURE =  {}
ENDELSE
H5F_CLOSE, file_id
;print, "return flag", RETURN_STRUCTURE
RETURN, RETURN_STRUCTURE
END

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
FUNCTION read_nexus, file, filter=filter, channel=channel, region = region, silent = silent, $
         verbose=verbose, group=group, header_only=header_only, one_image = one_image, $
         no_display=no_display, no_save=no_save, xy_correct=xy_correct,_extra=e
on_error,2
@stack_process_com
@axis_com
@nexus_com

;COMMON volume_data, image_stack
;COMMON sdformat, starnum

;if float(!version.release) LT 8.0 then begin
;axis_log, 'NeXus data files cannot be read'
;axis_log, 'by IDL versions before IDL 8.0'
;return, 0
;endif

if keyword_set(no_save) then Nexus_no_save = 1 else Nexus_no_save = 0

if NOT keyword_set(group) then group = 0 else begin
  if n_elements(axis_ID) EQ 1 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

if n_elements(file) eq 0 then begin  ;popup file dialog box
   if not keyword_set(filter) then fltr='*.hdf5' 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 file and read metadata -------
; first force to be a header
t = findfile(file)
IF strlen(t(0)) EQ 0 then begin
     axis_log, 'read_nexus: Header file not found'
     return, s
ENDIF

nexus_split_name = ax_name(file)
nexus_lastfile = nexus_split_name(1)+'.'+nexus_split_name(2)     ; define for later use
nexus_path = nexus_split_name(0)


nexus_regions = get_nexus_regions(file)
nexus_region_info = get_nexus_region_info(file, nexus_regions(0))
;print, "nexus region info ", TAG_NAMES(nexus_region_info), nexus_region_info, nexus_region_info.Energy
StackAxis = create_struct('Points', nexus_region_info.Energy)
;print, "StackAxis ", StackAxis
ScanDef = create_struct('Type',nexus_region_info.Type,'Channel_labels',nexus_region_info.Channels,'n_regions',n_elements(nexus_regions),'StackAxis',StackAxis)


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

; NB the region and channel numbers are identified using a 'number_from_1' format but IDL wants 'number_from_0' format
if NOT keyword_set(region) then region=0 else region = region - 1		; (27-mar-16 aph - add 'NOT' and #-shifter)
if NOT keyword_set(channel) then channel=0 else channel = channel - 1

; if keyword_set(region)  then region=0			; B. Watts version of 14-Mar-16
; if keyword_set(channel) then channel=0

; if n_elements(xy_correct)EQ 0 then xy_correct = 0		; KLUGE !

print, 'file ',file
print, 'nexus_regions(region) = ', nexus_regions(region)
print, 'nexus_region_info.Channels(channel) ', nexus_region_info.Channels(channel)
print, 'xy_correct ', xy_correct_flag
; print, 'one_image ', one_image

DDD = read_NXdata(file, nexus_regions(region),nexus_region_info.Channels(channel),xy_correct=xy_correct_flag,one_image=one_image)
return, DDD

END
