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'}
% -- or -- if there is only one variable, it can be a string
% 'vars', 'all'
%
% (unlike instclass "os" you'll never find navigation here)
% 
% 'yearbase', [yearbase]             % REQUIRED
% 'range',  [ens_start ens_stop]     % range of indices to use
% 'step', [step]                     % use only every "step" ensemble
% 'ens_list', [ensemble index list]  % explicit ensemble list overrides step
% 'ends', [0 or 1]                   % 0 (default)
%                                    % 1: get only first and last pair
%                                    % use of 'ends' overrides range, step,
%                                    %    and ens_list
%
% "ens_list" is a misnomer: it refers to vmdas ".ENS" suffix, not "ensemble"
%      this is a raw ping reader.
%

% 00/01/31 EF
% 01/01/26 JH - add instclass, work arguments to match the os template

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


global out yearbase

start = -1;          
stop = 1000000000;



%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.range       = [];
cfg.step        = 1;
cfg.ens_list    = [];
cfg.ends        = 0; %get just the first and last pair of ens
vars            = {'vel', 'sw', 'amp', 'pg', 'ts', 'bt', 'status' 'all'};


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


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

ens_list = cfg.ens_list ;
step     = cfg.step     ; 
if (cfg.ends == 1), out.ends = 1; else, out.ends = 0; 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


% which ensembles?
%try range
if ~(isempty(cfg.range))
   start = cfg.range(1);
   stop = cfg.range(2);
end
% ens_range and out.ends will get used later



fid = fopen(file, 'r', 'ieee-be');  % the OS processor writes little-endian
[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]);

if out.ends & n_profiles > 4,
   ens_list = [first_profile, (first_profile + 1), (last_profile - 1), last_profi
le];
elseif isempty(ens_list)
   ens_list = start:step:stop;
end


NP = floor(n_profiles/step);
NC = config.ncells;
NB = 4; % beams

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;

ipp = 1;
for ip = 1:n_profiles
   [config, buf] = read_buf(fid, config);
   if ischar(buf),
      d.error = buf;
      break
   end
   if rem(ip-1, step) == 0,   %%very inefficient
      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.temperature(ipp) = buf.temperature;
         d.HiV(ipp) = buf.HiV;
         d.XmtAmps(ipp) = buf.XmtAmps;
      end
      if out.vel & config.n.vel, d.vel(:, :, ipp) = buf.vel.'; end
      if out.sw & config.n.sw, d.sw(:, :, ipp) = buf.sw.'; end
      if out.amp & config.n.amp, d.amp(:, :, ipp) = buf.amp.'; end
      if out.pg & config.n.pg, d.pg(:, :, ipp) = buf.pg.'; end
      if out.status & config.n.status, d.status(:, :, ipp) = buf.status.'; end
      if out.bt
         d.bt.vel(:, ipp) = buf.bt.vel;
         d.bt.range(:, ipp) = buf.bt.range;
         d.bt.pg(:, ipp) = buf.bt.pg;
      end
      ipp = ipp + 1;
   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
    %% maybe need to change this depending on beam coords or geographic coords
    %% should be OK for shipboard use
   d.bt.vel = config.vel_multiplier * to_signed(d.bt.vel); 
   d.bt.pg = 6.67 * d.bt.pg;
end


fclose(fid)
d.error = [];


%=========================================================================
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 = [2000]; end  % avoids warning upon early error return
                 %%%% temporary year hack

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 == 1,
   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(config.yearbase, rec(ofs + (1:5))).';

buf.dday = to_day( config.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)

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)
s = u;
ineg = u > 32767;
s(ineg) = s(ineg) - 65536;
ibad = s == -32768;
s(ibad) = NaN;


%----------------------------------------------
function td = decode_time(yearbase, 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 = [yearbase; 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;

