function [mod06data]=modis_mod06_read(mod06file, varToRead, area)

%FUNCTION [MOD06_STRUCT]=MODIS_MOD06_READ(MOD06FILE, KEYWORDS, AREA)
% NAME: 
%    MODIS_MOD06_READ
%
% PURPOSE:
%    Read the product information according to varToRead from a MODIS MOD06 HDF
%    file
%
%    For IR Cloud Phase, the output image has 8 possible values -
%        0 = clear  (from cloud mask)
%        1 = water
%        2 = ice
%        3 = mixed phase
%        4 = ice
%        5 = water
%        6 = uncertain
%      127 = missing data
%
%
%    This procedure uses only HDF calls (it does *not* use HDF-EOS),
%    and only reads from SDS and Vdata arrays 
%    (it does *not* read ECS metadata).
%
% CATEGORY:
%    MODIS utilities.
%
% INPUTS:
%    MOD06FILE (string) The full-path name of the MOD06 HDF file
%    VARTOREAD (string) The variable to read from the MOD06 HDF file
%              'IRphase', 'Cloud Top Temperature', 'Cloud Top Pressure', 
%              'Effective Cloud Emissivity', 'Effective Radius', 
%              'Optical Thickness', 'Retrieval Phase'
%
% OPTIONAL INPUTS:
%    AREA - 4 element array containing information about a selected rectangular box
%
% OUTPUTS:
%    MOD06DATA   (struct) mod06data.pimage  (data)
%                         mod06data.pimage_string (name of data variable)
%
% EXAMPLES:
%    phaseproduct = mod06Reader('MOD06_L2.A2002345.1915.hdf','Cloud Phase');
%
%    Shaima Nasiri           shaima.nasiri@ssec.wisc.edu
%    UW-Madison
%    CIMSS/SSEC
%    05 June, 2003
%


varToRead = lower(varToRead);
switch varToRead
  case {'cloud_phase','cloud_phase_infrared', 'irphase', 'cphsae'}
      varInFile='Cloud_Phase_Infrared';
      areadims = '5by5';
      
  case {'cloud_top_temperature', 'ctt', 'cloud_top_temperature'}
      varInFile='Cloud_Top_Temperature';
      areadims = '5by5';
      
  case {'cloud_top_pressure', 'ctp', 'pressure', ...
	'cloud_top_pressure_infrared'}
      varInFile='Cloud_Top_Pressure_Infrared';
      areadims = '5by5';
      
  case 'effective_cloud_emissivity'
      varInFile='Cloud_Effective_Emissivity';
      areadims = '5by5';
      
  case {'re', 'effective_particle_radius', 'effective radius'}
      varInFile = 'Effective_Particle_Radius';
      areadims = '1by1';
      
  case {'tau', 'cloud_optical_thickness', 'optical thickness'}
      varInFile = 'Cloud_Optical_Thickness';
      areadims = '1by1';
      
  case {'wp', 'water_path'}
      varInFile = 'Water_Path';
      areadims = '1by1';
      
  case {'cr', 'cirrus_reflectance'}
      varInFile = 'Cirrus_Reflectance';
      areadims = '1by1';
      
  case {'retphase', 'retrieval_phase'}
      varInFile = 'Quality_Assurance_1km';
      areadims = '1by1';
end

% Check inputs
 nargs = nargin; 
 if (nargs < 1)
    error(['Usage:  MODIS_MOD06_READ(MASKFILENAME)']);
 end
 if (nargs == 3)
    if (~isnumeric(area) | (length(area) ~= 4))
       error(['AREA argument  must be 4 element vector of the form' ...
	      ' [XO YO NX NY]']);
    end
 end



% Check for valid MOD06 HDF file
 if (~exist(mod06file,'file'))
    error(sprintf('Filename : %s was not found',fullFilename));
 end

 len_filename = length(mod06file);
 
 if (~strcmp(mod06file(len_filename-3:len_filename), '.hdf'))
    error(sprintf('Filename: %s is not an HDF file',mod06file));
 end

% open HDF file for reading
 SD_id = hdfsd( 'start', mod06file, 'read' );
 if (SD_id < 0)
    error(sprintf('HDF file %s was not opened.',mod06file));
 end

% data wanted is specified in varInFile string
 file_ex = hdfsd('nametoindex', SD_id, varInFile);
 file_id = hdfsd('select',SD_id, file_ex);
 [name,rank,dimsizes,data_type,nattrs,status] = hdfsd('getinfo',file_id  );
%npixels_across = dimsizes(1);
%npixels_along =  dimsizes(2);
%start = [0 ; 0];
% count = [1 ; 1];
% edge = [npixels_across ; npixels_along];
 ndims = length(dimsizes);
 start = zeros(ndims,1);
 count = ones(ndims,1);
 edge = dimsizes;
 % changes required if area argument is specified
 if (nargin >= 3 & ~isempty(area))
    if strcmp(areadims,'5by5')
    % arrays read are only 406x270, so need to divide area arguments by five
       start(1) = min( [max([area(1)/5 0]) (npixels_across - 1)] );
       start(2) = min( [max([area(2)/5 0]) (npixels_along - 1)] );
       edge(1) = min( [max([area(3)/5 0]) (npixels_across - start(1))] );
       edge(2) = min( [max([area(4)/5 0]) (npixels_along - start(2))] );
    else  %1 by 1 values
       start(1) = min( [max([area(1) 0]) (npixels_across - 1)] );
       start(2) = min( [max([area(2) 0]) (npixels_along - 1)] );
       edge(1) = min( [max([area(3) 0]) (npixels_across - start(1))] );
       edge(2) = min( [max([area(4) 0]) (npixels_along - start(2))] );
    end    
 end
 
 [data,status] = hdfsd('readdata',file_id,start,count,edge);
 image_temp  = data;

 attr_ex = hdfsd('findattr', file_id, 'scale_factor');
 [scale, status] = hdfsd('readattr', file_id,attr_ex) ;
% following line was originally used for testing
% disp(['scale factor is:  ' num2str(scale)]);    
 
 attr_ex = hdfsd('findattr', file_id, 'add_offset');
 [offset, status] = hdfsd('readattr', file_id,attr_ex) ;

% stop accessing data 
 status = hdfsd('endaccess',file_id);

% close HDF file
 status = hdfsd('end', SD_id);
 if (status < 0)
    warning(['HDF file' filename ' was not closed.']);
 end

 if strcmp(varToRead,'retphase')
%   byte3 =     flipud(rot90( squeeze(image_temp(3,:,:)) ));
    byte3 =     squeeze(image_temp(3,:,:));
 % have to add 256 because matlab bit functions assume positive
 % integers
    byte3 = double(byte3) + 256; 
    pimage = bitand(byte3, 7);
    
 else
       pimage = double(scale) * (double(image_temp)  - double(offset));
%       pimage = flipud(rot90(pimage));
 end
 

mod06data.pimage = pimage;
mod06data.pimage_string = varInFile;
