; Copyright (c) 1998-2019 A.P. Hitchcock, B. Watts  All rights reserved
;+
;NAME:
; READ_NEXUS
;
;LAST CHANGED: -----------------------------------  06 Jan 2019 (aph)
;  PRIOR CHANGES:   01-Oct-18 (aph) from 11-Apr-18 (aph) from aph, 24-Jul-17; aph, from 19-Feb-17;
; (aph from 09-Oct-16 (aph), from 09-Sep-16 (bw)
;
;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: FILE   filename required
;
; KEYWORDS:

; NEXUS_FLIP flip images (needed for present co-ordinate syste definition of py-STXM image files)
; 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)
; 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 data
; SILENT    do not print any tracking information
; ONE_IMAGE if defined, index of image to be read (first image is number 1)
; XY_CORRECT - flag - if set, read interferometer (x,y) co-ordinates and interpolate data
; _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
; (25-Apr-16 aph) convert DOUBLES to FLOATs for all variables in stacks
; (08-May-16 aph) adapt to read multi-channel data (examples: Soleil: Sample_image_16-03-31_018.hdf5)
;                 defined norm_flag if not yet done (stand_alone)
; (25-Aug-16 bw ) make preview loading more robust
; (26-Aug-16 bw ) stop square pixels being applied to linescans
; (05-Sep-16 aph) use '$' continuation to reduce length of line 683 (pre-IDL7.0 limited to 255 characters)
;				  remember channel and region choices
; (09-Sep-16 bw ) ensure that recalled region and channel numbers are not greater than region and channel list lengths
; (09-Oct-16 aph) add FLIP option (needed for present co-ordinate system definition of py-STXM image files)
; (19-Feb-17 aph) correct y-axis in aborted image scans
; (24-Jul-17 aph) add KLUGE to correctly read-in single image files INCORRECTLY labelled as stacks by pySTXM
; (11-Apr-18 aph) fix read of point spectra; get single image working
; (01 Oct 18 aph) use h5_Close to release LUN after reading
; (16-Dec-18 aph) read part stacks  which are missing Tag name END_TIME; truncate aborted stacks to acquired images
; (29-Dec-18 aph) set default filename to the last 3 characters of theNeXus name (assumed to be a number)
; (06-Jan-19 aph) modify check for images with 0 content in a stack or stack map (from average = 0 to sum = 0)
;-

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_polarization, file, region_path
@nexus_com
polarization = 0 ; use this default if can't find a value in the file
file_id = H5F_OPEN(file)
group_list = list_nexus_groups(file_id, region_path)
FOR I=0, N_ELEMENTS(group_list)-1 DO BEGIN
  IF group_list[I].Class EQ 'NXcollection' THEN BEGIN
    collection_list = list_nexus_groups(file_id, region_path+'/'+group_list[I].Name)
    FOR J=0, N_ELEMENTS(collection_list)-1 DO BEGIN
      IF STRLOWCASE(collection_list[J].Name) EQ 'polarization' THEN BEGIN
        field_list = list_nexus_datasets(file_id, region_path+'/'+group_list[I].Name+'/'+collection_list[J].Name)
        FOR K=0, N_ELEMENTS(field_list)-1 DO BEGIN
          IF TOTAL(STRLOWCASE(field_list[K].Name) EQ ['value']) THEN BEGIN
            D_id = H5D_OPEN(file_id, region_path+'/'+group_list[I].Name+'/'+collection_list[J].Name+'/'+field_list[K].Name)
            polarization = float(H5D_READ(D_id))
            H5D_CLOSE, D_id
          ENDIF
        ENDFOR
      ENDIF
    ENDFOR
  ENDIF
ENDFOR
RETURN, polarization
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)
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)
Ch_field_list = list_nexus_datasets(file_id, region_path+'/'+channel_list[0])
FOR I=0, N_ELEMENTS(Ch_field_list)-1 DO BEGIN
  IF TOTAL(STRLOWCASE(Ch_field_list[I].Name) EQ ['count_time']) THEN BEGIN
    D_id = H5D_OPEN(file_id, region_path+'/'+channel_list[0]+'/'+Ch_field_list[I].Name)
    count_time = H5D_READ(D_id)
    H5D_CLOSE, D_id
    Region_info = CREATE_STRUCT(Region_info,'count_time',count_time)
  ENDIF
  IF TOTAL(STRLOWCASE(Ch_field_list[I].Name) EQ ['data']) THEN BEGIN
    D_id = H5D_OPEN(file_id, region_path+'/'+channel_list[0]+'/'+Ch_field_list[I].Name)
    data = H5D_READ(D_id)
    H5D_CLOSE, D_id
    Region_info = CREATE_STRUCT(Region_info,'data_size',N_ELEMENTS(data))
  ENDIF
ENDFOR
field_list = list_nexus_datasets(file_id, region_path)
FOR I=0, N_ELEMENTS(field_list)-1 DO BEGIN
  IF TOTAL(STRLOWCASE(field_list[I].Name) EQ ['definition']) AND TOTAL(STRLOWCASE(field_list[I].String) EQ ['nxstxm']) THEN BEGIN
    attr_list = get_nexus_attributes(file_id, region_path+'/'+field_list[I].Name)
    IF TOTAL(STRLOWCASE(TAG_NAMES(attr_list)) EQ ['version']) THEN BEGIN
      Region_info = CREATE_STRUCT(Region_info,'NXstxm_version',attr_list.version)
    ENDIF
  ENDIF
  IF TOTAL(STRLOWCASE(field_list[I].Name) EQ ['start_time']) THEN BEGIN
    D_id = H5D_OPEN(file_id, region_path+'/'+field_list[I].Name)
    start_time = H5D_READ(D_id)
    H5D_CLOSE, D_id
    Region_info = CREATE_STRUCT(Region_info,'start_time',start_time)
  ENDIF
  IF TOTAL(STRLOWCASE(field_list[I].Name) EQ ['end_time']) THEN BEGIN
    D_id = H5D_OPEN(file_id, region_path+'/'+field_list[I].Name)
    end_time = H5D_READ(D_id)
    H5D_CLOSE, D_id
    Region_info = CREATE_STRUCT(Region_info,'end_time',end_time)
  ENDIF
ENDFOR
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
OLD_LINESCAN_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)
; print, "Region_info = ", TAG_NAMES(Region_info)

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
        AXIS_FLAG = 0
        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
            AXIS_FLAG = 1
            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[N_ELEMENTS(Axes_list)-1] = Dim_Array ; set pointer to data array
            H5D_CLOSE, data_id
            BREAK
          ENDIF
        ENDFOR
        IF NOT AXIS_FLAG THEN BEGIN; This is needed for the very old style (where dims [E,Y,X] was forced) linescans where the axis attributes are: E=1,Y=3,X=3
          OLD_LINESCAN_FLAG = 1
        ENDIF
      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
  IF TOTAL(STRLOWCAsE(TAG_NAMES(Region_info)) EQ 'nxstxm_version') THEN BEGIN
    IF Region_info.NXstxm_version EQ '0.1' THEN BEGIN
      print, "NXstxm v0.1: Remove extra dimensions."
      Data_Array = REFORM(Data_Array,/OVERWRITE)
      FOR I=0, N_ELEMENTS(Axes_list)-1 DO BEGIN
        IF N_ELEMENTS(*Axes_list[I]) GT 1 THEN BEGIN
          PUSHLIST, temp_Axes_list, Ptr_New([0], /No_Copy) ; add pointer to list
          *temp_Axes_list[N_ELEMENTS(temp_Axes_list)-1] = *Axes_list[I] ; set pointer to data array
        ENDIF
      ENDFOR
      Axes_list = temp_Axes_list
    ENDIF
  ENDIF

; =============  KLUGE, KLUGE, KLUGE ========= 24-Jul-17 (aph) ================
;  pySTXM frequently sets nexus_head.Type to 'sample image stack' when only one image is acquired
;  here I make changes to make it read as a single image
  print, 'Data array dimensions ', SIZE(Data_Array)
  IF stxm_scan_type EQ 'sample image stack' AND n_elements(energy) EQ 1 then begin
       	print, 'CLS pySTXM special readin - forcing 1-image stack to read as a normal sample image'
        t = ax_name(file)
        print, 'Number of axes in ', t(1), '  is ', SIZE(Axes_list, /N_ELEMENTS),  '  reducing to 2 '
	    energy = energy[0]
	    count_time = count_time[0]
	    Data_Array = Data_Array[*,*,0]
	    Axes_list = Axes_list[1:*]
	    stxm_scan_type = 'single image (*)'
  ENDIF
  ; =============  END  KLUGE, KLUGE, KLUGE ========= 24-Jul-17 (aph) ================

  print, "Reading ", STRING(stxm_scan_type)+" with "+STRING(format='(i1)',SIZE(Data_Array, /N_DIMENSIONS))+ " dimensions."
  print, " Number of axes is ", SIZE(Axes_list, /N_ELEMENTS)
  IF keyword_set(one_image) AND SIZE(Data_Array, /N_DIMENSIONS) GT 2 AND N_ELEMENTS(energy) GT 1 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
  nexus_Dwell = count_time[0]
  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, dn: 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)'
          	RETURN_STRUCTURE = {t:'2d', x: ax_dts(*Axes_list[0]), y:ax_dts(*Axes_list[1]), $
               	d:TRANSPOSE(ax_dts(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)'
            if keyword_Set(nexus_flip) then begin							; ******************
            	t_img = rotate(ax_dts(Data_Array),7)				; implementation of FLIP for CLS pySTXM images
            	print, 'image flipped'								; FLIP parameter passed through nexus_com
            endif else t_img = ax_dts(Data_Array)
; ================================
; APH added: 19-Feb-2017
;  Test if Y-axis stalls (aborted image); if so replace with continuation of non-stalled y-axis points
; NB for CLS c-STXM the Y-axis in some circumstances is interferometer values which
;    are not accurate, so I am using the average of the y-axis step before the stall to define y-step
		  	y_now = ax_dts(*Axes_list[0])
;		  	print, '# of y-points ', ny
		  	if ny GT 2 then begin					; traps out errors when only 1 y-value
		  	 	if (y_now(ny-1) - y_now(ny-2)) EQ 0 then begin
			  		axis_log, "Read_NEXUS: aborted image scan - correcting y-axis"
			  		t = ny
			  		while (y_now(t-1) - y_now(t-2)) EQ 0 do begin
	;		  			print, t-1, '  ', y_now(t-1)
			  			t = t - 1
			  		endwhile
			  		yd = fltarr(t)
			  		for i = 0, t-2 do yd(i) = y_now(i+1) - y_now(i)
			  		del_y = moment(yd, sdev=y_sdev)
			  		axis_log, 'average step ' + string(del_y(0)) + ' +/- ' + string(y_sdev(0))
					for i = t-1, ny-1 do begin
						y_now(i)=y_now(i-1)+del_y(0)
					endfor
				endif
		  	endif
; ================================

            tmp =  {t:'2d', x:ax_dts(*Axes_list[1]), y:y_now, d:t_img, e:float(energy(0)), xl:X_string, yl:Y_string, dl:Data_File_Name+'   '+channel_name }
            IF nexus_head.Type NE 'osa focus' AND  nexus_head.Type NE 'sample focus' THEN BEGIN
            	if n_elements(tmp.y) GE 2 then begin
            		tmp = ax_interp_square(tmp)
            	endif else begin
            		axis_log, 'Only 1 line in data file'
				endelse
              	RETURN_STRUCTURE = tmp
            ENDIF ELSE BEGIN
              RETURN_STRUCTURE =  tmp
            ENDELSE
        ENDELSE
      END
    3:BEGIN
;        print, '(3D scan)'
        nx = fix(n_elements(*Axes_list[2]))
        ny = fix(n_elements(*Axes_list[1]))
        n_e = fix(n_elements(*Axes_list[0]))
;        help, nx, ny
        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
  ; -------- set default file name to just the number, assumed to be last 3 characters before '.'  (aph, 29-dec-2018)
  		  t = ax_name(Data_File_Name)
		  filename = strmid(t(1), 2,3,/reverse_offset)
		  filename = pickfile2(filter='*.ncb', file=filename, title='binary stack file ', WRITE=write)
          IF strlen(filename) GT 0 THEN BEGIN
              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))		; (24-Apr-16 aph)  syntax of filename_ev_msec_list not consistent with stack_process
              ; 	    filename_ev_msec_list(i) = filenm  + ' : '+$
	          ;			strtrim(string(ev(i),format='(f10.2)'),2)+' eV '+$
	          ;			strtrim(string(dwell,format='(f7.2)'),2)+' msec'

              FOR i = 0, n_elements(ev)-1 DO BEGIN
                filename_ev_msec_list(i) = file_name_split(1) + ' ' + string(ev(i), format = '(F10.2)') + $
                                           ' eV ' +  string(count_time(i), format='(F7.2)') + ' msec'
              ENDFOR
              image_stack = ax_dts(image_stack)	; convert doubles to single floats
              if keyword_Set(nexus_flip) then begin							; ******************
            	for i = 0, n_elements(ev)-1 do begin
            		image_stack(*,*,i) = rotate(image_stack(*,*,i),7)				; implementation of FLIP
            	endfor
            	print, 'stack flipped'								; FLIP parameter passed through nexus_com
              endif

;  ------- check for null images due to aborted stack and remove them ------------ ; (added by APH  16-Dec-18, changed 6-Jan-18)
				FOR i = 0, n_elements(ev)-1 DO BEGIN
					test = total(image_stack(*,*,i),/double, /nan)
					IF test EQ 0 then begin
						axis_log, 'ABORTED STACK'
						axis_log, 'last image is ' + strtrim(string(i),2) + ' of ' + strtrim(string(n_elements(ev)-1),2)
						tmp_stack = image_stack(*,*,0:i-1)
						image_stack = tmp_stack
						ev = ev(0:i-1)
						n_e = i-1
						goto, all_good
					ENDIF
				ENDFOR
				all_good:
              	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:ax_dts(dindgen(n_elements(*Axes_list[2]))), y: ax_dts(dindgen(n_elements(*Axes_list[1]))), d: ax_dts(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,flip=flip, _extra=e
on_error,2
@stack_process_com
@axis_com
@nexus_com

;COMMON volume_data, image_stack
;COMMON sdformat, starnum


if n_elements(norm_flag) LE 0 then norm_flag=0  ; so can use stand alone

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 = ax_name(file)
if t(2) NE 'hdf5' then file = t(0) + t(1) +'.hdf5'
t = findfile(file)
IF strlen(t(0)) EQ 0 then begin
     if NOT keyword_set(silent) then axis_log, 'read_nexus: ' + file + '   Header file not found'
     h5_close	; (01-Oct-18 aph) to release file after opening & reading
     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))
polarization = get_polarization(file, nexus_regions(0))
StackAxis = create_struct('Points', nexus_region_info.Energy)
;print, "StackAxis ", StackAxis
;
; ------- kluge to avoid error "Tag name END_TIME is undefined for structure <Anonymous>"   16-Dec-18
test = tag_names(nexus_region_info)
for i = 0, n_elements(test)-1 do begin
	if strlowcase(test(i)) EQ 'end_time' then goto, continue
endfor
end_time = nexus_region_info.start_time  ; if do not find an end time,  set it to the start time
nexus_region_info = CREATE_STRUCT(nexus_region_info, 'end_time', end_time)

continue:

ScanDef = create_struct('Type',nexus_region_info.Type,'Channel_labels',nexus_region_info.Channels, 'n_regions', $
          n_elements(nexus_regions),'StackAxis',StackAxis,'count_time',nexus_region_info.Count_time,'data_size', $
           nexus_region_info.Data_size,'start_time',nexus_region_info.start_time,'end_time',nexus_region_info.end_time,'polarization',polarization )

if keyword_set(header_only) then begin
; ------- print header information
; XDisplayFile,file_hdr
;PRINT, ScanDef
  h5_close	; (01-Oct-18 aph) to release file after opening & reading
  return, ScanDef
endif

; NB the region and channel numbers are identified using a 'number_from_1' format but IDL wants 'number_from_0' format
; (aph 5-Sep-16)  add remembering last selected region and channel
if NOT keyword_set(region) then begin
	if n_elements(nexus_region) GT 0 then begin
		region = nexus_region
	endif else region=0
endif else region = region - 1		; (27-mar-16 aph - add 'NOT' and #-shifter)
nexus_region = region

if NOT keyword_set(channel) then begin
	if n_elements(nexus_ch) GT 0 then begin
		channel = nexus_ch
	endif else channel=0
endif else channel = channel - 1
nexus_ch = channel

; if n_elements(xy_correct)EQ 0 then xy_correct = 0		; KLUGE !
 print, ' '
 print, 'file ',file
 print, 'name of nexus_regions = ', nexus_regions(MIN([N_ELEMENTS(nexus_regions)-1,region]))
 print, 'name of data channels(channel) ', nexus_region_info.Channels(MIN([N_ELEMENTS(nexus_region_info.Channels)-1,channel]))
; print, 'xy_correct ', xy_correct_flag
; print, 'one_image ', one_image
 print, ' '

DDD = read_NXdata(file, nexus_regions(MIN([N_ELEMENTS(nexus_regions)-1,region])),nexus_region_info.Channels(MIN([N_ELEMENTS(nexus_region_info.Channels)-1,channel])),xy_correct=xy_correct_flag,one_image=one_image)
h5_close	; (01-Oct-18 aph) to release file after opening & reading
return, DDD

END
