function [imagedata, Xplot, Yplot, Zi, settings] = ...
   showfan(fname, plottype, settings );
% SHOWFAN - Plot Imagenex imaging sonar data
% [imagedata, headangle, prange, settings] = showfan(fname, plottype, settings)
%
% s.fname = file name
% s.imagedata = the image as [angle, point]
% s.headangle = the head position in degrees (from the Sonar)
% s.fakeangle = what we think it is, in true
% s.prange = profile range in cm
% s.dn = Matlab datenum time
% metadata:
% s.settings:
% settings.SectorSweep 
% settings.StepSize 
% settings.DataPoints
% settings.Height 
% settings.Datayear
% plottype = 'square' | 'polar' | 'prange'

%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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 Marinna Martini 8/28/02
% for the U.S. Geological Survey, WOods Hole Field Center
% updated 2/18/03 on the PASTA turn around

% - added slant range, angle, and fake angle to return values;
% put them in a structure.
% - add polar plot
% Notes to Marinna and Peter Howd:
%  There are a couple of weird things:
%  1) imagedata(:,449) = 252 all the time. I just blank it out...is that
%  one point too far out the range? 
%  2) the real head angle goes retrograde sometimes. This is either a weird reading
%  in the file, or maybe an integer math thing with the angle calcs. I
%  make fake angles. 
% Last revised by Chris Sherwood, USGS, 3/5/04
% Last revised by Charlene Sullivan, USGS, 9/15/05

% Test values:

if(nargin == 0),
  fname = 'E:\CSULLIVAN\SOUTH_CAROLINA\DATAFILES\752\7523im\I0213\S1013120.F20';
  %fname = 'E:\CSULLIVAN\SOUTH_CAROLINA\DATAFILES\732\7323im\I1113\S1111214.F50';
  %fname = 'E:\EuroStrat\701im\images_for_movies\I1117\S1111608.F20'; % don't delete these
  %fname = 'G:\Feb03EuroStrat\701-Particle\701im\I1124\S1112313.F50' 
  %fname = 'C:\crs\data\EuroStrat_0203\701-Particle\701im\I1124\S1112313.F50'
  plottype = 'polar';
  %plottype = 'square';
end

if exist('fname') ~= 1,
    [fname, fpath] = uigetfile('c:\projects\euro\701im\*.f*','Fan beam data');
    if fname == 0, return, end
else
    [fpath,fname,ext,ver] = fileparts(fname);
    fname = [fname ext];
    fpath = [fpath '\'];
end

if exist('plottype') ~= 1,
    plottype = 'square';
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     = {'348','0','1','500','0.55','2002'};
    % settings for Myrtle Beach deployment are same as Eurostrat except for
    % the year
    def     = {'348','0','1','500','0.55','2003'};
    dlgtitle   = 'FAN 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
% fan beam specific things
%settings.DegPerStep = 0.3;
settings.DegPerStep = 0.15;

pingsperrec = ((settings.SectorSweep./(settings.StepSize.*settings.DegPerStep))+1)*2;

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));
if ~strcmp(head_ID,'10'),
    disp('This is not a Model 881-000-100 360 degree scanning head sonar!')
    disp('This program cannot read this data.')
    return
else
%    disp('Verified that data is from a Model 881-000-100 360 degree scanning head sonar')    
end
% 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'},   % switch data command for data points = 200
    switch(range),
    case 5,
        if ndatabits == 4, nreturnbytes = 263; end  % 250 Data Byes, 500 Points
        if ndatabits == 8, nreturnbytes = 513; end  % 500 Data Byes, 500 Points
        if ndatabits == 14, nreturnbytes = 1013; end  % 1000 Data Byes, 500 Points
    case 10,
        if ndatabits == 4, nreturnbytes = 513; end  % 500 Data Byes, 1000 Points
        if ndatabits == 8, nreturnbytes = 1013; end  % 1000 Data Byes, 1000 Points
        if ndatabits == 14, nreturnbytes = 2013; end  % 2000 Data Byes, 1000 Points
    case 15,
        if ndatabits == 4, nreturnbytes = 763; end  % 750 Data Byes, 1500 Points
        if ndatabits == 8, nreturnbytes = 1513; end  % 1500 Data Byes, 1500 Points
        if ndatabits == 14, nreturnbytes = 3013; end  % 3000 Data Byes, 5500 Points
    otherwise,
        if ndatabits == 4, nreturnbytes = 1013; end  % 1000 Data Byes, 2000 Points
        if ndatabits == 8, nreturnbytes = 2013; end  % 2000 Data Byes, 2000 Points
        if ndatabits == 14, nreturnbytes = 4013; end  % 4000 Data Byes, 2000 Points
    end
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

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

frewind(fid);

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

% guess at a dataholdersize
imagedata = zeros(nangles, nreturnbytes-nheaderbytes);
headangle = zeros(nangles, 1);
prange = zeros(nangles,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);
disp('Done reading file.')
settings
slantrange = settings.SonarRange/npoints:(settings.SonarRange/npoints): ...
    settings.SonarRange;
sectorswept = 1:settings.DegPerStep:settings.SectorSweep+settings.DegPerStep+1;

% 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 = 52;
farfieldcutoff = npoints;
imagedata = imagedata(:,nearfieldcutoff:farfieldcutoff);
slantrange = slantrange(nearfieldcutoff:farfieldcutoff);
% now get the new size for plotting
[nangles, npoints] = size(imagedata);

% other calculations
%horizrange = slantrange.*acos(asin((settings.Height./slantrange)));
%horizrange = settings.Height.*cos(asin(slantrange./settings.Height));
horizrange = real((slantrange.^2 - settings.Height^2).^(1/2));

if exist('plottype') == 1,
  switch plottype
   case 'square'
    % plot linear 
    set(gcf,'name','Square')
    set(gcf,'numbertitle','off')
    clims = [0 190];
    imagesc(headangle,horizrange,imagedata',clims)
    colormap(bone);
    colorbar
    ylabel('Horizontal Range, meters')
    xlabel('Head Position, degrees')
    text(0.01,0.05,[fname,' Fan beam image taken ',...
		    datestr(datenum(gdate))],...
	 'units','normalized','color','w')
   case 'polar'
    settings.fillval = 4; % sets value of no-data region
		 % median value is 8, mean is 12, range is 0 - 32
    % rotating by 180-17 will point EuroSTRATAFORM ADCP (beam 3) to north
    % rot = 180-17;
    % during first deployment, ADCP heading was 41 deg M = 43 deg T:
    % rot = 180-17+43
    % during second deployment, ADCP heading was 274 deg M = 276 deg T:
    %rot = 180-17+(276-360)
    %rot = settings.rot;
    % we don't yet know direction of north for Myrtle Beach deployment, so
    % for now set rot = 0;
    %rot = 0;
    % fall Myrtle Beach deployment ADCP heading was 149 deg M = 141 deg T
    rot = settings.rot;
    % spring Myrtle Beach deployment ADCP heading was 253 deg M = 245 deg T
    %rot = settings.rot;
    %
    
    settings.max_range = 5;
    settings.dxy = 0.01;    % Key setting...determines image resolution at cost of
                   % speed (reasonable range 0.02 to 0.005)

    Xplot = (- settings.max_range:settings.dxy:settings.max_range);
    Yplot = Xplot;
    [Xi, Yi] = meshgrid(Xplot,Yplot);
    [thi, ri]=cart2pol(Xi,Yi);
    thi = thi*180/pi;
    thi = thi+rot;
    thi(thi>180)=thi(thi>180)-360;
    thi(thi< -180)=thi(thi< -180)+360;
    fakeheadangle = linspace(headangle(1), headangle(nangles), ...
			     nangles)';
    fakeheadangle = -fakeheadangle; % because sonar scans from -160 to +170
    % pad ranges
    slantrange = [0,slantrange(1,1)-.1,...
		  slantrange(1,:), slantrange(1,npoints)+.1, 9];
    % pad angles to get complete circle
    lopad = [180; fakeheadangle(1)+.01];
    hipad = [fakeheadangle(nangles)-.01; -179.999];
    nlo = length(lopad);
    nhi = length(hipad);
    fakeheadangle = [lopad(:); fakeheadangle; hipad(:)];
    % pad data
    Zpad = (settings.fillval)*ones(nangles+nlo+nhi,npoints+4);       
    imagedata(imagedata<(1))=(settings.fillval);
    % kludge...last range values always seems to be 252
    imagedata(:,449)=(settings.fillval);
    Zpad(nlo+1:nangles+nlo,3:npoints+2)=imagedata;
    Zpad = log10(Zpad);
    % interpolate
    [ripad,thpad] = meshgrid(slantrange,fakeheadangle);
    Zi = interp2( ripad,thpad,Zpad, ri, thi);
%    Zi(Zi<1)=(8); % produces weird contrasty plot
    imagesc(Xplot,Yplot,Zi)
    set(gca,'ydir','normal') % because imagesc is sidedownup
    colormap gray;
    set(gca,'Fontsize',12);
    set(gca,'Xtick',[-5:1:5])
    set(gca,'Ytick',[-5:1:5])
    %axis equal tight
    if(1), % plot range circles
      hold on
      for icirc = 1:5,
	h = circle(icirc,0,0);
	set(h,'color',[.7 .7 .6])
      end
    end
    text(0.01,0.03,fname,...
	 'units','normalized','color','y','horizontalalignment','left','fontsize',10)
    text(0.8,0.95,'\uparrow North',...
	 'units','normalized','color','y','fontsize',10)
    ts = datestr(datenum(gdate));
    text(0.99,0.03,ts,...
	 'units','normalized','color','y','horizontalalignment','right','fontsize',10)
    
   case 'glacial' % this is really slow
    set(gcf,'name','Polar')
    set(gcf,'numbertitle','off')
    [th, r] = meshgrid(sectorswept.*(pi/180), horizrange);
    [X, Y] = pol2cart(th, r);
    %Z = imagedata';
    Z = log10(imagedata+eps)';
    surf(X, Y, Z)
    view([0 90]); % az, el
    shading flat
    colormap jet
    text(0.01,0.05,[fname,' Fan beam image taken ',datestr(datenum(gdate))],...
                'units','normalized','color','k','fontsize',9)
    %axis equal tight
    %set(gca,'xlim',[0 settings.Height+0.5]);
   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
    %return values
    Xplot = slantrange;
    Yplot = headangle;
    Zi = imagedata;
  end
end

keep3 imagedata Xplot Yplot Zi settings