function [a, s] = read_bin(f, varargin)
%
% Read binary files of the type written by serasc2bin.py.
%
% [a, s] = read_bin(filename, varargin)  or
% [a, s] = read_bin(s, varargin)
%
% filename is the full file name of the binary file to be read.
% s is a structure including s.fid, the file id of an open file,
%      and s.irange, the index range of the records in a.
% a is the output data array; each column is a record.
%
% If only one output argument is requested, the file is closed
% on exit.
%
% With no keyword arguments, the entire file is read. Otherwise,
% use one of the following keyword argument pairs to select
% the records to be returned.
%
%     'last', 5            return last 5 records
%     'irange', [1 20]     return records 1-20
%
% EF, JH 2001/04/11; modified 04/15
%

if ischar(f),
   filename = f;
   s.filename = f;
   s.fid = fopen(filename, 'rb');
   if s.fid <= 0,
      fprintf('returning empty: Cannot open %s\n', filename)
      a = [];
      return
   end
   [s.nrows, s.offset, s.rows] = read_header(s.fid);
   s.last = 0;
   s.irange = [];
   s = fillstruct(s, varargin);
elseif isstruct(f),
   if isempty(varargin),
      s = f;
   else
      f.last = 0;
      f.irange = [];
      s = fillstruct(f, varargin);
   end
else
   error(['First argument must be a filename, or a structure\n', ...
          'from a previous call to read_bin_incr'])
end

% Ensure that we read only complete records, if the file is
% still being written.
fseek(s.fid, 0, 'eof');
offset_end = ftell(s.fid);
% It seems that Matlab file handling is unreliable with
% seeks from varying locations, so we will close and
% reopen the file to be on the safe side.
fclose(s.fid);
s.fid = fopen(s.filename);
ncol_all = floor((offset_end - s.offset) / (8 * s.nrows));
offset_end = s.offset + ncol_all * 8 * s.nrows;

if s.last == 0,
   if isempty(s.irange)
      offset = s.offset;
      ncol = ncol_all;
   else
      offset = s.offset + (s.irange(1) - 1) * s.nrows * 8;
      ncol = 1 + diff(s.irange);
      ncol = min([ncol, (ncol_all - s.irange(1) + 1)]);
   end
else
   offset = offset_end - 8 * s.nrows * s.last;
   offset = max([offset, s.offset]);
   ncol = s.last;
end


fseek(s.fid, offset, 'bof');
[a, count] = fread(s.fid, [s.nrows, ncol], 'float64');
s.irange(1) = 1 + (offset - s.offset) / (8 * s.nrows);
s.irange(2) = s.irange(1) - 1 + count/s.nrows;

if nargout <= 1,
   fclose(s.fid);
end

%---------------------------------------------------
function [nrows, offset, rownames] = read_header(fid)

line1 = fgetl(fid);
[A, count, errmsg, nextindex] = sscanf(line1, '%d', 2);
nrows = A(1);
nlines = A(2);
rownames = cell(nrows, 1);
source = sscanf(line1(nextindex:end), '%s');
for ii = 2:nlines,
   line = fgetl(fid);
   rownames{ii-1} = sscanf(line, '%s');
end
offset = ftell(fid);

