function [ABS_Info,ABS_Data]=read_aquatec_v4(fname);
%
%  [ABS_Info, ABS_Data]=READ_AQUATEC_V4(fname)
%
%  Function to read ABS files from the Aquatec New Version (v3) system
%
%  Input:    
%           fname = filename of abs data including extention
%
%  Output:
%           ABS_Data   : a structure w/ the following fields
%               sdatenum   : Burst starting time (in Matlab serial number, 
%                            use datevec to make yyyy,mo,dd,hh,min,ss). 
%               freq       : Profile recording frequency
%               mdata      : Amplitude data(3-D array) mdata[1:maxbin, time , freq]
%               r          : Matrix of distances from the ABS head in cm (one for
%                            each transducer as these might have been set to be diffferent)
%           ABS_Info   : a structure with ABS settings and system information
%    
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Use of this program is self described.
% Program written in Matlab v7.1.0 SP3
% 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."
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
% Written by George Voulgaris, CPSD Lab., Univ. of South Carolina, 01/30/04
% if you find a bug please email me at: gvoulgaris@geol.sc.edu
% 
% Modified from a function written by Oleg Mouraenko, UF, 10/16/2003
% that was based on an original version created by Peter Traykovski, WHOI
%
% Updates V.1
% 
% 02/04/04 by GV to include all Transducer Serial Numbers
% 02/11/04 by GV to export (i) the number of profiles read and other parameters
%                          (ii) the number of abs freqs/units (i.e., transducers)
%                          (iii) Pulse interval and Pulse frequency
%
% Updates V.3
% 2/27/04 Included the method by Jingping Xu, U.S. Geological Survey, Menlo Park, CA
%         for fast reading (assumes that no other blocks exist after case 22)
% 
%
% Updates V.4
% 8/22/05 by C. Sullivan to suppress display of File Type and Version, and
% to output mdata, sdatenum, freq, and z as the structure ABS_data
% 3/21/06 by C. Sullivan to output 'z' (distance from the ABS) in
% millimeters for ABSTOOLS.

finfo=dir(fname);
flength=finfo.bytes;
%
fid = fopen(fname,'r');
%
No_of_Case_22=0;
%
while ftell(fid)< flength
  try
      command=freadchk(fid,1,'uint8');
      %
      %
     switch command;
         %-----------------------------------------------------------------
         % command=19 - File Type
         %-----------------------------------------------------------------
      case 19;
         junk=freadchk(fid,1,'uint8');
         rec_size=freadchk(fid,1,'uint16');
         chksum=freadchk(fid,1,'uint16');
         sdata=freadchk(fid,rec_size/2,'uint16');
         % Commented out by C. Sullivan
         %disp(['File Type (0=ABS): ',num2str(sdata(1))])
         %disp(['Version          : ',num2str(sdata(2)),'.',num2str(sdata(3))])
         %-----------------------------------------------------------------
         % command=20 - ABS Header
         %-----------------------------------------------------------------
      case 20; 
         junk=freadchk(fid,1,'uint8');
         rec_size=freadchk(fid,1,'uint16');
         chksum=freadchk(fid,1,'uint16');
         %
         InstrumentType=freadchk(fid,1,'uint16');
         InstrumentUnitNo=freadchk(fid,1,'uint16');
         burstnum=freadchk(fid,1,'uint32');
         %
         sdata=freadchk(fid,8,'uint8');
         second=str2num(dec2hex(sdata(2)));
         minute=str2num(dec2hex(sdata(3)));
         hour=str2num(dec2hex(sdata(4)));
         month=str2num(dec2hex(sdata(5)));
         day=str2num(dec2hex(sdata(6)));
         year=str2num(dec2hex(sdata(8)))*100+str2num(dec2hex(sdata(7)));
         sdatenum=datenum(year,month,day,hour,minute,second);
         %
         AuxChannelsBurstInterval=freadchk(fid,1,'uint16');
         reserved=freadchk(fid,1,'uint16');
         wakeupsource=freadchk(fid,1,'uint16');
         %-----------------------------------------------------------------
         % command=21 - ABS Parameters
         %-----------------------------------------------------------------
      case 21;
         junk=freadchk(fid,1,'uint8');
         rec_size=freadchk(fid,1,'uint16');
         chksum=freadchk(fid,1,'uint16');
         %
         BinaryFlag=freadchk(fid,1,'uint8');
         junk=freadchk(fid,1,'uint8');
         %
         ProfileRate=freadchk(fid,1,'uint16');  % Pinging rate in Hz
         Averaging1=freadchk(fid,1,'uint8');    % No of profiles averaged
         freq=ProfileRate/Averaging1;           % Sampling frequency in Hz
         Averaging2=freadchk(fid,1,'uint8');
         Averaging2=2^Averaging2;
         %
         ControlFlagBits=freadchk(fid,1,'uint16');
         futureuse=freadchk(fid,3,'uint16');
         Number_Channels=freadchk(fid,1,'uint16');  % Number of Channels Used
         %  
           for i=1:Number_Channels
               PulseLength(i)=freadchk(fid,1,'uint16');
               Decimation(i)=freadchk(fid,1,'uint16');
               StartDistance_cm(i)=freadchk(fid,1,'uint16');
               StopDistance_cm(i)=freadchk(fid,1,'uint16');
               Pulse_Int(i)=freadchk(fid,1,'uint16');    % 02/11/04 Changed into Pulse_Interval
               Tranducer_Angle(i)=freadchk(fid,1,'int16');
               Scaled_Gain(i)=freadchk(fid,1,'uint8');
               Gain_step_cm(i)=freadchk(fid,1,'uint8');  % 02/13/04 This is not used anymore (per Aquatec)
               Fixed_gain_db(i)=freadchk(fid,1,'int16');
               Gain_step_db(i)=freadchk(fid,1,'float');  % units: db/cm   
               reserved=freadchk(fid,2,'uint16');
               %
               BinSize(i)=1.25*2^(Decimation(i)); % in mm
               numbins(i)=(StopDistance_cm(i)-StartDistance_cm(i))/(BinSize(i)*0.1);
               %
           end
         %  Frequencies used
           tmp=fliplr(dec2bin(65536+BinaryFlag))';
           freqs_used=str2num(tmp([1:Number_Channels]'))'; 
           numfreqs=sum(freqs_used); % number of frequencies 
        %-----------------------------------------------------------------
        % command=40 - ABS Transducer Information
        %-----------------------------------------------------------------
      case 40; % command =40 - ABS Transducer Information
         junk=freadchk(fid,1,'uint8');
         rec_size=freadchk(fid,1,'uint16');
         chksum=freadchk(fid,1,'uint16');
         No_of_Channels=freadchk(fid,1,'uint16');
           for i=1:No_of_Channels
             Transducer_Type(i)=freadchk(fid,1,'uint16');
             Shape_Constant(i)=freadchk(fid,1,'uint16');
             Transducer_diameter(i)=freadchk(fid,1,'uint16');
             Transducer_horizontal_size(i)=freadchk(fid,1,'uint16');
             Transmit_Frequency(i)=freadchk(fid,1,'float');
             Receive_Frequency(i)=freadchk(fid,1,'float');
             TsN=char(freadchk(fid,16,'uchar'))';
             Transducer_serial_No(i,:)=TsN(1:12);      % Changed 2/4/4, GV
           end
        %-----------------------------------------------------------------
        % command=22 - ABS Data 16 - FAST Method by Jingping Xu
        %              It assumes no other data exist after this block
        %-----------------------------------------------------------------
      case 22; 
         No_of_Case_22 = No_of_Case_22 + 1;
         junk=freadchk(fid,1,'uint8');
         rec_size=freadchk(fid,1,'uint16');
         chksum=freadchk(fid,1,'uint16');
         [amp, count]=fread(fid,inf,'uint16');
         mm=sum(numbins(1:numfreqs));
         len=mm+3;
         nn=floor(count/len);
         amp=amp(1:nn*len);
         amp=reshape(amp,[len,nn]);
         amp=amp(1:mm,:);
         for i=1:numfreqs 
            mdata(:,:,i)=amp(1:mm/numfreqs,:);
            amp(1:mm/numfreqs,:)=[];
         end
      break
        %-----------------------------------------------------------------
        % command=45 - Details
        %-----------------------------------------------------------------
     case 45 %details
        fread(fid,3,'uint16');
        fread(fid,3,'uint32');
        fread(fid,14,'char');
        fread(fid,2,'ubit16');
        fread(fid,3,'uint16');
        bn=fread(fid,1,'uint16');
       for i=1:bn
         fread(fid,2,'uint16');
         fread(fid,12,'char');
       end
        %-----------------------------------------------------------------
        % command=29 - Errors Information
        %-----------------------------------------------------------------
     case 29 %errors info
        fread(fid,2,'int16');
        %-----------------------------------------------------------------
        % command=41 - Complex Data
        %-----------------------------------------------------------------
     case 41 %complex data
        for i=1:numfreqs;
           fread(fid,[2,channelsample(i)],'int16');
        end;
        %-----------------------------------------------------------------
        % command=44
        %-----------------------------------------------------------------
     case 44
        fread(fid,3,'uint16');
        %-----------------------------------------------------------------
        % command=46 
        %-----------------------------------------------------------------
     case 46
        sdata=fread(fid,4,'uint16');
        if isequal(sdata(2),10)
          fread(fid,7,'uint16');
          fread(fid,1,'int16');
        end;
        %-----------------------------------------------------------------
        % command=49
        %-----------------------------------------------------------------
     case 49
          sdata=fread(fid,2,'uint16');
          if isequal(sdata(1),1)
            fread(fid,1,'char');
          else
            fread(fid,1,'long')
            fread(fid,3,'double');
          end;
        %-----------------------------------------------------------------
        % other commands no recognized are skipped
        %----------------------------------------------------------------- 
      otherwise
         disp(sprintf('WARNING >> Command "%i" is not supported. Data omitted.',command));
         junk=freadchk(fid,1,'uint8');
         rec_size=freadchk(fid,1,'uint16');
         chksum=freadchk(fid,1,'uint16');
         sdata=freadchk(fid,rec_size,'uint8')';
      end;
   catch
      if No_of_Case_22==0;
         error('ERROR   >> Error detected. No Profile Data (case 22) were found');
      else 
         disp('WARNING >> Error detected. Data may be corrupted');
      end
   end;   
end; 
fclose(fid);
%
[as,lastprof]=size(mdata(:,:,1));
numprofs=lastprof;
%
%   Create the range array for each frequency (these could be
%   set to be different in V.3. The range is set in the middle of
%   the bin.
%
z=zeros(max(numbins),numfreqs)*NaN;
for i=1:numfreqs
    z(1:numbins(i),i)=StartDistance_cm(i) + [1:numbins(i)]' * (BinSize(i)/10) - BinSize(i)/10/2 ;  % Changed GV
    % The range has been estimated in th middle of each bin -> thus the - half bin size. 
    Pulse_Interval(i)=Pulse_Int(i);         % 02/11/04 Changed into Pulse_Interval
    Pulse_Frequency(i)=1000/Pulse_Int(i);   % 02/11/04, GV, new added
end
%
%
ABS_Info = struct('InstrumentType',InstrumentType,...
   'InstrumentUnitNo',InstrumentUnitNo,...    
   'burstnum',burstnum,...
   'year',year,...
   'month',month,...
   'day',day,...
   'hour',hour,...
   'minute',minute,...
   'second',second,...
   'ProfileRate',ProfileRate,...
   'Averaging1',Averaging1,...
   'wakeupsource',wakeupsource,...
   'numprofiles',numprofs,...                 % changed on 02/11/04 GV (was numprofiles)
   'No_of_Freqs',numfreqs,...                 % added on 02/11/04 GV 
   'Pulse_Interval',Pulse_Interval,...        % added on 02/11/04 GV 
   'Pulse_Frequency',Pulse_Frequency,...      % added on 02/11/04 GV 
   'Scaled_Gain',Scaled_Gain(1:numfreqs),...
   'Gain_step_cm',Gain_step_cm(1:numfreqs),...
   'Fixed_gain_db',Fixed_gain_db(1:numfreqs),...
   'Gain_step_db',Gain_step_db(1:numfreqs),...
   'Transducer_Type',Transducer_Type,...
   'Shape_Constant',Shape_Constant(1:numfreqs),...
   'Transducer_diameter_mm',Transducer_diameter(1:numfreqs),...
   'Transducer_horizontal_size_mm',Transducer_horizontal_size(1:numfreqs),...
   'Transmit_Frequency_MHz',Transmit_Frequency(1:numfreqs),...  % Changed from HZ to MHz 2/4/4, by GV
   'Receive_Frequency_MHz',Receive_Frequency(1:numfreqs),...    % Changed from Hz to MHz 2/4/4, by GV
   'Transducer_serial_No',Transducer_serial_No(1:numfreqs,:));  % Changed 2/4/4, by GV
%
%
ABS_Data = struct('mdata',mdata,...
    'freq',freq,...
    'sdatenum',sdatenum,...
    'r',z*10);
%
%
function M=freadchk(fid,count,type);
% Modification of FREAD function to check if all
% needed data was read.

M = fread(fid,count,type);
if length(M(:))<count
   disp('WARNING >> Unexpected End of File');
   error('End of file');
end;
return;
