pro read_hdf4, file_name, sds_num, TOASCII=toascii

;- Script to read in an hdf4 file into an array or output to an ascii
;  file.
;- sds_num = 1 for density, 2 = internal energy, 3-5 = x-z velocity.

;- Read the data set into a 3D array
 
DFSD_SETINFO, /RESTART
DFSD_GETINFO, file_name, NSDS=NumSDS

IF (sds_num gt numSDS) THEN BEGIN
   print, ' read_hdf: Cannot read sds#', sds_num
   print, ' there are only ', numsds,' scientific data sets in ', file_name
   print, ' returning zero valued fltarr'
   data = fltarr(3,3,3)
   RETURN
ENDIF
  
FOR i=1,sds_num DO BEGIN
   r=0.
   l=0.
   HDF_DFSD_GETINFO, file_name, label = l, dims=d, type=t,range=r
ENDFOR
  
print, ' read ',l, ' data from ', file_name
hdf_DFSD_GETDATA, file_name, data

data_dim = size(data)

print, 'data set info: ',data_dim

min1 = MIN(data)
max1 = MAX(data)

print, 'min, max =', MIN(data), MAX(data)


IF KEYWORD_SET(TOASCII) THEN BEGIN

   openw, unit, "ascii_data.dat", /get_lun

   if (data_dim(0) EQ 1) then begin
      for i=0, data_dim(1)-1 do begin
         printf, unit, i, data(i)
      endfor
   endif

   if (data_dim(0) EQ 2) then begin
      for i=0, data_dim(1)-1 do begin
         for j=0, data_dim(2)-1 do begin
            printf, unit, FORMAT='(1(e16.8))',data(i,j)
         endfor
      endfor
   endif

   if (data_dim(0) EQ 3) then begin
      for i=0, data_dim(1)-1 do begin
         for j=0, data_dim(2)-1 do begin
            for k=0, data_dim(3)-1 do begin
               printf, unit, i, j, k, data(i,j,k)
            endfor
         endfor
      endfor
   endif

   free_lun, unit
   
ENDIF




stop
END