function phasedata = modis_phase_read(phasefilename,area)

%FUNCTION  PHASE_STRUCT = MODIS_PHASE_READ(PHASEFILENAME)
% NAME: 
%    MODIS_PHASE_READ
%
% PURPOSE:
%    Read the phase product information from a MODIS MOD06 HDF
%    product file
%
%    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:
%    MODIS_PHASE_READ(PHASEFILENAME)
%
% INPUTS:
%    PHASEFILENAME (string)  Name of MODIS MOD06 HDF file
%
% OPTIONAL INPUTS:
%    AREA - 4 element array containing information about a selected rectangular box
%
% OUTPUTS:
%    PHASEDATA    A struct with the format: 
%             .PIMAGE     a two dimensional array of phase product
%                         data
%             .LONGNAME   a string describing the date type
%             .PNAME      the HDF filename (e.g. 
%                         'MOD06CT_L2.A2000182.0435.hdf'
%
% OPTIONAL OUTPUTS:
%     none
%
% 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:
%    phasepath = '/modisnfs2/MODIS/2000/182shaima/';
%    phaseproduct = modis_phase_read( ...
%                  [phasepath 'MOD06CT_L2.A2000182.0435.hdf'] );
%
% 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
%    Now includes an area argument to select a region


% Check for valid MOD06 HDF file
 if (~exist(phasefilename,'file'))
    error(['Filename : ' filename ' was not found']);
 end

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

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

% data we want is 'Cloud_Phase_Infrared'
 phase_ex = hdfsd('nametoindex', SD_id, 'Cloud_Phase_Infrared');
 phase_id = hdfsd('select',SD_id, phase_ex);
 [name,rank,dimsizes,data_type,nattrs,status] = hdfsd('getinfo',phase_id  );
 npixels_across = dimsizes(1);
 npixels_along =  dimsizes(2);
 start = [0 ; 0];
 count = [1 ; 1];
 edge = [npixels_across ; npixels_along];
 % Kevin's changes required by addition of area argument
 if (nargin == 2)
    % phase array is 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))] );
 end
 
 [data,status] = hdfsd('readdata',phase_id,start,count,edge);
 image_temp  = data;

 attr_ex = hdfsd('findattr', phase_id, 'scale_factor');
 [scale, status] = hdfsd('readattr', phase_id,attr_ex) ;
 
 attr_ex = hdfsd('findattr', phase_id, 'long_name');
 [longname, status] = hdfsd('readattr', phase_id,attr_ex) ;

 attr_ex = hdfsd('findattr', phase_id, 'Cell_Along_Swath_Sampling');
 [sampling, status] = hdfsd('readattr', phase_id,attr_ex) ;

 attr_ex = hdfsd('findattr', phase_id, 'add_offset');
 [offset, status] = hdfsd('readattr', phase_id,attr_ex) ;

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


% close HDF file
 status = hdfsd('end', SD_id);
 if (status < 0)
    warning(['HDF file' phasefilename ' 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));

  phase_s = struct('pimage',[]);
  phase_s.pimage = pimage;
  phase_s.longname = longname;
  phase_s.pname = phasefilename;

% return phase data struct
  phasedata = phase_s;

