function [imagedata, headangle, prange] = showpencil(fname, datayear);
%
% function [imagedata, headangle, prange] = showpencil(fname, datayear);
%
% display Imagenex Pencil Beam Sonar images taken with Robin Singer's data logger
%
% fname = file name
% datayear = the year the data was taken
% imagedata = the image as [angle, point]
% headangle = the head position in degrees
% prange = profile range in cm
%

% Written by Marinna Martini 8/28/02
% for the U.S. Geological Survey, WOods Hole Field Center
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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."
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%

if exist('fname') ~= 1,
    [fname, fpath] = uigetfile('*.p*','Pencil beam data');
else
    fpath = [];
end

if 0,
if exist('datayear') ~= 1,
    prompt  = {'Enter the year the data was recorded:'};
    windowtitle   = 'Missing data';
    lines= 1;
    def     = {'2003'};
    answer  = inputdlg(prompt,windowtitle,lines,def);
    datayear = str2num(answer{1});
end
end
datayear = 2002;
disp('assuming we are in the year 2002, if not, use the command line call')

fid = fopen([fpath, fname],'r');
if fid == -1,
    disp(['error reading file ',fpath, fname]);
    return
else
    disp(['Now reading file ', fpath, fname]);
end

% standard filename format is S1MMDDHH.PMM = S1082715.p30
gdate = [datayear str2num(fname(3:4)) str2num(fname(5:6)) str2num(fname(7:8)) str2num(fname(11:12)) 0];
jdate = julian(gdate);

nheaderbytes = 13;

% skip the 61 bytes of file header
[fheader, numread] = fread(fid,61,'uchar');

% decipher the first 13 bytes of sonar header so we know what we are reading
[sheader, numread] = fread(fid,nheaderbytes,'uchar');
disp('The header')
sheader'

% imagenex return data header will tell us the total number of return data bytes
% conver the first three bytes into an ASCII string
return_data_header = char(sheader(1:3,:))';
disp(sprintf('Return Data Header = %s',return_data_header));
head_ID = dec2hex(sheader(4));
disp(sprintf('Header ID = %s',head_ID))
% we can skip the serial status for now
% head position
headposhi = bitand(63,sheader(7));
headposlo = bitand(127,sheader(6));
headpos = headposhi*127+headposlo;
stepdir = bitget(sheader(7),6);
% translate from 0 to 2800 meaning -210 to 210 in 0.15 degree steps
headangle = 0.15*(headpos-1400);
disp(sprintf('Head angle = %f', headangle))
% sonar head range
sonar_head_range = sheader(8);
disp(sprintf('Sonar Head Range setting = %d m',sonar_head_range));
% profile range
prangehi = bitand(127,sheader(10));
prangelo = bitand(127,sheader(9));
prange = prangehi*127+prangelo;
disp(sprintf('First Digitized Range Value Above Threshold = %d cm', prange))
% Number of echo data bytes returned
ndatabyteshi = bitand(127,sheader(12));
ndatabyteslo = bitand(127,sheader(11));
ndatabytes = ndatabyteshi*127+ndatabyteslo;
disp(sprintf('Number of data bytes = %d', ndatabytes))

% 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)
ndatabits = 8;

switch(return_data_header),
case {'IMX'},   % switch data command for data points = 25
    if ndatabits == 4, nreturnbytes = 141; end  % 128 Data Byes, 256 Points
    if ndatabits == 8, nreturnbytes = 265; end  % 252 Data Byes, 252 Points
    if ndatabits == 14, nreturnbytes = 513; end % 500 Data Byes, 250 Points
case {'IGX'},   % switch data command for data points = 50
    if ndatabits == 4, 
        nreturnbytes = 265; % 252 Data Byes, 504 Points
    else
        nreturnbytes = 513; % 500 Data Byes, 500 Points
                            % 500 Data Byes, 250 Points
    end
case {'IQX'}, 
    % there is no IQX for the pencil beam sonar
case {'IPX'},  % Profile = ON
    nreturnbytes = 13;  % 0 Data Byes, 0 Points
end
disp(sprintf('Number of return bytes per ping = %d',nreturnbytes));
npoints = nreturnbytes - nheaderbytes;

% now figure out how many pings of data are in this file from the file size and header info.
% 61 bytes of file header
fseek(fid,0,'bof'); 
filebeg = ftell(fid);
fseek(fid,0,'eof'); % go to the end of the file
filebytes = ftell(fid)-filebeg-118; % byte from beginning of file

npings = round(filebytes./nreturnbytes);
disp(sprintf('Number of useable pings = %d (maybe)',npings))

frewind(fid);

% skip the 61 bytes of file header
[fheader, numread] = fread(fid,61,'uchar');

% guess at a dataholdersize
imagedata = zeros(npings, nreturnbytes-nheaderbytes);
headangle = zeros(npings, 1);
prange = zeros(npings,1);
pingidx = 0;

while 1,
    pingstart = ftell(fid);
    if pingstart < 0, break; end
    pingidx = pingidx+1;
    
    % read header
    [sheader, numread] = fread(fid,nheaderbytes,'uchar');
    if numread < nheaderbytes, break; end
    % grab the position
    % head position
    headposhi = bitand(63,sheader(7));
    headposlo = bitand(127,sheader(6));
    headpos = headposhi*127+headposlo;
    stepdir = bitget(sheader(7),6);
    % translate from 0 to 2800 meaning -210 to 210 in 0.15 degree steps
    headangle(pingidx) = 0.15*(headpos-1400);
    % profile range
    prangehi = bitand(127,sheader(10));
    prangelo = bitand(127,sheader(9));
    prange(pingidx) = prangehi*127+prangelo;
        
    % read image data
    [pingdata, numread] = fread(fid,nreturnbytes-nheaderbytes,'uint8');
    if numread < nreturnbytes-nheaderbytes, break; end
    imagedata(pingidx,:) = pingdata';
        
end
fclose(fid);

[npings, npoints] = size(imagedata);
headangle = headangle(1:npings);    % trim it if there's an extra point
prange = prange(1:npings);

% plot linear 
figure
clf;
imagesc(headangle,0:(sonar_head_range/npoints):sonar_head_range,imagedata')
colormap(jet);
colorbar
ylabel('Range, meters')
xlabel('Head Position, degrees')
title(['Fan beam image taken ',datestr(datenum(gdate))])

% plot headposition and prange
figure
[ax, h1, h2] = plotyy(1:npings,headangle,1:npings,prange);
title(['Fan beam image taken ',datestr(datenum(gdate))])
xlabel('Head position and profile range by ping')
ylabel('Deg')