function d = read_bb(file, varargin)
%function d = read_bb(file, range, step, 'vel', ...)
%file is 'file0123.000';
%range is number of ensembles
%step is discretization of profiles to read
%'vel','cor','amp','ts','bt','all'


% 00/01/31 EF
% 00/02/22 JH (bb from read_os)
% 00/07/05 JH (fixed ipp bug (here and in read_os))  
%  *** NOTE: fix was insufficient - check read_os and fix the fix

step = 1;
start = -1;
 stop = 1000000000;
out.vel = 0; out.cor = 0; out.amp = 0; out.ts = 0; out.bt = 0;

for iarg = 1:(nargin-1),
   arg = varargin{iarg};
   if ischar(arg),
      switch lower(arg)
         case 'vel', out.vel = 1;
         case 'cor', out.cor = 1;
         case 'amp', out.amp = 1;
         case 'ts', out.ts = 1;
         case 'bt', out.bt = 1;
         case 'all', out.vel = 1; out.cor = 1; out.amp = 1;
                     out.ts = 1; out.bt = 1;
         otherwise
            error('invalid data type requested')
      end
   else
      nn = length(arg);
      if nn == 1,
         step = arg;
      elseif nn == 2,
         start = arg(1);
         stop = arg(2);
      else
         error('Invalid size of numeric argument')
      end
   end
end

fid = fopen(file, 'r', 'ieee-le');
[config, buf] = read_buf(fid);  % First time: read config.
if ischar(buf)
   d.error = buf;
   return
end
fseek(fid, 0, 'eof');
len_file = ftell(fid);
nbytes = config.nbytes + 2;
n_profiles = floor(len_file / nbytes);
first_profile = buf.ens_num;
last_profile = first_profile + n_profiles - 1;

start = max([start, first_profile]);
stop = min([stop, last_profile]);

NP = 1 + floor((stop - start) / step);
NC = config.ncells;
NB = 4; % beams

if out.vel, d.vel = zeros(NC, NB, NP); end
if out.cor, d.cor = zeros(NC, NB, NP); end
if out.amp, d.amp = zeros(NC, NB, NP); end
if out.ts,
   d.heading = zeros(NP, 1);
   d.pitch =      d.heading;
   d.roll =       d.heading;
   d.soundspeed = d.heading;
   d.temperature =d.heading;
end
if out.bt,
   d.bt.vel = zeros(NB, NP);
   d.bt.range = d.bt.vel;
   d.bt.cor =   d.bt.vel;
   d.bt.amp =   d.bt.vel;
%   d.bt.rssi =  d.bt.vel;
end

% Always include dday and ens_num.
d.dday = zeros(NP, 1);
d.ens_num = zeros(NP, 1);
d.config = config;

%NOTE: this loop assumes the configuration does not change in the file
%  being read
ipp = 1;
for ip = start:step:stop
   fseek(fid, nbytes * (ip - first_profile), 'bof');
   [config, buf] = read_buf(fid, config);
   if ischar(buf),
      d.error = buf;
      break
   end
   d.dday(ipp) = buf.dday;
   d.ens_num(ipp) = buf.ens_num;
   if out.ts,
      d.heading(ipp) = buf.heading;
      d.pitch(ipp)   = buf.pitch;
      d.roll(ipp)    = buf.roll;
      d.soundspeed(ipp)  = buf.soundspeed;
      d.temperature(ipp) = buf.temperature;
   end
   if out.vel, d.vel(:, :, ipp) = buf.vel; end
   if out.cor, d.cor(:, :, ipp) = buf.cor; end
   if out.amp, d.amp(:, :, ipp) = buf.amp; end
   if out.bt,
      d.bt.vel(:, ipp) = buf.bt.vel;
      d.bt.cor(:, ipp) = buf.bt.cor;
      d.bt.amp(:, ipp) = buf.bt.amp;
%      d.bt.rssi(:, ipp) = buf.bt.rssi;
   end
   ipp = ipp+1;
 end

if out.ts,
   d.heading = 0.01 * d.heading;
   d.pitch = 0.01 * to_signed(d.pitch);
   d.roll = 0.01 * to_signed(d.roll);
   d.temperature = 0.01 * to_signed(d.temperature);
end

if out.vel, d.vel = 0.001 * to_signed(d.vel); end

if out.bt,
   d.bt.vel = 0.001 * to_signed(d.bt.vel);
   d.bt.range = 0.01 * d.bt.range;  % BAD is 0; leave it.
end


fclose(fid)
d.error = [];



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% functions %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [config, buf] = read_buf(fid, config)

% This function is inefficient; among other things, it always
% reads the first few bytes and then goes back and reads the
% whole block.  It also calculates config whether it is likely
% to have changed or not.  (I think that the config items should
% always be constant in a single file.) It seems unlikely that
% it is worth optimizing, however.
%


if nargin == 1, config.yearbase = []; end  % avoids warning upon early error return
start = ftell(fid);
[ss, count] = fread(fid, 2, 'uint16');
if count ~= 2,
   buf = sprintf('error 1, only %d bytes at start of buffer\n', count);
   return
end
if ss(1) ~= 32639 % 7F7F;  %2 bytes of 7Fh each
   buf = 'error 2, Header ID/Data Source not found in read_buf';
   return
end
config.nbytes = ss(2);  % Total bytes up to, not including, checksum.
fseek(fid, start, 'bof'); %go back to beginning of file
[rec, count] = fread(fid, config.nbytes, 'uint8');
if count ~= config.nbytes,
   buf = sprintf('error 3, only %d bytes read, %d expected\n', ...
                     count, config.nbytes);
   return
end
test_cs = mod(sum(rec), 65536);
cs = fread(fid, 1, 'uint16');
if test_cs ~= cs,
   buf = ['error 4, Checksum failure; CS is ', int2str(cs), ' versus calculated ', ...
            int2str(test_cs)];
   return
end

if nargin == 1,
   config.ndata = rec(6);
   ii = 7 + 2 * (0:(config.ndata-1));
   config.offsets = rec(ii) + 256 * rec(ii+1);
   config.ids = rec(config.offsets + 1) + rec(config.offsets + 2) * 256;

%   %%% fixed leader:
%   lsb_base =                 5; 1           %least significant base
%   msb_base =                 6; 1           %most significant base
%   nbeams_used_base =         9; 1
%   ncells_base =             10; 1
%   pings_per_ens_base =      11; 2
%   cell_base =               13; 2
%   blank_base =              15; 2
%   mode_base =               17; 1
%   code_reps_base =          19; 1
%   coord_transf_base =       26; 1
%   heading_align_base =      27; 2
%   heading_bias_base =       29; 2
%   EZ_base =                 31; 1
%   bin1distance_base =       33; 2
%   pulse_base =              35; 2
%   lag_base =                41; 2
      
   ofs = config.offsets(find(config.ids == 0));
	
   config.nbeams_used =        rec(ofs + 9);
   config.ncells =             rec(ofs + 10);
   config.pings_per_ens_base = 0.01 * (rec(ofs + 11) + 256 * rec(ofs + 12));
   config.cell =               0.01 * (rec(ofs + 13) + 256 * rec(ofs + 14));
   config.blank =              0.01 * (rec(ofs + 15) + 256 * rec(ofs + 15));
   config.mode =               rec(ofs + 17);
   config.code_reps =          rec(ofs + 19);
   config.blank =              0.01 * (rec(ofs + 15) + 256*rec(ofs + 16));
   config.coord_transf =       rec(ofs + 26);
   config.heading_align =      0.01 * (rec(ofs + 27) + 256 * rec(ofs + 28));
   config.heading_bias =       0.01 * (rec(ofs + 29) + 256 * rec(ofs + 30));
   config.EZ =                 rec(ofs + 31);
   config.bin1distance =       0.01 * (rec(ofs + 33) + 256 * rec(ofs + 34));
   config.pulse =              0.01 * (rec(ofs + 35) + 256 * rec(ofs + 36));
   config.lag =                0.01 * (rec(ofs + 41) + 256 * rec(ofs + 42));

   lsb = rec(ofs + 5);
   msb = rec(ofs + 6);
   config.sysconfig = decode_sysconfig(lsb, msb);
end

%%% variable leader:
ofs = config.offsets(find(config.ids == 128));
buf.ens_num = rec(ofs + 3) + 256 * rec(ofs + 4) + 65536 * rec(ofs + 12);

tt = rec(ofs + (5:11));
tt(1) = tt(1) + 1900;
if (tt(1) < 1980)  tt(1) = tt(1) + 100; end
% Fix later: substitute yearbase for tt(1)

if (isempty(config.yearbase)), config.yearbase = tt(1); end
buf.dday = to_day( config.yearbase, tt(1:6).') + tt(7) / 8640000;
buf.ymdhmsh = tt;

buf.soundspeed =     rec(ofs + 15) + 256 * rec(ofs + 16);
buf.heading =        rec(ofs + 19) + 256 * rec(ofs + 20);
buf.pitch   =        rec(ofs + 21) + 256 * rec(ofs + 22);
buf.roll    =        rec(ofs + 23) + 256 * rec(ofs + 24);
buf.temperature =    rec(ofs + 27) + 256 * rec(ofs + 28);

% In the following blocks, the "min" should not be needed
% when finding the offsets; we may eliminate it.
% The empty check is needed, however.

%%% velocity:
ofs = min(config.offsets(find(config.ids == 256)));
ii = 2 + (1:2:(config.ncells*8));
if isempty(ofs),
   a = NaN * ii;
else
   a = rec(ofs + ii) + 256 * rec(ofs + 1 + ii);
end
buf.vel = reshape(a, 4, config.ncells).';
% We are deferring the handling of negative numbers and NaNs.


%%% correlation:
ofs = min(config.offsets(find(config.ids == 512)));
ii = 2 + (1:(config.ncells*4));
if isempty(ofs)
   buf.cor = reshape(ii, 4, config.ncells).';
else
   buf.cor = reshape(rec(ofs + ii), 4, config.ncells).';
end

%%% amp:
ofs = min(config.offsets(find(config.ids == 768)));
ii = 2 + (1:(config.ncells*4));
if isempty(ofs)
   buf.amp = reshape(ii, 4, config.ncells).';
else
   buf.amp = reshape(rec(ofs + ii), 4, config.ncells).';
end


%%% percent good:
ofs = min(config.offsets(find(config.ids == 1024)));
ii = 2 + (1:(config.ncells*4));
if isempty(ofs)
   buf.pg = reshape(ii, 4, config.ncells).';
else
   buf.pg = reshape(rec(ofs + ii), 4, config.ncells).';
end

%%% status starts at 1280


%% BT:
ofs = min(config.offsets(find(config.ids == 1536)));
if isempty(ofs)

   config.bt.BP = NaN;
   config.bt.BC = NaN;
   config.bt.BA = NaN;
   config.bt.BX = NaN;

   ii = (0:2:6);
   buf.bt.range = NaN*ii;
   buf.bt.vel =   NaN*ii;
   buf.bt.cor =   NaN*ii;
   buf.bt.amp =   NaN*ii;
   buf.bt.rssi =  NaN*ii;

else
   config.bt.BP = rec(ofs + 3) + 256 * rec(ofs + 4);
   config.bt.BC = rec(ofs + 7);
   config.bt.BA = rec(ofs + 8);
   config.bt.BX = rec(ofs + 71) + 256 * rec(ofs + 72);

   ii = (0:2:6); iii = 0:3;
   buf.bt.range = rec(ofs + ii + 17) + 256 * rec(ofs + ii + 18) + ...
                  65536 * rec(ofs + iii + 78);

   buf.bt.vel = rec(ofs + ii + 25) + 256 * rec(ofs + ii + 26);
   buf.bt.cor = rec(ofs + iii + 33);
   buf.bt.amp = rec(ofs + iii + 37);
%   buf.bt.rssi = rec(ofs + iii + 73);
end


%-----------------------------------------------------------------------------
function c = decode_sysconfig(l, m)

%%% was in read_bb.m
  
freqs = [75, 150, 300, 600, 1200, 2400, 38];
freqbits = bitand(l, 7);
c.frequency = freqs(1 + freqbits);
c.up =  bitand(bitshift(l, -7), 1);
c.convex = bitand(bitshift(l, -3), 1);

angles = [15, 20, 30, NaN];
anglebits = bitand(m, 3);
c.angle = angles(1 + anglebits);

