function [varargout]=mod06Reader(varToRead,area,fullFileName)

%FUNCTION  STRUCT = MOD06READER(VARTOREAD,AREA)
% NAME: 
%    MOD06READER
%
% PURPOSE:
%    Read the product information according to varToRead from a MODIS MOD06 HDF
%    file
%
%    For 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.
%
% CALLING SEQUENCE:
%    MOD06READER(VARTOREAD,AREA)
%
% INPUTS:
%    VARTOREAD (string) The variable to read from the MOD06 HDF file
%
% OPTIONAL INPUTS:
%    AREA - 4 element array containing information about a selected rectangular box
%    FULLFILENAME - if available, then don't open again - used in textOutput routine
%
% OUTPUTS:
%    Plots of figures showing the data
% OPTIONAL OUTPUTS:
%    OUTDATA    image of the data
%    fullFileName   used for textOutput routine to indicate to use the first selected
%                   mod06 files for all subsequent products selected
%
% RESTRICTIONS:
%    Only tested on Matlab version 5.3.0 (R11) - performance under
%    other versions of Matlab is unknown
%    Kevin Baggett - Tested on Matlab 6 - works fine
%
% EXAMPLES:
%    phaseproduct = mod06Reader('Cloud Phase',h.rect);
%
% Adapted from code by:
%    Shaima Nasiri           shaima.nasiri@ssec.wisc.edu
%    UW-Madison
%    CIMSS/SSEC
%    16 July, 2000
% UPDATED BY:
%    Kevin Baggett
%    UW-SSEC
%    Tue Feb 26 19:01:04 UTC 2002

switch varToRead
case 'Cloud Phase'
    varInFile='Cloud_Phase_Infrared';
case 'Cloud Top Temperature'
    varInFile='Cloud_Top_Temperature';
case 'Cloud Top Pressure'
    %varInFile='Cloud_Top_Pressure_Infrared';
    varInFile='Cloud_Top_Pressure';
case 'Effective Cloud Emissivity'
    varInFile='Cloud_Effective_Emissivity';
end

%try
    if nargin ~=3
        if isunix
            [fileName, pathName]=uigetfile('../Data/M*D06*.hdf',sprintf('%s File',varToRead));
        else
            [fileName, pathName]=uigetfile('..\Data\*',sprintf('%s File',varToRead));
        end
		
        fullFileName=[pathName fileName];
    end
%       catch
%           msgbox('File not loaded!');
%       end

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

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

% open HDF file for reading
 SD_id = hdfsd( 'start', fullFileName, 'read' );
 if (SD_id < 0)
    error(sprintf('HDF file %s was not opened.',fullFileName));
 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];
 % changes required if area argument is specified
 if (nargin >= 2 & ~isempty(area))
    % arrays read are only 406x270, so need to divide area arguments by five
    start(1) = min( [max([area(1) 0]) (npixels_across - 1)] );
    %start(1) = min( [max([area(1)/5 0]) (npixels_across - 1)] );
    start(2) = min( [max([area(2) 0]) (npixels_along - 1)] );
    %start(2) = min( [max([area(2)/5 0]) (npixels_along - 1)] );
    edge(1) = min( [max([area(3) 0]) (npixels_across - start(1))] );
    %edge(1) = min( [max([area(3)/5 0]) (npixels_across - start(1))] );
    edge(2) = min( [max([area(4) 0]) (npixels_along - start(2))] );
    %edge(2) = min( [max([area(4)/5 0]) (npixels_along - start(2))] );
 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) ;
 
 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

% adjust image array to work with matlab - 
% convert values to double precision and rotate and flip the image
pimage = double(scale) * (double(image_temp)  - double(offset));
pimage = flipud(rot90(pimage));

% if output argument specified, output data; otherwise, plot the data
if (nargout==1)
    varargout{1}=pimage;
elseif (nargout==2)
    varargout{1}=pimage;
    varargout{2}=fullFileName;
else
    mod06Plotter(pimage,varToRead,fileName);
end

  
