function d = read(instclass, file, varargin)
% d = read(instclass, file, varargin)
%
% assume varargin is a collection of name, values pairs
% aquisition of single data types  in a cellarray called 'vars'
% include any or all in the cellarray to turn them on
% 'vars',{'vel', 'sw', 'amp', 'pg', 'ts', 'bt', 'status' 'all', 'nav'}
% -- or -- if there is only one variable, it can be a string
% 'vars', 'ts'
%
% (unlike instclass "os" navigation does not exist here, but we can make
%  the structure have the fields
% 
% 'yearbase', [yearbase]             % REQUIRED  (so 'yearbase' field is ok)
%                                    % dday is calculated with RDIyearbase=1997
%                                    % (you have to convert it)
% 'irange',  [istart istop]          % range of indices to use
% 'ilist', [ index list]             % explicit index list overrides step
% 'ens_range'  [ens_start ens_end]    % ens_num range
% 'ens_list', [ens_num list]         % explicit list of ens_num
% 'step', [step]                     % use only every "step" index
%                                    % (do not use with ens_list or ilist)
% 'ends', [choose from 0,1,2,4]      % 0 (default - don't do it)
%                                    % 1: get only first profile
%                                    % 2: get first and last profile
%                                    % 3: not implimented
%                                    % 4: get first and last pair of profiles
% 'add_to_ens' [add to ens_num(1)]   % add this ens_num before matching
% 'RDI_ens'     [0 1]                % output original RDI ens (default = 0;
%                                    %   i.e. default is to unwrap at 65536)
% 'repl_RDIyearbase', [yearbase]     % use thus yearbase for dday calculation
%                                    % instead of 1997 (RDIyearbase).  Only 
%                                    % useful if adcp clock set well enough
%                                    % enough to use as a time base (rare)
% 'verbose', [0 or 1]                % be verbose. default = 0 (not verbose)
%
%  NOTE on ranges: default is 1:end
%       You may only use one of the following choices.  If you
%       have specified too much (i.e. have conflicting information) read will
%       exit.  Choices are:
%       -  irange (step is optional, defaults to 1)
%       -  ilist
%       -  ens_range (step is optional, defaults to 1)
%       -  ens_list 
%       -  ends  
% 
%  
%  NOTE on ens_num: this is a nondecreasing ensemble count based on the
%       ensemble number (RDIens_num) read directly from the raw file, but 
%       takes into account the 2^16 wrap that may occur within the file in 
%       RDIens_num.  It is assumed that RDIens_num will only decrease within
%       a file because of the 2^16 wrap, and it decreases between files if and
%       only if it resets to 1 because of a possible configuration change.
%       


% 00/01/31 EF
% 01/01/26 JH - add instclass, work arguments to match the os template
%             - adding "add_to_ens" to varargin and as an output structure
%                 so now ens_num will always increase (usually by 1) and 
%                 will start with ens_num(1)+add_to_ens 

%% undocumented features:
%% 'vars' defaults to none, turn on pieces or 'all' exlicitly
%% range defaults to all; choose specific subsets explicitly

% should make it so all cfg fields are initialized - then can use fillstruct(x,x,0)
%    to help with debugging (i.e. fail if a new field is introduced)

if (nargin < 2)
   help read 
   error('must set instclass and file')
end


start = -1;          
stop = 1000000000;
RDIyearbase = 1997;  %a non-leap year for RDI's calendar
                     %NOTE: we're going to use this to guarentee a 1s/s
                     %dt in the RDI dday; the config.yearbase will
                     %simply fill the data.yearbase slot

%set defaults.  all except field "vars" are present in the original
%  cfg structure; fillstruct will update the values listed
%  Field "vars" will be initialized by fillstruct.  set those chosen to 1
cfg.RDI_ens   = 0;  % do not unwrap ens_num (maybe use in scan?)
cfg.add_to_ens= 0;  % add this value (integer >= 0) to all ens_num
cfg.yearbase  = 0;  % it it's still 0 after varargin, abort.  
cfg.repl_RDIyearbase  = [];  % use RDIyearbase = 1997 
cfg.step      = 1;
cfg.verbose   = 0;  % not verbose
cfg.irange    = []; % default all range specifications to [] (see function
cfg.ilist     = []; % default all range specifications to [] test_range_request
cfg.ens_range = []; % default all range specifications to [] for error checking)
cfg.ens_list  = []; % default all range specifications to []
cfg.ends      = []; % default all range specifications to []
cfg.vars      = []; % default variables to []
vars          = {'vel', 'sw', 'amp', 'pg', 'ts', 'bt', 'status' 'all'};
out.nav       = 0; %initialized "out" as a nonempty structure AND
                   % separate 'nav' from the rest of 'all' (since 
                   % there may not be any)

% structure "out" fields from "vars" cellarray initialized as zero
for cii = 1:length(vars)
   out = setfield(out, vars{cii}, 0); 
end


% update cfg.
if (nargin > 2)%   we'll assume the rest is a set of name,value pairs
               %override input structure with name,value pairs
   cfg = fillstruct(cfg, varargin, 0 );
end


if (length(cfg.irange) > 2)
   fprintf('---> ''irange'' was a list.')
   fprintf('using index endpoints [%d %d] for range.\n', ...
           cfg.irange(1), cfg.irange(end))
   cfg.irange = cfg.irange([1 end]);
end
if (length(cfg.ens_range) > 2)
   fprintf('---> ''ens_range'' was a list.')
   fprintf('using ensemble endpoints [%d %d] for range.\n', ...
           cfg.ens_range(1), cfg.ens_range(end))
   cfg.ens_range = cfg.ens_range([1 end]);
end
   


%take cfg it apart into variables
% first look at the variables that already exist
if (cfg.yearbase == 0)
   error('must set yearbase')
end
if (cfg.yearbase < 70)  %they used a 2-digit year
   cfg.yearbase = cfg.yearbase + 1900;
   warning('you should use a 4-digit yearbase for clarity')
end

if ~isempty(cfg.repl_RDIyearbase)
   RDIyearbase = cfg.repl_RDIyearbase;
end

%now check whether cellarray 'vars' has been added - use the elements of
% cfg.vars (if it exists) to turn on aquisition of the listed fields
% NOTE: cellarray "vars" does not have "nav", whereas cfg.vars could.
if (isfield(cfg, 'vars'))
   if (ischar(cfg.vars))  %if there is only variable chosen, allow a string
      cfg.vars = {cfg.vars};
   end
   
   for cii = 1:length(cfg.vars)
      out = setfield(out, cfg.vars{cii}, 1); 
   end
   
   % out.all (override individual chosen data types and output everything)
   if ismember({'all'}, cfg.vars)
      %they turned it on so get everything
      for cii = 1:length(vars)
         out = setfield(out, vars{cii}, 1);
      end
   end
   out = rmfield(out, 'all');  %now remove: there shouldn't be an out.all
end



%fail if range is decreasing or if more than one option is chosen
doit = test_range_request(cfg);

if (~isstruct(doit)) % otherwise it would be a structure with doit.irange
   error('failed range request')
end

%this is quick
if (~exist(file, 'file'))
   error(sprintf('cannot open %s\n', file))
end

params.yearbase = cfg.yearbase;
params.RDIyearbase = RDIyearbase;
params.out = out;

fid = fopen(file, 'r', 'ieee-be');  % the OS processor writes little-endian
[config, buf] = read_buf(fid, params);  % 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_indexes = floor(len_file / nbytes);
first_index = 1; %(was) buf.ens_num;
first_RDIens = buf.ens_num;
last_index = first_index + n_indexes - 1;
fclose(fid); %we might read it below in read_rawcfg


% revise what the last good index
fid = fopen(file, 'r', 'ieee-be');  % the OS processor writes little-endian
for last_good_index = last_index:-1:first_index
   buf = [];
   fseek(fid, nbytes * (last_good_index - 1), 'bof'); %use ilist
   if (cfg.verbose), last_good_index, end
   try
      [config, buf] = read_buf(fid, params, config);
   catch
      if (cfg.verbose)
         fprintf('exited read early at index %d (first: %d, last: %d):\n', ...
                 last_good_index, first_index, last_index)
         fprintf('%s\n',d.error)
      end
   end 
   if isstruct(buf)
      break
   end
end
fclose(fid);
% update these things
last_index = last_good_index;
n_indexes = last_index - first_index + 1; 
last_RDIens = buf.ens_num;
if (last_RDIens < first_RDIens)
   last_RDIens = last_RDIens + 65536;  
end

%nothing was chosen; do it all (doit.sum should be 0 for this)
% set the default anyway so we can subsample it when we do ens_list
% override if any of the other options is chosen
% NOTE: ilist is an INDEX list, not an ens_num list (conversion takes
%   place later)

ilist = first_index:cfg.step:last_index;
%use ens_num this to find indices: assumes ens_num always goes up by 1
ens_num = first_RDIens + [0:n_indexes-1] + cfg.add_to_ens;
if (cfg.verbose & (cfg.add_to_ens ~= 0 ))
   fprintf('adding %d to ensembles before finding matching indices\n')
end



%this was an explicit index list
if (doit.ilist)
   ilist = cfg.ilist;
   badi = find(ilist > last_index);
   if (~isempty(badi))
      ilist(badi)=[];
      if (cfg.verbose)
         fprintf('ilist excedes last index: truncating after %d\n',last_index)
      end
   end
end

% this was an explicit index range
if (doit.irange)
   ilist = [cfg.irange(1) : cfg.step : min(cfg.irange(end), last_index)]';
   if (cfg.verbose & (cfg.irange(end)>last_index))
      fprintf('irange excedes last index: truncating after %d\n',last_index)
   end
end
 
% just do the ends
if (doit.ends) 
   out.ends = cfg.ends;
   istart = max([start, first_index]);
   istop = min([stop, last_index]);
   
   if (out.ends & (n_indexes >= out.ends))
      if (out.ends == 4)
         ilist = [first_index, (first_index + 1), ...
                  (last_index - 1), last_index];
      elseif (out.ends == 2)
         ilist = [first_index last_index];
      else %it has to be 1 because we checked it before
         ilist = [first_index];
      end
   elseif isempty(ilist)
      ilist = istart:step:istop;
   end
end

%now we're left with ens_range or ens_list
%start by getting the actual final ensemble list (unless it's specified)

if (doit.ens_list)
   ilist = ilist(ismember(ens_num, cfg.ens_list));
end

if (doit.ens_range)
   %want to use the appropriately wrapped ensemble numbers if possible
   
   cfg.ens_range(1) = max(cfg.ens_range(1), ens_num(1));
   cfg.ens_range(2) = min(cfg.ens_range(2), ens_num(end));
   first_index = find(cfg.ens_range(1) == ens_num);
   last_index =  find(cfg.ens_range(end) == ens_num);
   ilist = [first_index : cfg.step : last_index]';
   if (cfg.verbose)
      fprintf('starting ens (index) %d (%d); ending ens (index) %d (%d)\n',...
              cfg.ens_range(1), first_index, cfg.ens_range(end), last_index)
   end
end


%initialize everything; may have to stop prematurely if the end of the file
%was badly terminated (may have to change NP later)
NP = length(ilist); 
NC = config.ncells;
NB = 4; % beams


%initialize structure
if out.vel,    d.vel = zeros(NC, NB, NP);    end
if out.sw,     d.sw = zeros(NC, NB, NP);     end
if out.amp,    d.amp = zeros(NC, NB, NP);    end
if out.pg,     d.pg = zeros(NC, NB, NP);     end
if out.status, d.status = zeros(NC, NB, NP); end
if out.ts,
   d.heading = zeros(NP, 1);
   d.pitch =      d.heading;
   d.roll =       d.heading;
   d.temperature =d.heading;
   d.HiV =        d.heading;
   d.XmtAmps =    d.heading;
end
if out.bt,
   d.bt.vel = zeros(NB, NP);
   d.bt.range = d.bt.vel;
   d.bt.pg =   d.bt.vel;
end
if out.nav,
   d.nav.sec_pc_minus_utc = zeros(NP, 1);
   d.nav.txy1 = zeros(3, NP);
   d.nav.txy2 = d.nav.txy1;
end


d.add_to_ens = cfg.add_to_ens; %initialize  from default or  varargin 

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

%changes here:
% (1) r.ens_num will always be (unwrapped) r.RDIens_num + add_to_ens
% (2) if r.ens_num == 65535 is about to wrap 
% (3) there is no checking for other drops in RDIens_num 


fid = fopen(file, 'r', 'ieee-be');  % the OS processor writes little-endian
fseek(fid,0, 'bof'); % start at the beginning of the file
[config, buf] = read_buf(fid, params, config);

add_for_wrap = 0;
for ii = 1:length(ilist);
fseek(fid, nbytes * (ilist(ii) - 1), 'bof'); %use ilist
   [config, buf] = read_buf(fid, params, config);
   
   d.dday(ii) = buf.dday;
   RDIens_num = buf.ens_num;
   %unwrap it:
   if (buf.ens_num + d.add_to_ens < d.ens_num(1)) %then we passed a wrap
      add_for_wrap = (~cfg.RDI_ens * 65536);
   end
   d.ens_num(ii) = RDIens_num + d.add_to_ens + add_for_wrap;  
   
   if out.ts,
      d.heading(ii) = buf.heading;
      d.pitch(ii)   = buf.pitch;
      d.roll(ii)    = buf.roll;
      d.temperature(ii) = buf.temperature;
      d.HiV(ii) = buf.HiV;
      d.XmtAmps(ii) = buf.XmtAmps;
   end
   if out.vel & config.n.vel, d.vel(:, :, ii) = buf.vel.'; end
   if out.sw & config.n.sw, d.sw(:, :, ii) = buf.sw.'; end
   if out.amp & config.n.amp, d.amp(:, :, ii) = buf.amp.'; end
   if out.pg & config.n.pg, d.pg(:, :, ii) = buf.pg.'; end
   if out.status & config.n.status, d.status(:, :, ii) = buf.status.'; end
   if out.bt
      d.bt.vel(:, ii) = buf.bt.vel;
      d.bt.range(:, ii) = buf.bt.range;
      d.bt.pg(:, ii) = buf.bt.pg;
   end
end


%if (brokeout) %we broke out early
%   fprintf('exited read early at %d (first = %d, last = %d):\n', ...
%           ilist(ii-1), ilist(1), ilist(end))
%   fprintf('%s\n',d.error)
%   d = cutadcp(instclass, d, [1:ii-1]);
%   d.error = sprintf('%s: nb read exited read early %d (out of %d)',...
%                     d.error,  ilist(ii-1), ilist(end));
%end


if out.ts,
   d.heading = d.heading * 360/65536;
   d.pitch = d.pitch * 360/65536;
   d.roll = d.roll * 360/65536;
   d.temperature = 45 - (50 * d.temperature / 4096);
   d.HiV = d.HiV * 0.17;
   d.XmtAmps = d.XmtAmps * 0.02;
end

if out.vel, d.vel = config.vel_multiplier * d.vel; end

if out.bt,
    %% hardwired for beam--temporary
    %%this is beam coordinates so the config.vel_multiplier is hardwired
    %%       for beam coord      %% should be OK for shipboard use
    %% maybe need to change this depending on beam coords or geographic coords
   d.bt.vel = config.vel_multiplier * to_signed(d.bt.vel); 
   d.bt.pg = 6.67 * d.bt.pg;
end


%%%% add things which will make our lives easier later:
d.config.bin1distance = config.blank + (config.cell + config.pulse)/2;
d.config.tr_depth = 0;        % unknown in nb data; set to 0
d.num_pings = ones(size(d.dday));% only one ping per "ensemble"
d.soundspeed = 1530 + zeros(size(d.dday));  %any field would do

%gol'durnit - I'm going to do this here too.
index = [1:length(d.dday)]';
d.dday = d.dday + index*10^(-7); %now it's monotonic


fclose(fid);
d.error = [];


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


RDIyearbase = params.RDIyearbase;
yearbase = params.yearbase;
out = params.out;


if nargin == 2, config.yearbase = yearbase; end % initializes config

start = ftell(fid);
[hdr, count] = fread(fid, 7, 'uint16');

if count ~= 7,
   buf = sprintf('error 1, only %d bytes at start of buffer\n', count);
   return
end
config.nbytes = hdr(1);
config.n.leader = hdr(2);
config.n.vel = hdr(3);
config.n.sw = hdr(4);
config.n.amp = hdr(5);
config.n.pg = hdr(6);
config.n.status = hdr(7);

fseek(fid, start, 'bof');
[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 == 2,
   config.ofs.leader = 14;
   config.ofs.vel = config.ofs.leader + config.n.leader;
   config.ofs.sw = config.ofs.vel + config.n.vel;
   config.ofs.amp = config.ofs.sw + config.n.sw;
   config.ofs.pg = config.ofs.amp + config.n.amp;
   config.ofs.status = config.ofs.pg + config.n.pg;
   config.ndata = 5;

   %%% leader:
   ofs = config.ofs.leader;
   config.nbeams_used = 4;
   config.ncells = rec(ofs + 11);
   config.cell = 2 .^ (rec(ofs + 12));
   config.blank = rec(ofs + 14);
   config.pulse = rec(ofs + 13);
   config.sysconfig = decode_sysconfig(rec(ofs + 19));
   if config.sysconfig.hirange,   x = 0.25;
   else x = 0.125;
   end
   if config.sysconfig.earthref, x = x * 2; end
   if config.sysconfig.frequency == 75 & config.sysconfig.hirange == 0,
      x = x * 2;
   end
   config.vel_multiplier = x * 0.01;  % m/s per count
end


%%% variable part of leader:
ofs = config.ofs.leader;
buf.ens_num = rec(ofs + 17) + 256 * rec(ofs + 16);

% Use 1997 to get a 1s/s dday from RDI
tt = decode_time(RDIyearbase, rec(ofs + (1:5))).';

%??? is this wrong?
buf.dday = to_day(RDIyearbase, tt);
buf.ymdhms = tt;

buf.heading = (rec(ofs + 27) + 256 * rec(ofs + 26));
buf.pitch   = (rec(ofs + 23) + 256 * rec(ofs + 22));
buf.roll    = (rec(ofs + 25) + 256 * rec(ofs + 24));

buf.temperature = rec(ofs + 29) + 256 * rec(ofs + 28);
buf.HiV = rec(ofs + 30);
buf.XmtAmps = rec(ofs + 31);

% bottom track
buf.bt.vel = decode_vel(rec(ofs + (42:47)));
ibad = find(buf.bt.vel == 2048);
buf.bt.vel(ibad) = NaN;
buf.bt.range = rec(ofs + (48:2:54)) * 256 + rec(ofs + (49:2:55));
buf.bt.range(ibad) = NaN;
a = rec(ofs + (62:63));
buf.bt.pg = zeros(4, 1);
buf.bt.pg([1, 3]) = bitshift(a, -4);
buf.bt.pg([2, 4]) = bitand(a, 15);

%%% velocity:
ofs = config.ofs.vel;
buf.vel = decode_vel(rec(ofs + (1:config.n.vel)));

if config.n.status ~= 0,
   ofs = config.ofs.status;
   a = rec(ofs + (1:config.n.status));
   buf.status = zeros(4, config.ncells);
   buf.status(1,:) = bitshift(a(1:2:end), -4);
   buf.status(2,:) = bitand(a(1:2:end), 15);
   buf.status(3,:) = bitshift(a(2:2:end), -4);
   buf.status(4,:) = bitand(a(2:2:end), 15);
   ibad = find(buf.status ~= 0);
   buf.vel(ibad) = NaN;
else
   ibad = find(buf.vel == 2048);
   buf.vel(ibad) = NaN;
end

%%% amp:
if config.n.amp ~= 0,
   ofs = config.ofs.amp;
   a = rec(ofs + (1:config.n.amp));
   buf.amp = reshape(a, 4, config.ncells);
end

%%% spectral width:
if config.n.sw ~= 0,
   ofs = config.ofs.sw;
   a = rec(ofs + (1:config.n.sw));
   buf.sw = reshape(a, 4, config.ncells);
end


%%% percent good:
if config.n.pg ~= 0,
   ofs = config.ofs.pg;
   a = rec(ofs + (1:config.n.pg));
   buf.pg = reshape(a, 4, config.ncells);
end


%--------------------------------------------------------
function c = decode_sysconfig(b)
%specific to nb

freqs = [75, 150, 300, 600, 1200, 115, 999];
freqbits = bitshift(bitand(b, 16+32+64), -4);
c.frequency = freqs(1 + freqbits);
c.hirange = bitand(b, 1);
c.earthref = bitand(b, 2);
c.up = bitand(b, 4);
c.concave = bitand(b, 8);


%---------------------------------------------
% function s = to_signed(u)
%----------------------------------------------
function td = decode_time(yrbase, t)
% t is an array with 5 "packed binary coded decimal"
% numbers: month, day, hour, minute, second
upper = bitshift(t, -4);
lower = bitand(t, 15);
td = upper * 10 + lower;
td = [yrbase; td(:)];


%-----------------------------------------------
function v = decode_vel(a)
n = length(a)/6; % 6 bytes for 4 beams
v = zeros(4,n);
v(1,:) = (a(1:6:end) * 16 + bitshift(a(2:6:end), -4)).';
v(2,:) = (bitand(a(2:6:end), 15) * 256 + a(3:6:end)).';
v(3,:) = (a(4:6:end) * 16 + bitshift(a(5:6:end), -4)).';
v(4,:) = (bitand(a(5:6:end), 15) * 256 + a(6:6:end)).';
% Figure out negative versus positive:
%ineg = find(bitshift(buf.vel, -11) == 1);
%buf.vel = bitand(buf.vel, 2047);  % Chop off sign bit.
ineg = find(v > 2048);
v(ineg) = v(ineg) - 4096;


%--------------------------------------------------
function doit = test_range_request(cfg)

%test whether everything is consistent: bomb out otherwise
%
% ... from help:
%
%  NOTE on ranges: you may only use one of the following choices.  If you
%       have specified too much (i.e. have conflicting information) read will
%       exit.  Choices are:
%       -  irange (step is optional)
%       -  ilist
%       -  ens_range (step is optional)
%       -  ens_list 
%       -  ends  
% 

%defaults were [] except cfg.step (which defaulted to 1)
fnames = {'irange', 'ilist', 'ens_range', 'ens_list', 'ends'};

doit.sum = 0; 
for ii = 1:length(fnames)
   val = getfield(cfg, fnames{ii});
   if (~isempty(val))
      doit = setfield(doit,  fnames{ii}, 1);
      doit.sum = doit.sum + 1;
   else
      doit = setfield(doit,  fnames{ii}, 0);
   end
end

%now doit is a structure so we can say "if (doit.ends)"  
% if there is an error, rewrite doit as the number 0

if (doit.sum > 1);
   fprintf('more than one range option chosen. avoid ambiguity: use only one\n')
   doit = 0;
end

   
if (~ismember(cfg.ends, [0 1 2 4]))
   fprintf(' option "ends" must be 0,1,2, or 4')
   doit == 0;
end

if (diff(cfg.ens_range)<=0)
   fprintf('ens_range must increase\n')
   doit = 0;
end


if (diff(cfg.irange)<=0)
   fprintf('irange must increase\n')
   doit = 0;
end




