PRO HDF_SD_CONCATENATE, FILE_LIST, FILE_OUT

;+
; Concatenate all SDS arrays from a list of input files to an output file
;
; Usage: HDF_SD_CONCATENATE, FILE_LIST, FILE_OUT
;
; where
;
; FILE_LIST is an array of HDF4 files to be concatenated
; FILE_OUT is the name of the output file
;-

;- Check arguments
if (n_elements(file_list) eq 0) then message, 'Argument FILE_LIST is undefined'
if (n_elements(file_out) eq 0) then message, 'Argument FILE_OUT is undefined'

;------------------------------------
; GET SDS NAMES
;------------------------------------

;- Get list of variable names from first input file
hdfid = hdf_sd_start(file_list[0])
result = hdf_sd_varlist(hdfid)
hdf_sd_end, hdfid
if (result.nvars eq 0) then message, 'No SDS variables were found in first input file'
var_list = result.varnames

;- Open output file
hdfid_out = hdf_sd_start(file_out, /create)

;- Loop over variable names
for var_index = 0, n_elements(var_list) - 1 do begin

  ;- Get this SDS name
  var_name = var_list[var_index]
  
  ;------------------------------------
  ; COPY SDS DATA
  ;------------------------------------

  ;- Loop over input files
  for file_index = 0, n_elements(file_list) - 1 do begin

    ;- Get this input file name
    file_name = file_list[file_index]
    if (var_index eq 0) then print, file_name
        
    ;- Open this input file
    hdfid_in = hdf_sd_start(file_name)

    ;- Read this SDS
    hdf_sd_varread, hdfid_in, var_name, data

    ;- Get number of dimensions
    dims = size(data, /dimensions)
    ndims = n_elements(dims)
    
    ;- Concatenate the data
    if (file_index eq 0) then begin
      data_out = temporary(data)
    endif else begin
      case ndims of
        1 : data_out = [temporary(data_out), temporary(data)]
        2 : data_out = [[temporary(data_out)], [temporary(data)]]
        3 : data_out = [[[temporary(data_out)]], [[temporary(data)]]]
      endcase
    endelse

    ;- Close this input file
    hdf_sd_end, hdfid_in

  endfor

  ;- Write concatenated data to output file
  hdf_sd_varwrite, hdfid_out, var_name, temporary(data_out)

  ;------------------------------------
  ; COPY ATTRIBUTES
  ;------------------------------------
  
  ;- Open first input file
  hdfid = hdf_sd_start(file_list[0])

  ;- Get list of attributes for this variable from first input file
  result = hdf_sd_attlist(hdfid, var_name)
  att_list = result.attnames

  ;- Loop over attributes  
  for att_index = 0, n_elements(att_list) - 1 do begin

    ;- Get this attribute name
    att_name = att_list[att_index]

    ;- Copy attributes for this variable from first input file to output file  
    if (att_name ne '') then begin

      ;- Get the attribute from input file
      varid = hdf_sd_select(hdfid, hdf_sd_nametoindex(hdfid, var_name))
      attid = hdf_sd_attrfind(varid, att_name)
      hdf_sd_attrinfo, varid, attid, data=att_data
      hdf_sd_endaccess, varid

      ;- Write the attribute to output file
      varid = hdf_sd_select(hdfid_out, hdf_sd_nametoindex(hdfid_out, var_name))
      hdf_sd_attrset, varid, att_name, att_data
      hdf_sd_endaccess, varid

    endif
  
  endfor

  ;- Close first input file
  hdf_sd_end, hdfid
  
endfor

;- Close output file
hdf_sd_end, hdfid_out

END