% Function to read data and attributes (slope,intercept) from a satellite HDF file.
% Copyright by Jasmine S Nahorniak, Oregon State University.
% jun 2001 
% [data, slope, intercept, equation, name]=readMODISHDF(pathfile);

function [data, slope, intercept, equation, name]=readMODISHDF(pathfile);

%[file, path]=uigetfile('*.*','Open MODIS HDF File'); 

SD_id = hdfsd('start',pathfile,'read');   % open HDF file

%[ndatasets,nglobal_attr,status1] = hdfsd('fileinfo',SD_id); % general info about file


% NOTE: to check which reference number is applicable, loop through this next
% command for 40 or so reference numbers, letting "name" from the
% IloadHDFdata routine print to the screen.  The ref number is the
% order that the desired quantity appears.
%for i=1:40, 
%[data, slope, intercept, name]=IloadHDFdata(SD_id,i);
%i
%name
%end

[data, slope, intercept, equation, name]=IloadHDFdata(SD_id,2);  % MODIS product (chlorophyll, nLw_443, etc)
%[data, slope, intercept, equation, name]=IloadHDFdata(SD_id,11);  % MODIS l3 mapped sum product


status5 = hdfsd('end',SD_id);   % close HDF file

% if didn't find any scaling equations, look in main metadata
if equation==-1,
  file_id=hdfh('open',pathfile,'read',50);  % the number 50 is ignored 
  status = hdfv('start',file_id); 
  vdata_ref=hdfvs('find',file_id,'Scaling Equation');
  if vdata_ref==0,
     vdata_ref=hdfvs('find',file_id,'Scale_type');
  end
  if vdata_ref==0,
     vdata_ref=hdfvs('find',file_id,'scale_type');
  end
  if vdata_ref~=0,
     vdata_id=hdfvs('attach',file_id,vdata_ref,'r');
     [fields,count] = hdfvs('getfields',vdata_id);
     status=hdfvs('setfields',vdata_id,fields); 
     [vdata,junk]=hdfvs('read',vdata_id,1);
     equation=vdata{1}'; 
     status=hdfvs('detach',vdata_id);
  end
  % find slope and intercept if appropriate 
  if equation~=-1,
   vdata_ref=hdfvs('find',file_id,'Slope');
   if vdata_ref==-1,
      vdata_ref=hdfvs('find',file_id,'slope');
   end
   if vdata_ref~=-1,
      vdata_id=hdfvs('attach',file_id,vdata_ref,'r');
      [fields,count] = hdfvs('getfields',vdata_id); 
      status=hdfvs('setfields',vdata_id,fields);  
      [vdata,junk]=hdfvs('read',vdata_id,1);    
      slope=vdata{1}'; 
      status=hdfvs('detach',vdata_id);
   end
   vdata_ref=hdfvs('find',file_id,'Intercept'); 
   if vdata_ref==-1,
      vdata_ref=hdfvs('find',file_id,'intercept');
   end
   if vdata_ref~=-1,
      vdata_id=hdfvs('attach',file_id,vdata_ref,'r');
      [fields,count] = hdfvs('getfields',vdata_id); 
      status=hdfvs('setfields',vdata_id,fields);  
      [vdata,junk]=hdfvs('read',vdata_id,1);    
      intercept=vdata{1}'; 
      status=hdfvs('detach',vdata_id);
   end
  end
  status = hdfv('end',file_id);
  status=hdfh('close',file_id);
  equation='y= Slope * x + Intercept;';
end


% *************************************************************************
% **********
function [data, slope, intercept, equation, name]=IloadHDFdata(SD_id,ref);

sds_index = hdfsd('reftoindex',SD_id,ref);  % create an index number from the reference number

sds_id = hdfsd('select',SD_id,sds_index); % create an id number

[name,rank,dimsizes,data_type,nattrs,status2] = hdfsd('getinfo',sds_id); % info about data

% lots of variations below for compatability with old MODIS oceans files (first is current format 06/01)
slopeid=hdfsd('findattr',sds_id,'Slope');
if slopeid==-1,
  slopeid=hdfsd('findattr',sds_id,'slope');
end
interceptid=hdfsd('findattr',sds_id,'Intercept');
if interceptid==-1,
   interceptid=hdfsd('findattr',sds_id,'intercept');
end
equationid=hdfsd('findattr',sds_id,'Scale_type');
if equationid==-1,
  equationid=hdfsd('findattr',sds_id,'Scaling Equation');
end
if equationid==-1,
  equationid=hdfsd('findattr',sds_id,'scale_type');
end

[slope,slopestat] = hdfsd('readattr',sds_id,slopeid);
[intercept,interceptstat] = hdfsd('readattr',sds_id,interceptid);
[equation,equationstat] = hdfsd('readattr',sds_id,equationid);
if slopestat==-1,
   slope=1;
   intercept=0;
   equation=-1;
end

%fprintf(1,'%s\n',['Parameter: ' name]);

start=0*ones(1,rank);               % starting data point
stride=1*ones(1,rank);              % step between data points in each direction
edge=dimsizes;               % number of data points read in each direction
[data,status3] = hdfsd('readdata',sds_id,start,stride,edge);  % read dataset
%data=double(data);   % convert from int16 to double precision
status4 = hdfsd('endaccess',sds_id);  % close dataset


return