function varargout = adp(varargin)
% adp.m: a general function for extracting variables from Sontek
%        adp data files.  This version is specific to the 4-beam "type 2"
%        files, but could easily be modified for "type 1" or others.
% Usage:
%    [ts, ts2, v, s, a, n, dd] = adp(fileroot, 'ts', 'ts2', ...
%                    'vel', 'std', 'amp', 'noise', step, range);
%    Inputs:
%        fileroot: name of adp data file, with or without
%              extension ('.adp' or '.ADP')
%        'ts' ... 'amp' is a list of the desired variables.
%              The output arrays will be in the order of these input
%              strings, and only the ones listed will be output.
%        step: optional scalar for subsampling; 1 (default) means take
%              every sample, 2 means take every other one, etc.
%        range: optional 2-element vector specifying the first and last
%              profile in the file to be extracted.  The profiles
%              taken are range(1):step:range(2); default range is
%              entire file.
%     Dimensions:
%        NP: number of profiles extracted
%        NC: number of depth cells
%        NB: number of beams (hardwired here to 4!)
%     Outputs:
%        (dd is always available; the others are available only if
%        requested via the input arguments.)
%        ts: NPx5 array.  Columns are heading, pitch, roll, temperature,
%              and pressure.  Units are degrees, degrees C, and decibars.
%        ts2: NPx5 array.  Columns are Xflux, Yflux, Zflux,
%              compass temperature, and voltage.  Units are
%              microteslas, degrees C, and Volts.
%        v: NCxNBxNP array of velocities in m/s
%        s: NCxNBxNP array of velocity standard deviations in m/s
%        a: NCxNBxNP array of signal/noise in db (counts * 0.43)
%              based on the noise level given in the header.
%              If "rawamp" is requested, then it is just raw counts.
%        n: 1xNB array of noise levels from header.
%        dd: NPx1 time of each profile in decimal days since the
%              start of the year; noon on January 1 is 0.5, NOT 1.5.
%     Note: the squeeze function has to be used to plot slices of
%     these 3-D arrays.  For example, to plot the time series of
%     velocity in beam 1, cells 2 and 6, use:
%         plot(dd, squeeze(v([2 6],1,:)).')
%
%     Sun  98/03/22 Eric Firing
%     Mon  98/04/27 EF added ts2 array.
%     2000/03/17 EF allow input filename with extension
%     2000/04/19 EF add noise output

% The following are hardwired here, but could be read from
% the file.
nbeams = 4;
fileheadlength = 416;
profheadlength = 80;

% Defaults:
step = 1;
profs = 0;
raw = 0;    % Default: snr rather than raw amplitude.

out.ts = 0;
out.ts2 = 0;
out.vel = 0;
out.amp = 0;
out.std = 0;
out.noise = 0;

% Read the input arguments:
fileroot = varargin{1};
nout = 0;
for iarg = 2:nargin,
   arg = varargin{iarg};
   if ischar(arg),
      switch lower(arg)
         case 'ts',  nout = nout + 1; out.ts = nout;
         case 'ts2', nout = nout + 1; out.ts2 = nout;
         case 'vel', nout = nout + 1; out.vel = nout;
         case 'amp', nout = nout + 1; out.amp = nout;
         case 'rawamp', nout = nout + 1; out.amp = nout; raw = 1;
         case 'noise', nout = nout + 1; out.noise = nout;
         case 'std', nout = nout + 1; out.std = nout;
         otherwise
            error('invalid data type requested')
      end
   else
      nn = length(arg);
      if  nn == 1,
         step = arg;
      elseif nn == 2,
         profs = arg;
      else
         error('invalid size of numeric argument')
      end
   end
end

if exist(fileroot) == 2,
   filename = fileroot;
else
   filename = [fileroot, '.adp'];
   if exist(filename) ~= 2,
      filename = [fileroot, '.ADP'];
   end
end
fid = fopen(filename, 'r', 'l'); % Little-endian data file.
if fid == -1,
   error(['Can''t open input file: ', fileroot, ', with or without .adp or .ADP']);
end

% Read the pressure calibration factors:
fseek(fid, 70, 'bof');
p_params = fread(fid, 2, 'int32');
pscale = 1e-8*p_params(1);
poffset = 1e-5*p_params(2);
fseek(fid, 82, 'bof');
p_params = fread(fid, 1, 'int16');
pscale2 = 1e-12*p_params;

% Read the number of cells:
fseek(fid, 178, 'bof');
ncells = fread(fid, 1, 'int16');

proflength = profheadlength + 4*nbeams*ncells + 2; % 2 for vel, 1 each std, amp
%% The extra two bytes are the checksum.

% How many profiles are in the file?
fseek(fid, 0, 'eof');
file_nbytes = ftell(fid);
nprofiles = fix((file_nbytes - fileheadlength)/proflength);


if size(profs(:),1) ~= 2,
   irange = 1:step:nprofiles;
else
   irange =  max([1 profs(1)]):step:min([nprofiles profs(2)]);
end

np = length(irange);  % Number of profiles to be extracted.
tt = NaN + zeros(np,6);




ofs.time  = 20;
ofs.noise = 68;
ofs.ts    = 40;
ofs.ts2   = 60;
ofs.vel   = profheadlength;
ofs.amp   = profheadlength + 3*nbeams*ncells;
ofs.std   = profheadlength + 2*nbeams*ncells;


if out.ts,   A = zeros(np,5);               end
if out.ts2,  B = zeros(np,5);               end
if out.ts2,  B1 = zeros(np,2);              end
if out.ts2,  B2 = zeros(np,1);              end
if out.ts2,  B3 = zeros(np,2);              end
if out.vel,  v = zeros(ncells, nbeams, np); end
if out.amp,  a = zeros(ncells, nbeams, np); end
if out.std,  s = zeros(ncells, nbeams, np); end

cs_offset = hex2dec('0xa596');
bad_checksum = 0;
iout = 1;
for ii = irange,
   p_offset = fileheadlength+(ii-1)*proflength;

   fs = fseek(fid, p_offset, 'bof');
   buf = fread(fid, proflength-2, 'uint8');
   checksum = fread(fid, 1, 'uint16');
   cs = rem(sum(buf) + cs_offset, 65536);
   if checksum ~= cs,
      bad_checksum = 1;
      break
   end

   if iout == 1,
      fs = fseek(fid, p_offset + ofs.time - 2, 'bof');
      year = fread(fid, 1, 'int16');
   end

   fs = fseek(fid, p_offset + ofs.time, 'bof');
   tt(iout,:) = fread(fid, [1,6], 'int8');

   if out.ts,
      fs = fseek(fid, p_offset + ofs.ts, 'bof');
      A(iout,:) = fread(fid, [1,5], 'int16');
   end

   if out.ts2,
      fs = fseek(fid, p_offset + ofs.ts2, 'bof');
      B1(iout,:) = fread(fid, [1,2], 'int16');
      fs = fseek(fid, 11, 'cof');
      B2(iout,:) = fread(fid, [1,1], 'int8');
      B3(iout,:) = fread(fid, [1,2], 'int16');
   end

   if out.vel,
      fs = fseek(fid, p_offset + ofs.vel, 'bof');
      v(:,:,iout) = fread(fid, [ncells,nbeams], 'int16') * 1e-3;
   end

   if out.std,
      fs = fseek(fid, p_offset + ofs.std, 'bof');
      ss = fread(fid, [ncells,nbeams], 'int16');
      ineg = find(ss<0);
      ss(ineg) = ss(ineg) + 256.0;
      s(:,:,iout) = ss * 1e-3;
   end

   if iout == 1 & (out.noise ~= 0 | raw == 0)
      fs = fseek(fid, p_offset + ofs.noise, 'bof');
      noise = fread(fid, [1,4], 'uint8');
      N = noise(ones(ncells,1),:);
   end

   if out.amp,
      fs = fseek(fid, p_offset + ofs.amp, 'bof');
      amp = fread(fid, [ncells,nbeams], 'uint8');
      if raw == 0,
         a(:,:,iout) = 0.43 * (amp - N);
      else
         a(:,:,iout) = amp;
      end
   end

   iout = iout+1;
end

fclose(fid);

if bad_checksum,
   imiss = ii:np;
   tt(imiss,:) = [];
   if out.ts,   A(imiss, :) = [];   end
   if out.ts2,  B(imiss, :) = [];   end
   if out.ts2,  B(imiss, :) = [];   end
   if out.ts2,  B2(imiss, :) = [];  end
   if out.ts2,  B3(imiss, :) = [];  end
   if out.vel,  v(:, :, imiss) = []; end
   if out.amp,  a(:, :, imiss) = []; end
   if out.std,  s(:, :, imiss) = []; end
end

if out.ts,
   iwrap = find(A(:,5) < 0);
   if ~isempty(iwrap), A(iwrap,5) = A(iwrap,5) + 65536; end
   A(:,5) = poffset + pscale * A(:,5) + pscale2 * A(:,5).^2;
   A(:,1:3) = 0.1*A(:,1:3);
   A(:,4) = 0.01*A(:,4);
end

if out.ts2,
   B(:,1:2) = B1 * 0.01;               % Xflux, Yflux; microteslas
   B(:,3)   = B3(:,1) * 0.01;          % Zflux
   B(:,4)   = B3(:,2) * 0.1;          % CompassTemp
   B(:,5)   = B2 * 0.2;               % Voltage
end

tn = datenum(year,tt(:,2),tt(:,1),tt(:,4),tt(:,3),tt(:,6)+tt(:,5)/100);
dday = tn - datenum(year,1,1,0,0,0);

varargout{nout + 1} = dday;

if out.ts,   varargout{out.ts}  = A; end
if out.ts2,  varargout{out.ts2} = B; end
if out.vel,  varargout{out.vel} = v; end
if out.std,  varargout{out.std} = s; end
if out.amp,  varargout{out.amp} = a; end
if out.noise, varargout{out.noise} = noise; end



