function [imagedata, headangle, prange, settings] = showpencil(fname, plottype, settings);
% SHOWPENCIL.M     DISPLAY IMAGENEX PENCIL BEAM DATA COLLECTED BY USGS PERSISTOR DATA LOGGER
% function [imagedata, headangle, prange, settings] = showpencil(fname, plottype, settings);
%
% fname = file name
% imagedata = the image as [angle, point]
% headangle = the head position in degrees
% prange = profile range in cm
% metadata:
% settings.SectorSweep 
% settings.StepSize 
% settings.DataPoints 
% settings.Height 
% settings.datayear
% plottype = 'square' | 'polar' | 'prange'

% Written by Marinna Martini 8/28/02
% for the U.S. Geological Survey, WOods Hole Field Center
% updated 2/18/03 on the PASTA turn around

% a test file from Eurostrat mooring 701
%fname = 'C:\projects\euro\701im\I1204\S1120301.F50'

if exist('fname') ~= 1,
    [fname, fpath] = uigetfile('*.p*','Pencil beam data');
    if fname == 0, return, end
else
    [fpath,fname,ext,ver] = fileparts(fname);
    fname = [fname ext];
    if ~isempty(fpath), fpath = [fpath '\']; end
end

if exist('plottype') ~= 1,
    plottype = 'polar';
end

if exist('settings') ~= 1,
    % get the metadata from the user
    prompt  = {'Enter the sector that was swept:',...
            'Enter the center to sweep around:',...
            'Enter the Step size used:',...
            'Enter the data points setting (points/10):',...
            'Enter the sonar''s height off the bottom, m:',...
            'Enter the year the instrument was deployed:',...
        };
    % settings for Eurostrat fall deployment, mooring 701
    def     = {'132','0','1','500','1.06','2002'};
    dlgtitle   = 'PENCIL BEAM Input metadata for the setup';
    lineNo  = 1;
    dlgresult  = inputdlg(prompt,dlgtitle,lineNo,def);
    settings.SectorSweep = str2num(dlgresult{1});
    settings.AngleSweepAround = str2num(dlgresult{2});
    settings.StepSize = str2num(dlgresult{3});
    settings.DataPoints = str2num(dlgresult{4});  % number of 'ranges' recorded in each ping
    settings.Height = str2num(dlgresult{5});
    settings.datayear = str2num(dlgresult{6});
end
% pencil beam specific things
settings.DegPerStep = 0.3;
%settings.DegPerStep = 0.15;

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 = [settings.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
settings.SonarRange = sheader(8);
%disp(sprintf('Sonar Head Range setting = %d m',settings.SonarRange));
% 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);

slantrange = settings.SonarRange/npoints:(settings.SonarRange/npoints):settings.SonarRange;
%sectorswept = 1:settings.DegPerStep:settings.SectorSweep+settings.DegPerStep+1;
startangle = settings.AngleSweepAround-settings.SectorSweep/2;
endangle = settings.AngleSweepAround+settings.SectorSweep/2;
sectorswept = startangle:settings.DegPerStep:endangle+settings.DegPerStep;

% trim the image to useable area
[nangles, npoints] = size(imagedata);
% trim the redundant first and last ping
imagedata = imagedata(2:nangles-1,:);
headangle = headangle(2:nangles-1);    
prange = prange(2:nangles-1);
sectorswept = sectorswept(2:nangles-1);
% trim range to real area of data
nearfieldcutoff = 1;
farfieldcutoff = npoints;
imagedata = imagedata(:,nearfieldcutoff:farfieldcutoff);
vertrange = slantrange(nearfieldcutoff:farfieldcutoff);
% now get the new size for plotting
[nangles, npoints] = size(imagedata);
 
if exist('plottype') == 1,
    switch plottype
        case 'square'
            % plot linear 
            set(gcf,'name','Square')
            set(gcf,'numbertitle','off')
            %Z = log10(imagedata+eps)';
            %imagesc(headangle,vertrange,Z);
            %clims = [0 190];
            imagesc(headangle,vertrange,imagedata'); %,clims)
            colormap(bone);
            colorbar
            ylabel('Horizontal Range, meters')
            xlabel('Head Position, degrees')
            text(0.01,0.05,[fname,' Pencil beam image taken ',datestr(datenum(gdate))],...
                'units','normalized','color','w')
        case 'polar'
            set(gcf,'name','Polar')
            set(gcf,'numbertitle','off')
            set(gcf,'position',[1 29 1000 330])
            %[th, r] = meshgrid(sectorswept.*(pi/180), vertrange);
            [th, r] = meshgrid((sectorswept-90).*(pi/180), vertrange);
            [X, Y] = pol2cart(th, r);
            Z = imagedata';
            %Z = log10(imagedata+eps)';
            surf(X, Y, Z)
            shading flat
            colormap jet
            axis equal tight
            set(gca,'ylim',[-settings.Height-0.25 0]);
            set(gca,'xlim',[-3 3]);
            view([0 90]); % az, el
            set(gca,'position',[0.05 0.3 0.9 0.5])
            title([fname,' Pencil beam image taken ',datestr(datenum(gdate))])
            ylabel('Range, m')
            xlabel('Width, m')
        case 'prange'
            % plot headposition and prange
            [ax, h1, h2] = plotyy(1:nangles,headangle,1:nangles,prange);
            title([fname,' Fan beam image taken ',datestr(datenum(gdate))])
            xlabel('Head position and profile range by ping')
            ylabel('Deg')
            %set(ax(2),'ylabel',text(0,0,'string','Prange'))
        otherwise % no plot
    end
end