%======================================================================

%                    R D I T Y P E . M 

%                    doc: Wed Jan  7 16:49:09 2009

%                    dlm: Wed Jan  7 16:49:22 2009

%                    (c) 2009 A.M. Thurnherr

%                    uE-Info: 287 0 NIL 0 0 72 0 2 8 NIL ofnI

%======================================================================

% CHANGES BY ANT:

%   Jan  7, 2009: - tightened use of exist()

function [data,d] = rditype(file, verbose, select);

% function data = rditype(file, verbose, select);

%	Read the fixed leader data from a raw ADCP

%	data file.

%	Returns the contents of the fixed leader

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

%	the same length as the number of fields in the

%	record, currently 32.

%	Set verbose=1 for a text output.

% get fixed leader data

 if nargin<3, select = ones(1,31); end

 if nargin<2, verbose = 0; end

 if nargin<1, disp(' provide file name '), return, end

 fid = fopen(file,'r','l');

 hid = 127;  % header identification byte

 [id,n] = fread(fid,2,'uint8');

 status = fseek(fid,0,'bof');

 if id(1) ~= hid

  % Fix parameter for NB-ADCP Kiel

  data=zeros(1,32);

  d.nbytes=NaN;

  d.ntypes=NaN;

  d.Firm_Version=0;

  d.Frequency=150;

  d.Up=0;

  d.Beam_angle=20;

  d.Cell_length=1600;

  d.Pings_per_Ensemble=8;

  d.Pulse_length=1600;

  d.Blank=NaN;

  d.Mode=NaN;

  d.Time_Pings=NaN;

  d.Coordinates=3;

  d.NarrowBand=1;

  return

 end

 [d.nbytes, d.ntypes, offsets] = rdhead(fid, verbose);

 % ---- now at start of fixed leader data

% originally 

% Written by Marinna Martini

% for the U.S. Geological Survey

% Atlantic Marine Geology, Woods Hole, MA

% 1/7/95

data=zeros(1,32);

fld=1;  

% make sure we're looking at the beginning of

% the fixed leader record by testing for it's ID

data(fld)=fread(fid,1,'ushort');

if(data(fld)~=0),

	disp('Fixed Leader ID not found');

	data=[];

	return;

end

fld=fld+1;

% version number of CPU firmware

data(fld)=fread(fid,1,'uchar');

fld=fld+1;

% revision number of CPU firmware

data(fld)=fread(fid,1,'uchar');

d.Firm_Version=data(fld-1)+data(fld)/100;

if verbose, disp(sprintf('CPU Version %d.%d',data(fld-1),data(fld))); end;

fld=fld+1;

% configuration, uninterpreted

data(fld)=fread(fid,1,'uchar');

	b=dec2binv(data(fld));

	freqs=[75 150 300 600 1200 2400];

	junk=bin2decv(b(6:8));

d.Frequency=freqs(junk+1);

d.Up=str2num(b(1));

if verbose, 

	disp(sprintf('Hardware Configuration for LSB %d',data(fld))); 

	disp(sprintf('	System Frequency = %d kHz',freqs(junk+1))); 

	if b(5) == '0', disp('	Concave Beam'); end		

	if b(5) == '1', disp('	Convex Beam'); end		

	junk=bin2decv(b(3:4));

	disp(sprintf('Sensor Configuration #%d',junk+1)); 

	if b(2) == '0', disp('	Transducer head not attached'); end		

	if b(2) == '1', disp('	Transducer head attached'); end		

	if b(1) == '0', disp('	Downward facing beam orientation'); end		

	if b(1) == '1', disp('	Upward facing beam orientation'); end		

end;

fld=fld+1;

data(fld)=fread(fid,1,'uchar');

	b=dec2binv(data(fld));

	angles = [15 20 30 0];

	junk=bin2decv(b(7:8));

d.Beam_angle=angles(junk+1);

if verbose, 

	disp(sprintf('Hardware Configuration MSB %d',data(fld))); 

	disp(sprintf('	Beam angle = %d degrees',angles(junk+1))); 

	junk=bin2decv(b(1:4));

	if junk == 4, disp('	4-beam janus configuration'); end

	if junk == 5, disp('	5-beam janus configuration, 3 demodulators'); end

	if junk == 15, disp('	4-beam janus configuration, 2 demodulators'); end

end;

fld=fld+1;

% real (0) or simulated (1) data flag

data(fld)=fread(fid,1,'uchar');	fld=fld+1;

% undefined

data(fld)=fread(fid,1,'uchar');	fld=fld+1;

% number of beams

data(fld)=fread(fid,1,'uchar');	fld=fld+1;

% number of depth cells

data(fld)=fread(fid,1,'uchar');

d.Depth_Cells=data(fld);

if verbose, disp(sprintf('Number of depth cells %d',data(fld))); end;

fld=fld+1;

% pings per ensemble

data(fld)=fread(fid,1,'ushort');

d.Pings_per_Ensemble=data(fld);

if verbose, disp(sprintf('Pings per ensemble %d',data(fld))); end;

fld=fld+1;

% depth cell length in cm

data(fld)=fread(fid,1,'ushort');

d.Cell_length=data(fld);

if verbose, disp(sprintf('Depth cell size %d cm',data(fld))); end

fld=fld+1;

% blanking distance (WF)

data(fld)=fread(fid,1,'ushort');

d.Blank=data(fld);

if verbose, disp(sprintf('Blank after xmit distance %d cm',data(fld))); end

fld=fld+1;

% Profiling mode (WM)

data(fld)=fread(fid,1,'uchar');

d.Mode=data(fld);

if verbose, disp(sprintf('Profiling mode %d',data(fld))); end

fld=fld+1;

% Minimum correlation threshold (WC)

data(fld)=fread(fid,1,'uchar');

d.Min_Correlation=data(fld);

if verbose, disp(sprintf('Correlation threshold %d',data(fld))); end

fld=fld+1;

% number of code repetitions

data(fld)=fread(fid,1,'uchar');

d.Code_rep=data(fld);

fld=fld+1;

% Minimum percent good to output data (WG)

data(fld)=fread(fid,1,'uchar');

d.Min_Percgood=data(fld);

fld=fld+1;

% Error velocity threshold (WE)

d.Max_Errorvel=data(fld);

data(fld)=fread(fid,1,'ushort');

if verbose, disp(sprintf('Error Velocity Threshold %d mm/s',data(fld))); end

fld=fld+1;

% time between ping groups (TP)

data(fld)=fread(fid,1,'uchar');

d.Time_Pings=data(fld)*60;

fld=fld+1;

data(fld)=fread(fid,1,'uchar');

d.Time_Pings=d.Time_Pings+data(fld);

fld=fld+1;

data(fld)=fread(fid,1,'uchar');

d.Time_Pings=d.Time_Pings+data(fld)/100;

if verbose, disp(sprintf('Time between ping groups %d:%d.%d',...

	data(fld-2), data(fld-1), data(fld))); end

fld=fld+1;

% coordinate transformation (EX)

data(fld)=fread(fid,1,'uchar');

	b=dec2binv(data(fld));

	junk=bin2decv(b(4:5));

d.Coordinates=junk;

d.use_tilt=str2num(b(6));

if verbose, 

	disp(sprintf('Coordinate Transformation = %d',data(fld))); 

	if junk == 0, disp('	Data stored coordinates = Beam'); end

	if junk == 1, disp('	Data stored coordinates = Instrument'); end

	if junk == 2, disp('	Data stored coordinates = Ship'); end

	if junk == 3, disp('	Data stored coordinates = Earth'); end

	if b(6) == '1', disp('	Tilts used in transformation'); end

%	if b(7) == '1', disp('	3-beam solution used, this ensemble'); end

end

fld=fld+1;

% Heading Alignment (EA)

data(fld)=fread(fid,1,'uint16');

d.headin_alignment=data(fld);

fld=fld+1;

% Heading Bias (EB)

data(fld)=fread(fid,1,'uint16');

d.headin_bias=data(fld);

if verbose, disp(sprintf('Heading Bias: %d deg',data(fld)./100)); end

fld=fld+1;

% Sensor source (EZ)

data(fld)=fread(fid,1,'uchar');

	b=dec2binv(data(fld));

d.sensor_source=b;

if verbose,

	disp(sprintf('Sensor Source = %d',data(fld))); 

	if b(2) == '1', disp('	Sound speed computed from ED, ES, ET'); end

	if b(3) == '1', disp('	ED taken from depth sensor'); end	

	if b(4) == '1', disp('	EH taken from xducer heading sensor'); end	

	if b(5) == '1', disp('	EP taken from xducer pitch sensor'); end	

	if b(6) == '1', disp('	ER taken from xducer roll sensor'); end	

	if b(7) == '1', disp('	ES derived from conductivity sensor'); end	

	if b(8) == '1', disp('	ET taken from temperature sensor'); end	

end

fld=fld+1;

% Sensors available

data(fld)=fread(fid,1,'uchar');

	b=dec2binv(data(fld));

d.sensor_avail=b;

if verbose,

	disp(sprintf('Sensor Availability = %d',data(fld))); 

	if b(3) == '1', disp('	depth sensor'); end	

	if b(4) == '1', disp('	heading sensor'); end	

	if b(5) == '1', disp('	pitch sensor'); end	

	if b(6) == '1', disp('	roll sensor'); end	

	if b(7) == '1', disp('	conductivity sensor'); end	

	if b(8) == '1', disp('	temperature sensor'); end	

end

fld=fld+1;

% Bin 1 distance

data(fld)=fread(fid,1,'ushort');

if verbose, disp(sprintf('Distance to center of bin 1: %d cm',data(fld))); end

fld=fld+1;

% xmit pulse length

data(fld)=fread(fid,1,'ushort');

d.Pulse_length=data(fld);

fld=fld+1;

% starting depth cell

data(fld)=fread(fid,1,'uchar');

fld=fld+1;

% ending depth cell

data(fld)=fread(fid,1,'uchar');

fld=fld+1;

% false target reject threshold

data(fld)=fread(fid,1,'uchar');

d.target_max=data(fld);

fld=fld+1;

% spare

data(fld)=fread(fid,1,'uchar');

fld=fld+1;

% transmit lag distance

data(fld)=fread(fid,1,'ushort');

d.xmit_lag=data(fld);

if length(select) == length(data),

	data(find(select==0))=[];

end

%--------------------------------------

function [nb, nt, off] = rdhead(fid, verbose);

%function [nb, nt, off] = rdhead(fid, verbose);

%	Read the header data from a raw ADCP

%	data file opened for binary reading.

%	fid = file handle returned by fopen

%	nb = number of bytes in the ensemble

%	nt = number of data types

%	off = offset to the data for each type

%	Set verbose = 1 for a text output.

% Written by Marinna Martini

% for the U.S. Geological Survey

% Atlantic Marine Geology, Woods Hole, MA

% 1/7/95

data = zeros(1,2);

fld=1;

if ~exist('verbose','var'),

	verbose = 0;

end

nb=[]; nt=[]; off=[];

% make sure we're looking at the beginning of

% the header record by testing for it's ID

junk=fread(fid,2,'uint8');

if((length(junk)~=2) | (ftell(fid)<0)),

	disp('End of file found in rdhead.');

	return;

end

if ((junk(1)~=127) | (junk(2)~=127)),

	disp('Header ID not found');

	return;

end

% get the number of bytes this ensemble

nb = fread(fid,1,'uint16');

if verbose, disp(sprintf('Number of bytes per ensemble %d',nb)); end;

% get the number of data types

fseek(fid,1,'cof');	% skip spare byte position

nt=fread(fid,1,'uchar');

if verbose, disp(sprintf('Number of data types %d',nt)); end;

% get the type offset

off=zeros(nt,1);

for j=1:nt, off(j)=fread(fid,1,'uint16'); end

% ---------------------------------

function d = bin2decv(h)

%BIN2DEC  BIN2DEC('X') returns the decimal number corresponding to the

%        binary number in quotes.  

%		 For example, BIN2DEC([1 1 0 0]) returns 12.

n=length(h);

h=fliplr(h);

d=0;

for i=1:n

	if isstr(h(i)),

		p(i) = str2num(h(i)).*(2^(i-1));

	else

		p(i) = h(i).*(2^(i-1));

	end

end

d=sum(p);

% ---------------------------------

function h = dec2binv(d)

%DEC2BIN DEC2BIN(d) returns the binary number corresponding to the decimal

%        number d.  For example, DEC2BIN(202) returns '11001010'.

%

%	

h=dec2bin(d);

bits=length(h)/8;

if bits ~= fix(bits),

	% not an even multiple of 8, so pad

	npad = 8-((bits-fix(bits))*8);

    h0(1:npad)='0';

	h = [h0,h];

end