function readtrmm(filename)
%READTRMM read TRMM HDF file
%   READTRMM(FILENAME) reads the Scientific Data sets and V Data sets from the
%   TRMM HDF file FILENAME.  SDS and VDATA sets may be written to output files
%   based on user input.  The output file(s) are named FILENAME.<VDATA name> or
%   FILENAME.<SDS name>.
%
%   Example:
%      readtrmm('3A25.971201.5.HDF')
%
%   Notes:
%  
%   The original program used as the basis for this code was developed by the
%   NASA/Goddard Earth Sciences Distributed Active Archive Center.
%
%   This M-file is not an official component of MATLAB, has not been officially
%   tested, and so it is being provided "as is".  The user should be aware that
%   there is a minimal amount of error checking done in this function.
%

%  
%   Copyright 1984-2000 The MathWorks, Inc

%Open TRMM HDF file
file_id = hdfh('open',filename,'read',0);

% Open HDF SDS and VDATA Interfaces
sd_id = hdfsd('start',filename,'read');

status = hdfv('start',file_id);

[n_datasets,n_file_attrs, status] = hdfsd('fileinfo',sd_id);

disp(['Number of SDS arrays in file ' num2str(n_datasets) '.']);

for index = 0:(n_datasets-1)
  sds_id = hdfsd('select',sd_id, index);
  [sds_name, rank, dim_sizes, num_type, attributes, status] = hdfsd('getinfo', ...
						  sds_id);
  str = ['   %2d) %-15s dimensions = ' repmat('%7d ',1,rank) '\n'];
  disp(sprintf(str, index+1, sds_name,dim_sizes))
  status = hdfsd('endaccess',sds_id);
end

nparam = input('Number of parameters to write out or 0 to exit program. ');

if nparam == 0
   hdfml('closeall');
   return;
elseif nparam == n_datasets
  nindex = 0;
  for index = 1:n_datasets
    sds_index(index) = index;
  end
else
  temparm = input('Enter parameter numbers from above list as a MATLAB array. ');
  sds_index = temparm;
  nindex = length(temparm);
end

if nparam > 0 
  for index = 1:nindex
    num_element = 1;
    sds_idx = sds_index(index)-1;
    sds_id = hdfsd('select',sd_id,sds_idx);
    [sds_name, rank, dim_sizes, num_type, attributes, status] = hdfsd('getinfo',  sds_id);
    [gain,gain_err,offset,offset_err,cal_data_type,status] = ...
	hdfsd('getcal',sds_id);
    
    %Stopped at line 262 of read_trmm.c
    
    fprintf('\n ****************************************\n');
    fprintf('  SDS name  = %s\n',  sds_name);
    fprintf('  SDS type  = %4s\n', num_type);
    fprintf('  SDS rank  = %4d\n', rank);
    fprintf('  SDS scale = %f\n',  gain);
    fprintf('  SDS dims  = ');

    for j=1:rank
      fprintf('%6d',dim_sizes(j));
      num_element = num_element*dim_sizes(j);
    end
    fprintf('\n');
    
    stride = [];
    edges = [];
    start = [];
    for j = 1:rank
      edges(j) = dim_sizes(j);
      stride(j) = 1;
      start(j) = 0;
    end
    
    % Read TRMM sds arrays
    [data, status] =  hdfsd('readdata',sds_id,start,[],edges);
    
    %output sds to binary data file
    outfile = [filename '.' sds_name];
    sds_fp = fopen(outfile,'wb','b');
    if sds_fp == -1
      error('Failed to open output binary file.');
    end
    count = fwrite(sds_fp,double(data),num_type);
    fclose(sds_fp);
    status = hdfsd('endaccess',sds_id);
  end
end

status = hdfsd('end',sd_id);

% Now read in the Vdata tables and write to seperate output files.

vdata_ref = - 1;
vdata_ref = hdfvs('getid',file_id, vdata_ref);
if ( vdata_ref == -1 ) 
  error('\n No Vdatas found in data set\n');
end

fprintf('\n\n  List of Vdata name:               (rec_num   rec_size)\n\n');

total_vidx = 0;
while(vdata_ref ~= -1)
  vdata_id = hdfvs('attach',file_id, vdata_ref, 'r');

  [n_records, interlace, fields, vdata_size, vdata_name, status] = ...
      hdfvs('inquire', vdata_id);
  
  total_vref(total_vidx+1) = vdata_ref;
  total_vidx = total_vidx + 1;
  fprintf(' %2d)  %-30s   %5d     %5d \n', total_vidx, vdata_name, n_records, vdata_size);
  status = hdfvs('detach',vdata_id);
  vdata_ref = hdfvs('getid',file_id,vdata_ref);
end

% line 436 of read_trmm.c
nparm = 0;
nparm = input('Enter total number of vdata to write out or 0 to exit program ');
if nparm == 0
   hdfml('closeall');
   return;
elseif (nparm == total_vidx)
  nindex = 1;
  vdata_index(nindex) = index+1;
  nindex = ndindex +1;
else
  temparm = 0;
  temparm = input('Enter parameter numbers from above list as a MATLAB array ');
  nindex = 0;
  for index  = 0:(nparm-1)
    if temparm <=total_vidx
      vdata_index = temparm;
      nindex = length(temparm);
    end
  end
end

if nparm > 0
  disp('The valid vdata numbers you Enter are: ');
  disp(vdata_index);
else
  disp('No vdata selected -- # of vdata = 0');
end

% Now loop through vdatas and write selected data(s) to binary files
if nindex > 0
  for index = 1:nindex
    vdata_idx = vdata_index(index)-1;
    vdata_id = hdfvs('attach',file_id,total_vref(vdata_idx+1),'r');
    
    [n_records,interlace,fields, vdata_size, vdata_name, status] = hdfvs('inquire',vdata_id);
    
    if status ~= 0
      error(sprintf('VSinquire failed on vdata %s\n', vdata_name));
    end
    
    fprintf('\n\n *************************************************\n');
    fprintf('   Vdata name     :   %s\n', vdata_name);
    fprintf('   Vdata recs     :   %d\n', n_records);
    fprintf('   Vdata recsize  :   %d\n', vdata_size);
    
    status = hdfvs('setfields',vdata_id,fields);
    [vdatabuf, status] = hdfvs('read',vdata_id,n_records);
    
    if status ~= n_records
      error('Vdata recs read and num_rec do not match');
    end
    
    ifield = hdfvf('nfields',vdata_id);
    
    disp('Field Name list: ')
    for field_index = 0:(ifield-1)
      fieldname = hdfvf('fieldname',vdata_id,field_index);
      field_size = hdfvf('fieldisize',vdata_id,field_index);
      field_order = hdfvf('fieldorder',vdata_id,field_index);
      field_type = hdfvf('fieldtype',vdata_id,field_index);
      fprintf('         %-20s  (size: %5d, type: %3s )\n', fieldname,field_size, field_type);
    end
    
    % Now write out the Vdata information to a binary file
    outfile = [filename '.' vdata_name];
    vd_fp = fopen(outfile,'wb','b');
    if vd_fp == -1
      error(sprintf('Open vdata file failed  %s\n', vdata_name));
    end
    
    count = fwrite(vd_fp,double(vdatabuf{1}),'uint8');
    fclose(vd_fp);
    status = hdfvs('detach',vdata_id);
    
  end % end of for statement for index of total selected vdata
end    % end of if statement for nindex 

status = hdfv('end',file_id);
status = hdfh('close',file_id);
