function d = read(instclass, file, varargin)
% d = read(instclass, file, varargin)
%
% assume varargin is a collection of name, values pairs
% acquisition 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'}
% -- or -- if there is only one variable, it can be a string
% 'vars', 'ts'
%
% (unlike instclass "os" navigation does not exist here)
% 
% 'yearbase', [yearbase]             % - REQUIRED  (so 'yearbase' field is ok
%                                    %    and to calculate decimal day)
%
% %% accessing the data : default is [1 : step : end]
% 'irange',  [istart istop]          % - range of indices to use
% 'ilist', [ index list]             % - explicit index list overrides step
% 'step', [step]                     % - use only every "step" index
% '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
%
% '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:
%       -  ilist
%       -  irange (step is optional, defaults to 1) - overrides ilist
%       -  ends  - overrides ilist and irange
% 
%
%  NOTE ens_num is the raw ensemble value contained in the data.  There could
%       potentially be gaps in ens_num, but it should always increase (with 
%       the single exception of a wrap which occurs at 2^16)
%
%  NOTE a structure called "info" in the data structure contains
%       the base of the file read (info.rawfilebase) and the index 
%       list used to get the data (info.ilist)
%

% 00/01/31 EF
% 01/01/26 JH - add instclass, work arguments to match the os template,
%               remove references to ens_num except as value to output;
%               

%% 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 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


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

%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.yearbase  = 0;  % it it's still 0 after varargin, abort.  
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.ends      = []; % default all range specifications to []
cfg.vars      = []; % default variables to []
vars          = {'vel', 'sw', 'amp', 'pg', 'ts', 'bt', 'status' 'all'};
out.ts        = []; %have to initialize the structure

% 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

yearbase = fix_yearbase(cfg);
cfg.yearbase = yearbase;

%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
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

params.yearbase = cfg.yearbase;
params.out = out;
[ilist, config, nbytes, last_index] = get_ilist(file, cfg, params);


%initialize everything
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


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


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);

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;
   d.ens_num(ii) =  buf.ens_num;
   
   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 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 = 1470 + zeros(size(d.dday)); 

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

%augment d with some useful information
d.info.ilist = ilist;
[fpath, rawfilebase, ext] = fileparts(file);
d.info.rawfilebase = rawfilebase;

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

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);

tt = decode_time(yearbase, rec(ofs + (1:5))).';

%??? is this wrong?
buf.dday = to_day(yearbase, 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 [ilist, config, nbytes, last_index] = get_ilist(file, cfg, params)
% trap for obvious errors
% choosing 'ends' overrides ilist or irange; irange overrides ilist
% 'step' implimented for irange 

%(1)  look for errors
if (~ismember(cfg.ends, [0 1 2 4]))
   error(' option "ends" must be 0,1,2, or 4')
end


if (diff(cfg.irange)<=0)
   error('irange must increase\n')
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

% (2) read the file; look for errors

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;
last_index = first_index + n_indexes - 1;
fclose(fid); 


% revise 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))
      if (cfg.verbose)
         fprintf('buf returns %s\n',buf)
      end
   elseif isstruct(buf)
      break
   else
      warning('buf is not struct or string')
   end
end
fclose(fid);

% update these things
last_index = last_good_index;

n_indexes = last_index - first_index + 1; 

% default
ilist = first_index:cfg.step:last_index;

if (~isempty(cfg.ilist))  %this was an explicit index list
   ilist = cfg.ilist;
end

% this was an explicit index range
if (~isempty(cfg.irange))
   ilist = [cfg.irange(1) : cfg.step : min(cfg.irange(end), last_index)]';
end
 
% just do the ends
if (~isempty(cfg.ends) & (cfg.ends > 0))
   if (cfg.ends == 4)
      ilist = sort(unique([first_index, (first_index + 1), ...
                      (last_index - 1), last_index]));
   elseif (cfg.ends == 2)
      ilist = sort(unique([first_index last_index]));
   else %it has to be 1 because we checked it before
      ilist = [first_index];
   end
end


%toss out indices that are too small
badi = find(ilist < first_index);
if (~isempty(badi))
   ilist(badi)=[];
   if (cfg.verbose)
      fprintf('ilist smaller than first index: truncating before %d\n',...
              first_index)
   end
end

%toss out indices that are too large
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



if (isempty(cfg.ilist) & ...
    (isempty(cfg.ends)) & ...
    (isempty(cfg.irange)) & ...
    (isempty(ilist)))
   error('no indices found')
end


%----------------------------------------
function fixed_yearbase = fix_yearbase(cfg)
%take cfg it apart into variables
% first look at the variables that already exist

fixed_yearbase = cfg.yearbase;
% now override if necessary
if (cfg.yearbase == 0)
   error('must set yearbase')
end
if (cfg.yearbase < 70)  %they used a 2-digit year
   fixed_yearbase = cfg.yearbase + 1900;
   warning('you should use a 4-digit yearbase for clarity')
end
