function [FanData, FanHeader, FanSwitches] = readfan(fname, verbose)
% readfan - Read in data from Imagenex fan beam sonar
% [FanData, FanHeader, FanSwitches] = readfan(fname, verbose);
%
% fname = file name
% verbose = 0; turn off output
%
% FanData.imagedata: [ping x bin double], image data with header removed
%
% HeaderData is a dump of the header info sent with *each ping*
% HeaderData.ReturnDataHeaderType = IGX/IMX etc.
% HeaderData.HeadID = ID number
% HeaderData.HeadPosition = head position in counts
% HeaderData.HeadAngle = head rotation angle in degrees
% HeaderData.StepDirection = direction of rotation, -1 = CCW, 1 = CW
% HeaderData.Range = range setting
% HeaderData.ProfileRange = first digitized value above the threshold in cm
% HeaderData.NDataBytes = number of byte in the ping
% HeaderData.NReturnBytes = number of bytes in the ping + header
% HeaderData.FanTime = time of first ping for logger data,
%                   vector of times for 881B file
% HeaderData.StepSize = degrees per step of rotation
% HeaderData.FileName = name of the file read
% HeaderData.NPoints = number of data points returned for the ping
% HeaderData.SectorWidth = size of sector swept in degrees
%
% FanSwitches is a dump of the switches sent to the sonar for this sweep
% FanSwitches =
%            HeadID: 16
%             Range: 5
%          Reserved: [NaN 0 0 0 NaN 0 0]
%       RangeOffset: 0
%           RevHold: 'Reverse Step Direction'
%       MasterSlave: 'Send Data Slave Mode'
%         StartGain: 4
%              LOGF: 20
%        Absorption: 0.0600
%        TrainAngle: 0
%       SectorWidth: 348
%          StepSize: 0.6000
%       PulseLength: 10
%     StepDirection: 0
%      MoveRelative: 0
%        DataPoints: 500
%          DataBits: 8
%            UpBaud: 6
%           Profile: 'OFF'
%         Calibrate: 'Calibrate'
%       SwitchDelay: 2

% TODO force user to provide settings such as
% head orientation for logger files

% 16jul07 fix bad bit parsing of Profile Range and NDataBytes in header
% Written by Marinna Martini 5/22/07
% for the U.S. Geological Survey, Woods Hole Field Center

NRETURNHEADERBYTES = 12;
NWHSCFILEHEADERBYTES = 61;
N81BFILEHEADERBYTES = 100;
FanData = [];
FanHeader = [];
FanSwitches = [];

if ~exist('fname', 'var') || ~exist(fname, 'file'),
    [fname, pathname] = uigetfile( ...
        {'*.F*','WHSC logger files (S1*.F*)'; ...
        '*.81B','Win8812W output files (*.81B)'; ...
        '*.*',  'All Files (*.*)'});
    if isempty(fname), return, end
    fname = fullfile(pathname, fname);
end
if ~exist('verbose', 'var'),
    verbose = 1;
end
fid = fopen(fname,'r');
if fid == -1,
    disp(['error reading file ',fname]);
    return
else
    disp(['Now reading file ', fname]);
end
FanHeader.FileName = fname;

% detect which file we have
% read the first few bytes
tag = fread(fid,3,'uchar');
fseek(fid,0,'bof');

switch(char(tag)')
    case 'The' % this is USGS WHSC logger data
        % extract year info from file
        d = dir(fname);
        datayear = datevec(d.datenum);
        % first part is the file name, as written by the logger,
        % which is essentially the time stamp
        fheader = fread(fid,NWHSCFILEHEADERBYTES,'uchar');
        [FanTime, FanSwitches] = parseLoggerHeader(fheader, datayear(1), verbose);
        clear d

        % now at the data
        % decipher the first 13 bytes of sonar header so we know what we are reading
        sheader = fread(fid,NRETURNHEADERBYTES,'uchar');
        header = parseFanheader(sheader, FanSwitches, verbose); % display them once
        if verbose, disp(sprintf('Number of return bytes per ping = %d',header.NReturnBytes)); end
        npoints = header.NReturnBytes-NRETURNHEADERBYTES+1;
        if verbose, disp(sprintf('Number of points per ping = %d',npoints)); end
        FanHeader.FanTime = FanTime; % pass out the time stamp
        FanHeader.NPoints = FanSwitches.DataPoints;
        FanHeader.StepSize = FanSwitches.StepSize;
        FanHeader.SectorWidth = FanSwitches.SectorWidth;
        clear sheader

        % now figure out how many pings of data are in this file from the file size and header info.
        fileinfo = dir(fname);
        if verbose, disp(sprintf('File size = %d',fileinfo.bytes)); end
        npings = (fileinfo.bytes-NWHSCFILEHEADERBYTES)./header.NReturnBytes;
        if verbose, disp(sprintf('Number of pings in the file = %d',npings)); end

        fseek(fid, NWHSCFILEHEADERBYTES, 'bof'); % go to data start
        raw = fread(fid,'uchar'); % read in all the ping data
        imagedata = reshape(raw(1:floor(npings)*header.NReturnBytes),header.NReturnBytes, floor(npings));
        % let's check out the header for each ping
        for iping = 1:npings,
            header = parseFanheader(imagedata(1:13,iping), FanSwitches, 0);
            % rearrange the output for convenience
            FanHeader.ReturnDataHeaderType{iping} = header.ReturnDataHeaderType;
            FanHeader.HeadID(iping) = header.HeadID;
            FanHeader.HeadPosition(iping) = header.HeadPosition;
            FanHeader.HeadAngle(iping) = header.HeadAngle;
            FanHeader.StepDirection(iping) = header.StepDirection;
            FanHeader.Range(iping) = header.Range;
            FanHeader.ProfileRange(iping) = header.ProfileRange;
            FanHeader.NDataBytes(iping) = header.NDataBytes;
            FanHeader.NReturnBytes(iping) = header.NReturnBytes;
        end

        FanData.imagedata = imagedata(14:end,:); % pass back just data

    case '81B' % this is data from an 81B file
        fheader = fread(fid,N81BFILEHEADERBYTES,'uchar');
        [FanTime, FanSwitches] = parse81BHeader(fheader,1);

        % now at the data
        % decipher the first 13 bytes of sonar header so we know what we are reading
        sheader = fread(fid,NRETURNHEADERBYTES,'uchar');
        header = parseFanheader(sheader, FanSwitches, 1); % display them once
        FanHeader.NPoints = header.NDataBytes;
        FanHeader.StepSize = FanSwitches.StepSize;
        FanHeader.SectorWidth = FanSwitches.SectorWidth;
        clear sheader

        % now figure out how many pings of data are in this file from the file size and header info.
        fileinfo = dir(fname);
        disp(sprintf('The file size is %d',fileinfo.bytes));
        npings = (fileinfo.bytes)./(FanSwitches.TotalBytes);
        disp(sprintf('Number of pings in the file = %d',npings))

        % allocate memory
        fseek(fid, 0, 'bof'); % go to data start

        raw = fread(fid,'uchar'); % read in all the ping data at once
        imagedata = reshape(raw, FanSwitches.TotalBytes, npings);
        FanData.raw = imagedata;

        % let's check out the header for each ping
        for iping = 1:npings,
            [FanTime, switches] = parse81BHeader(imagedata(1:100,iping),0);
            FanHeader.FanTime(iping) = FanTime;
            header = parseFanheader(imagedata(101:112,iping), FanSwitches, 0);
            if ~isempty(header),
                % rearrange the output for convenience
                FanHeader.ReturnDataHeaderType{iping} = header.ReturnDataHeaderType;
                FanHeader.HeadID(iping) = header.HeadID;
                FanHeader.HeadPosition(iping) = header.HeadPosition;
                FanHeader.HeadAngle(iping) = header.HeadAngle;
                FanHeader.StepDirection(iping) = header.StepDirection;
                FanHeader.Range(iping) = header.Range;
                FanHeader.ProfileRange(iping) = header.ProfileRange;
                FanHeader.NDataBytes(iping) = header.NDataBytes;
                if ~isempty(header.NReturnBytes), FanHeader.NReturnBytes(iping) = header.NReturnBytes;
                else FanHeader.NReturnBytes(iping) = NaN;
                    disp(sprintf('Problem at ping #%d',iping))
                end
            else
                disp(sprintf('Problem at ping #%d',iping))
            end
        end

        FanData.imagedata = imagedata(113:end,:); % pass back just data

    case {'IGX','IMX','IQX','IPX'} % this is USGS WHSC logger *.tmp format
        % extract year info from file
        d = dir(fname);
        FanTime = julian(datevec(d.datenum));
        FanHeader.FanTime = FanTime;

        % need to know the data bits, 4, 8 or 16!
        button = questdlg('Select number of Data Bits','No SwitchInformation','4','8','16','8');
        FanSwitches.DataBits = str2num(char(button));

        % this file format has no header, it goes right to the data
        % decipher the first 13 bytes of sonar header so we know what we are reading
        sheader = fread(fid,NRETURNHEADERBYTES,'uchar');
        header = parseFanheader(sheader, FanSwitches, 1); % display them once
        if verbose, disp(sprintf('Number of return bytes per ping = %d',header.NReturnBytes)); end
        npoints = header.NReturnBytes-NRETURNHEADERBYTES+1;
        if verbose, disp(sprintf('Number of points per ping = %d',npoints)); end
        clear sheader

        % now figure out how many pings of data are in this file from the file size and header info.
        fileinfo = dir(fname);
        disp(fileinfo.bytes);
        npings = (fileinfo.bytes)./header.NReturnBytes;
        disp(sprintf('Number of pings in the file = %d',npings))

        fseek(fid, 0, 'bof'); % go to data start
        raw = fread(fid,'uchar'); % read in all the ping data
        imagedata = reshape(raw, header.NReturnBytes, npings);

        % let's check out the header for each ping
        for iping = 1:npings,
            header = parseFanheader(imagedata(1:13,iping), FanSwitches, 0);
            % rearrange the output for convenience
            FanHeader.ReturnDataHeaderType{iping} = header.ReturnDataHeaderType;
            FanHeader.HeadID(iping) = header.HeadID;
            FanHeader.HeadPosition(iping) = header.HeadPosition;
            FanHeader.HeadAngle(iping) = header.HeadAngle;
            FanHeader.StepDirection(iping) = header.StepDirection;
            FanHeader.Range(iping) = header.Range;
            FanHeader.ProfileRange(iping) = header.ProfileRange;
            FanHeader.NDataBytes(iping) = header.NDataBytes;
            FanHeader.NReturnBytes(iping) = header.NReturnBytes;
        end

        FanData.imagedata = imagedata(14:end,:); % pass back just data
        FanHeader.NPoints =  header.NReturnBytes-NRETURNHEADERBYTES+1;
        % TODO must fix this hardwire!
        FanHeader.StepSize = 0.15;

    otherwise
        disp('Data is unrecognized')
end

fclose(fid);

% put a few switch values into the header for convenience
if ~exist('File81BHeader','var'),
    FanHeader.Orientation = 'DOWN'; % assume down on all tripods
    if verbose,
        disp(FanHeader.Orientation);
    end
end

return

% read a WHSC logger file header from a S1*.F* file
function [FanTime, FanSwitches] = parseLoggerHeader(fheader, datayear, verbose)

% bytes 1:33 are:
% The file name is \S1\S1031414.p12
% it is possible that this can vary
% extract the time information from the file header
cfheader = char(fheader);
MM = sscanf(cfheader(24:25),'%d');
dd = sscanf(cfheader(26:27),'%d');
hh = sscanf(cfheader(28:29),'%d');
mm = sscanf(cfheader(32:33),'%d');
FanTime = julian([datayear, MM, dd, hh, mm, 0]);
clear MM dd hh mm cfheader
% next thing the logger writes out are the switches.  29 bytes.
% not sure why byte #34 is written.  It is ASCII LF = x0A
% pull out switch data for indexing convenience
FanSwitches = parseFanswitches(fheader(35:end), verbose);

return

function [FanTime, File81BHeader] = parse81BHeader(fheader, verbose)

File81BHeader.ReturnDataHeaderType = char(fheader(1:3))';
junk = [0 128 250 500 750 1000 1500 2000 3000 4000];
File81BHeader.NDataBytes = junk(fheader(4)+1);
File81BHeader.TotalBytes = fheader(5).*256+fheader(6);
File81BHeader.nToRead = fheader(7).*256+fheader(8);
% date format 11-May-2007.21:30:05
dstr = char(fheader(9:29)');
dstr(12) = ' ';
hsec = str2num(char(fheader(31:32))); % byte 30 is '.' and byte 33 is null
FanTime = julian(datevec(datenum(dstr)));
File81BHeader.xBytes = fheader(35).*128;
File81BHeader.TotalBytes = File81BHeader.TotalBytes + File81BHeader.xBytes;
File81BHeader.StepDirection = bitget(fheader(38),8); % 1 = CW
File81BHeader.Orientation = bitget(fheader(38),7); % 1 = Up
junk = {'Sector','Polar','Sidescan'};
n = bitshift(bitand(fheader(38), 56),-3);
File81BHeader.Mode = junk{n+1};
junk = [0 0.3 0.6 0.9 1.2 2.4];
File81BHeader.StepSize = junk(bitand(fheader(38), 7)+1);
File81BHeader.StartGain = fheader(39);
File81BHeader.SectorSize = fheader(40).*3;
File81BHeader.TrainAngle = (fheader(41).*3) - 210;
File81BHeader.RangeOffset = fheader(42);
File81BHeader.Absorption = fheader(43).*0.1;
File81BHeader.ProfileGrid = bitget(fheader(44),8); % 1 = ON
File81BHeader.Zero = bitget(fheader(44),7); % 1 = Dn
junk = [4 8 14];
n = bitshift(bitand(fheader(44), 56),-3);
File81BHeader.DataBits = junk(n+1);
File81BHeader.LogF = bitand(fheader(44), 7).*10;
File81BHeader.PulseLength = fheader(45);
junk = {'Off','PointsOnly','Lo-mix','Hi-Mix'};
File81BHeader.Profile = junk{fheader(46)+1};
if bitget(fheader(47),8),
    File81BHeader.SoundVelocity = (bitand(fheader(46),127).*256+fheader(47))/10;
else  File81BHeader.SoundVelocity = 1500;
end
File81BHeader.Comment = fheader(49:80);

if verbose,
    disp(File81BHeader)
end

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% read the header that comes with each ping.
% applies to logger and live data
function FanHeader = parseFanheader(sheader, switches, verbose)

% Marinna's method to read the return data header
FanHeader.ReturnDataHeaderType = char(sheader(1:3))';
FanHeader.HeadID = sheader(4); % always 0x10 for 881-000-100 360 degree scanning head
if FanHeader.HeadID ~= 16,
    % expected 0x10 for a 881-000-100 360 degree scanning head
    disp(sprintf('Wrong head ID %d',sheader(4)))
    FanHeader = [];
    return
end
b=dec2base(sheader(5),2,8); % get the bits
if b(7),  FanHeader.SerialStatus = 'Switches Accepted'; end
if b(8),  FanHeader.SerialStatus = 'Character Overrun'; end
% head position
FanHeader.HeadPosition = bitand(63,sheader(7)).*128 + bitand(127,sheader(6));
FanHeader.HeadAngle = ( FanHeader.HeadPosition-1400).*0.15; % degrees
if bitget(sheader(7),6),  FanHeader.StepDirection = 1; %CW
else FanHeader.StepDirection = -1; %CCW
end
FanHeader.Range = sheader(8); % meters
%FanHeader.ProfileRange = bitand(127,sheader(10)).*127 + bitand(127,sheader(9)); % meters
%FanHeader.NDataBytes = bitand(127,sheader(12)).*127 + bitand(127,sheader(11)); % meters
% MM 10jul07
% bitshift(A, k) % negative = shift right
PRHB = bitshift(bitand(sheader(10), 126), -1);
PRLB = bitor(bitshift(bitand(sheader(10), 1), 7),bitand(sheader(9), 127));
FanHeader.ProfileRange = bitor(bitshift(PRHB,8),PRLB);
PRHB = bitshift(bitand(sheader(12), 126), -1);
PRLB = bitor(bitshift(bitand(sheader(12), 1), 7),bitand(sheader(11), 127));
FanHeader.NDataBytes = bitor(bitshift(PRHB,8),PRLB);
% given the header info, figure out how many data point we're dealing with
% the number of data bits is hardwired to 8 in the logger source code (new_bits)
FanHeader.NReturnBytes = [];
FanHeader.NZeroFillBytes = [];
switch(FanHeader.ReturnDataHeaderType(2)),
    case {'M'},   % switch data command for data points = 25
        if switches.DataBits == 4, % 128 Data Byes, 256 Points
            FanHeader.NReturnBytes = 141;
            FanHeader.NPoints = 256;
            FanHeader.NZeroFillBytes = 15;
        elseif switches.DataBits == 8, % 252 Data Byes, 252 Points
            FanHeader.NReturnBytes = 265;
            FanHeader.NPoints = 252;
            FanHeader.NZeroFillBytes = 19;
        elseif switches.DataBits == 14, % 500 Data Byes, 250 Points
            FanHeader.NReturnBytes = 513;
            FanHeader.NPoints = 250;
            FanHeader.NZeroFillBytes = 27;
        end
    case {'G'},   % switch data command for data points = 50
        if switches.DataBits == 4, % 252 Data Byes, 504 Points
            FanHeader.NReturnBytes = 265;
            FanHeader.NPoints = 504;
            FanHeader.NZeroFillBytes = 19;
        elseif switches.DataBits == 8, % 500 Data Byes, 500 Points
            FanHeader.NReturnBytes = 513;
            FanHeader.NPoints = 500;
            FanHeader.NZeroFillBytes = 27;
        elseif switches.DataBits == 14, % 500 Data Byes, 250 Points
            FanHeader.NReturnBytes = 513;
            FanHeader.NPoints = 250;
            FanHeader.NZeroFillBytes = 27;
        end
    case {'Q'},  % switch data command for data points = 200
        switch(FanHeader.Range), % range is from 5 to 255 m in 5 m increments
            case 5,
                FanHeader.NPoints = 500;
                if switches.DataBits == 4, % 250 Data Byes, 500 Points
                    FanHeader.NReturnBytes = 263;
                    FanHeader.NZeroFillBytes = 21;
                elseif switches.DataBits == 8, % 500 Data Byes, 500 Points
                    FanHeader.NReturnBytes = 513;
                    FanHeader.NZeroFillBytes = 27;
                elseif switches.DataBits == 14, % 1000 Data Byes, 500 Points
                    FanHeader.NReturnBytes = 1013;
                    FanHeader.NZeroFillBytes = 39; % there was no explicit case for this
                end
            case 10,
                FanHeader.NPoints = 1000;
                if switches.DataBits == 4, % 500 Data Byes, 1000 Points
                    FanHeader.NReturnBytes = 513;
                    FanHeader.NZeroFillBytes = 27;
                elseif switches.DataBits == 8, % 1000 Data Byes, 1000 Points
                    FanHeader.NReturnBytes = 1013;
                    FanHeader.NZeroFillBytes = 39;
                elseif switches.DataBits == 14, % 2000 Data Byes, 1000 Points
                    FanHeader.NReturnBytes = 2013;
                    FanHeader.NZeroFillBytes = 63;
                end
            case 15,
                FanHeader.NPoints = 1500;
                if switches.DataBits == 4, % 750 Data Byes, 1500 Points
                    FanHeader.NReturnBytes = 763;
                    FanHeader.NZeroFillBytes = 33;
                elseif switches.DataBits == 8, % 1500 Data Byes, 1500 Points
                    FanHeader.NReturnBytes = 1513;
                    FanHeader.NZeroFillBytes = 51;
                elseif switches.DataBits == 14, % 3000 Data Byes, 1500 Points
                    FanHeader.NReturnBytes = 3013;
                    FanHeader.NZeroFillBytes = 87;
                end
            otherwise % >= 20 m
                FanHeader.NPoints = 2000;
                if switches.DataBits == 4, % 1000 Data Byes, 2000 Points
                    FanHeader.NReturnBytes = 1013;
                    FanHeader.NZeroFillBytes = 39;
                elseif switches.DataBits == 8, % 2000 Data Byes, 2000 Points
                    FanHeader.NReturnBytes = 2013;
                    FanHeader.NZeroFillBytes = 63;
                elseif switches.DataBits == 14, % 4000 Data Byes, 2000 Points
                    FanHeader.NReturnBytes = 4013;
                    FanHeader.NZeroFillBytes = 111;
                end
        end
    case {'P'},  % Profile = ON
        FanHeader.NReturnBytes = 13;  % 0 Data Byes, 0 Points
        FanHeader.NZeroFillBytes = 15;
        FanHeader.NPoints = 0;
    otherwise
        sprintf([FanHeader.ReturnDataHeaderType,' not recognized'])
        return
end
if verbose,
    disp(FanHeader)
end

return


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% read the settings sent to the sonar head
% these are saved at the beginning of each data file in the logger
function switches = parseFanswitches(switchdata, verbose)

% first two bytes are switch data header ID = 0xFE44
if (switchdata(1) ~= 254) || (switchdata(2) ~= 68),
    disp('Switch data header ID not found')
end
switches.HeadID = switchdata(3); % always 0x10
if (switchdata(3) ~= 16),
    disp('Switch data Head ID not 0x10, as it should be for an 881 Fan sonar')
end
switches.Range = switchdata(4); % 5 to 255 m in 5 m increments
switches.Reserved(1) = NaN;
switches.RangeOffset = switchdata(5); % 0 to 255 in 1 m increments
% return data is from RangeOffset to RangeOffset + Range
b=dec2base(switchdata(6),2,8); % get the bits
if b(1), switches.RevHold = 'Hold Head';
else switches.RevHold = 'Resume';
end
if b(7), switches.RevHold = 'Reverse Step Direction';
else switches.RevHold = 'Normal Operation';
end
b=dec2base(switchdata(7),2,8); % get the bits
% I'm not too sure how to interpret the manual on this one
% but this should not matter to logger data
if b(1),
    if b(7), switches.MasterSlave = 'Transmit Slave Mode';
    else switches.MasterSlave = 'Transmit Master Mode';
    end
end
if b(2),
    if b(7), switches.MasterSlave = 'Send Data Slave Mode';
    else switches.MasterSlave = 'Send Data Master Mode';
    end
end
switches.Reserved(2) = switchdata(8); % always 0
switches.StartGain = switchdata(9); % 0 to 40dB in 1 dB increments
switches.LOGF = switchdata(10).*10; % in dB
switches.Absorption = switchdata(11)./100; % dB/m
if switches.Absorption == 253,
    disp('Absorption value of 253 detected, this is not allowed')
end
switches.TrainAngle = (switchdata(12).*3)-210; % degrees from 0 to 140
switches.SectorWidth = switchdata(13).*3; % degrees from 0 to 120
% fan beam step increments are 0.15 degrees
switches.StepSize = switchdata(14).*0.15; % degrees per step up to 2.4
if switches.StepSize > 2.4,
    disp('Step size greater than 2.4 detected, this is not possible')
end
switches.PulseLength = switchdata(15).*10; % microseconds
switches.Reserved(3:4) = switchdata(16:17); % always 0
switches.StepDirection = bitget(switchdata(18), 8); % 0 = CCW
switches.MoveRelative = bitand(127, switchdata(18)); % 1 to 127 in 0.15 degree increments
switches.Reserved(5) = NaN;
switches.Reserved(6) = switchdata(19); % always 0
switches.DataPoints = switchdata(20).*10; % 250 for IMX and 500 for IGX and 2000 for IQX
% resolution of data as data width...
% 4 bits = 2 data points per byte
% 8 bits = 1 data point per byte
% 16 bits = 2 bytes per data point
switches.DataBits = switchdata(21); % in bits
switches.UpBaud = switchdata(22); % not yet used
if switchdata(23), switches.Profile = 'ON';
else switches.Profile = 'OFF'; % return data will have IPX header
end
if switchdata(24), switches.Calibrate = 'Normal Operation';
else switches.Calibrate = 'Calibrate'; % head is moved to 0 degrees
end
% the delay to give a control program time to serially upload data
switches.SwitchDelay = switchdata(25).*2; % milliseconds
if switches.SwitchDelay == 253,
    disp('Switch Delay value of 253 detected, this is not allowed')
end
switches.Reserved(7) = switchdata(26); % always 0
% termination byte
if (switchdata(27) ~= 253),
    disp('Switch data termination byte not found')
end
if verbose,
    disp(switches)
end
return