%---------------------------------------------------------------------
function [ds,varargout] = obtain_datafield(...
                                hdfID, dataName, doScaling, doNaN)
%---------------------------------------------------------------------
%-- 120219: Dont use HDF-EOS interface, it needs a seperate HDF view to
%   read attributes that is local to data field. Use directly the HDF
%   interface is simple and enough for my current application.
%   4 ways to read a HDFEOS file
%     1-- hdftool 
%     2-- hdfread !the high-level Matlab function. no way to read attributes.
%     3-- hdfeos low-level functions !need seperate HDF view to read attributres.
%     4-- hdf low-level functions !may be complicated for some HDFEOS files.
%--2022: the old HDF functions does not work. The use the matlab.io.hdf4.* package
      import matlab.io.hdf4.*

      dsIdx = sd.nameToIndex(  hdfID, dataName ); %pisoIdx=9, pvolIdx=19, pgeoIdx=29
      dsHdl = sd.select( hdfID, dsIdx );
      [dsname, dsdims, dstype, dsatts] =...
                                            sd.getInfo(dsHdl );
      sts = isempty(dsdims);                                      
      if sts~=0 disp(['---obtain_datafield.m: getinfo() sts~=0! ','dataName=',dataName]); keyboard; end


      %dsndims  = 2;
      %dsdims   = [NY,NX]; %hdf follows C, starts from 0 and is row first.
      dsndims = length(dsdims);
      dsstart  = zeros( 1,dsndims ); % Creates the vector [0 0]
      dsstride = []; 
      dsedges  = dsdims;
      [ds] = sd.readData(  dsHdl, dsstart, dsedges );
      sts = isempty(ds);
      if sts~=0 disp(['---obtain_datafield.m: read data sts~=0! ','dataName=',dataName]); keyboard; end

      attHdl       = sd.findAttr(  dsHdl, '_FillValue' );
      [misval] = sd.readAttr(  dsHdl, attHdl );
      sts = isempty(misval);
      if sts~=0 disp(['---obtain_datafield.m: read miss sts~=0! ','dataName=',dataName]); keyboard; end

      mis = ds==misval;

      %-- do scaling, if necessary. Only act on non-missing vlaues.
      %   pay attention to the usage of the sclfac and offset.
      %
      if strcmpi(doScaling,'doScl')

         attHdl       = sd.findAttr(  dsHdl, 'scale_factor' );
         [sclfac] = sd.readAttr(  dsHdl, attHdl );
         sts = isempty(sclfac);
         if sts~=0 disp(['---obtain_datafield.m: read scale sts~=0! ','dataName=',dataName]); keyboard; end
      
         attHdl       = sd.findAttr(  dsHdl, 'add_offset' );
         [offset] = sd.readAttr(  dsHdl, attHdl );
         sts = isempty(offset);
         if sts~=0 disp(['---obtain_datafield.m: read offset sts~=0! ','dataName=',dataName]); keyboard; end

         ds = single( ds );
         %ds = double(ds);
         ds(~mis) = sclfac*( ds(~mis) - offset );
         %ds = sclfac*ds+offset;
      end %if strcmpi(doScaling,'scl')


      %-- Matlab integer type CAN NOT be set to NaN! would be set to 0.
      %   but ds==misval is OK, although misval is a int16 integer.
      if strcmpi(doNaN,'doNaN')
         ds(mis) =NaN;
      end


      sd.endAccess(  dsHdl ); %need for WRITTING.
      
      varargout{1} = dsdims;

   %end %--for ip = 1:length(paraNames)

%keyboard
   
%{
   %Open swath
   swathID = hdfsw('attach', hdfID, swathName);

   %Reading Data from a Data Field
   %dataName='Retrieved_Moisture_Profile';

   [data, sts] = hdfsw('readfield', swathID, dataName, [], [], []);
   
   %Reading lat and lon data
   %[lon, status] = hdfsw('readfield', swathID, 'Longitude', [], [], []);
   %[lat, status] = hdfsw('readfield', swathID, 'Latitude', [], [], []);

   %Detaching from the Swath Object
   hdfsw('detach', swathID);

   
   %Convert the data to double type for plot
   %data=double(data);
   %lon=double(lon);
   %lat=double(lat);

   %Reading attributes from the data field
   %dataName='Retrieved_Moisture_Profile';

   sds_index = hdfsd('nametoindex', SD_id, dataName);
   sds_id = hdfsd('select',SD_id, sds_index);


   %Reading filledValue from the data field
   fillvalue_index = hdfsd('findattr', sds_id, '_FillValue');
   [fillvalue, status] = hdfsd('readattr',sds_id, fillvalue_index);

   %Reading units from the data field
   units_index = hdfsd('findattr', sds_id, 'units');
   [units, status] = hdfsd('readattr',sds_id, units_index);

   %Reading scale_factor from the data field
   scale_index = hdfsd('findattr', sds_id, 'scale_factor');
   [scale, status] = hdfsd('readattr',sds_id, scale_index);

   %Reading add_offset from the data field
   offset_index = hdfsd('findattr', sds_id, 'add_offset');
   [offset, status] = hdfsd('readattr',sds_id, offset_index);

   %Terminate access to the corresponding data set
   hdfsd('endaccess', sds_id);
%}

end

