function data = rdbtadcp(fid, verbose, select)
% rdbtadcp.m reads the bottom track data from a raw ADCP
%
% function data = rdbtadcp(fid, verbose, select);
%	Read the bottom track data from a raw ADCP
%	data file opened for binary reading
%	Returns the contents of the variable leader
%	as elements of the vector 'data' or an
%	empty matrix if the fixed leader ID is not
%	identified (error condition).
%	If the variable select is provided as a vector
%	of zeros and ones, the function will return
%	only the elements of data which correspond to
%	a one in the vector select.  Select must be the
%  same length as the number of fields in the
%	record, currently 59 (see NFIELDS).
%	Set verbose=1 for a text output.

% Updated 20-sep-07 to incorpoate with adcp_toolbox
% Written by Marinna Martini
% for the U.S. Geological Survey
% Woods Hole Field Center, Woods Hole, MA
% 3/7/01

%%% START USGS BOILERPLATE -------------%
% Use of this program is described in:
%
% Acoustic Doppler Current Profiler Data Processing System Manual 
% Jessica M. Côté, Frances A. Hotchkiss, Marinna Martini, Charles R. Denham
% Revisions by: Andrée L. Ramsey, Stephen Ruane
% U.S. Geological Survey Open File Report 00-458 
% Check for later versions of this Open-File, it is a living document.
%
% Program written in Matlab v7.1.0 SP3
% Program updated in Matlab 7.2.0.232 (R2006a)
% Program ran on PC with Windows XP Professional OS.
%
% "Although this program has been used by the USGS, no warranty, 
% expressed or implied, is made by the USGS or the United States 
% Government as to the accuracy and functioning of the program 
% and related program material nor shall the fact of distribution 
% constitute any such warranty, and no responsibility is assumed 
% by the USGS in connection therewith."
%
%%% END USGS BOILERPLATE --------------

NFIELDS = 59;				% as of RDI Workhorse  Bottom Tracking Addendum Sep 1998
data=zeros(1,NFIELDS);
fld=1;  

if exist('verbose','var') ~= 1,
	verbose = 0;
end
if exist('select','var') ~= 1,
	select = [];
end

% make sure we're looking at the beginning of
% byte 1,2 the bottom track record by testing for it's ID
data(fld)=fread(fid,1,'int16');
if(data(fld)~=1536),
	disp('Bottom Track ID not found');
	data=[];
	return;
end
fld=fld+1;
% byte 3,4 BP/BT pings per ensemble
data(fld)=fread(fid,1,'uint16');
fld=fld+1;
% byte 5,6 BD/BT delay before reaquire
data(fld)=fread(fid,1,'uint16');
fld=fld+1;
% byte 7 BC/BT Minimum Correlation Magnitude
data(fld)=fread(fid,1,'uchar');
fld=fld+1;
% byte 8 BA/BT minimum evaluation amplitude
data(fld)=fread(fid,1,'uchar');
fld=fld+1;
% byte 9 BG/BT minimum percent good BT pings to accept vel data
data(fld)=fread(fid,1,'uchar');
fld=fld+1;
% byte 10 BM/BT mode
data(fld)=fread(fid,1,'uchar');
fld=fld+1;
% byte 11,12 BE/BT maximum error velocity
data(fld)=fread(fid,1,'uint16');
fld=fld+1;
% byte 13,14 spare
data(fld)=fread(fid,1,'uint16');
fld=fld+1;
% byte 15, 16 spare
data(fld)=fread(fid,1,'uint16');
fld=fld+1;
% byte 17-24 BT Range for each beam
data(fld)=fread(fid,1,'uint16');
fld=fld+1;
data(fld)=fread(fid,1,'uint16');
fld=fld+1;
data(fld)=fread(fid,1,'uint16');
fld=fld+1;
data(fld)=fread(fid,1,'uint16');
fld=fld+1;
% byte 25-32 BT Velocity for each beam
data(fld)=fread(fid,1,'int16');
fld=fld+1;
data(fld)=fread(fid,1,'int16');
fld=fld+1;
data(fld)=fread(fid,1,'int16');
fld=fld+1;
data(fld)=fread(fid,1,'int16');
fld=fld+1;
% byte 33-36 BT Correlation (BTCM) for each beam
data(fld)=fread(fid,1,'uchar');
fld=fld+1;
data(fld)=fread(fid,1,'uchar');
fld=fld+1;
data(fld)=fread(fid,1,'uchar');
fld=fld+1;
data(fld)=fread(fid,1,'uchar');
fld=fld+1;
% BT Evaluation Amplitude (BTEA) for each beam
data(fld)=fread(fid,1,'uchar');
fld=fld+1;
data(fld)=fread(fid,1,'uchar');
fld=fld+1;
data(fld)=fread(fid,1,'uchar');
fld=fld+1;
data(fld)=fread(fid,1,'uchar');
fld=fld+1;
% BT Percent Good (BTPG) for each beam, bytes 41-44
data(fld)=fread(fid,1,'uchar');
fld=fld+1;
data(fld)=fread(fid,1,'uchar');
fld=fld+1;
data(fld)=fread(fid,1,'uchar');
fld=fld+1;
data(fld)=fread(fid,1,'uchar');
fld=fld+1;
% Ref Layer Min, Near, Far boundaries
data(fld)=fread(fid,1,'uint16');
fld=fld+1;
data(fld)=fread(fid,1,'uint16');
fld=fld+1;
data(fld)=fread(fid,1,'uint16');
fld=fld+1;
% Reference Layer Velocities, per beam
data(fld)=fread(fid,1,'uint16');
fld=fld+1;
data(fld)=fread(fid,1,'uint16');
fld=fld+1;
data(fld)=fread(fid,1,'uint16');
fld=fld+1;
data(fld)=fread(fid,1,'uint16');
fld=fld+1;
% Reference Layer Correlation (RLCM), per beam
data(fld)=fread(fid,1,'uchar');
fld=fld+1;
data(fld)=fread(fid,1,'uchar');
fld=fld+1;
data(fld)=fread(fid,1,'uchar');
fld=fld+1;
data(fld)=fread(fid,1,'uchar');
fld=fld+1;
% Reference Layer Echo Intensity (RLEI), per beam
data(fld)=fread(fid,1,'uchar');
fld=fld+1;
data(fld)=fread(fid,1,'uchar');
fld=fld+1;
data(fld)=fread(fid,1,'uchar');
fld=fld+1;
data(fld)=fread(fid,1,'uchar');
fld=fld+1;
% Reference Layer Percent Good (RLPG), per beam
data(fld)=fread(fid,1,'uchar');
fld=fld+1;
data(fld)=fread(fid,1,'uchar');
fld=fld+1;
data(fld)=fread(fid,1,'uchar');
fld=fld+1;
data(fld)=fread(fid,1,'uchar');
fld=fld+1;
% Maximum tracking depth 
data(fld)=fread(fid,1,'uint16');
fld=fld+1;
% Receiver Signal Strength (RSSI), per beam
data(fld)=fread(fid,1,'uchar');
fld=fld+1;
data(fld)=fread(fid,1,'uchar');
fld=fld+1;
data(fld)=fread(fid,1,'uchar');
fld=fld+1;
data(fld)=fread(fid,1,'uchar');
fld=fld+1;
% Gain level
data(fld)=fread(fid,1,'uchar');
fld=fld+1;
% BT range MSB
data(fld)=fread(fid,1,'uchar');
fld=fld+1;
data(fld)=fread(fid,1,'uchar');
fld=fld+1;
data(fld)=fread(fid,1,'uchar');
fld=fld+1;
data(fld)=fread(fid,1,'uchar');
%fld=fld+1;

if length(select) == length(data),
	data(find(select==0))=[];
end