function [PencilData, settings] = showpen07(fname, settings);
% SHOWPEN_07 - Display Imagenex Pencil Beam Sonar images taken with Robin Singer's data logger
% [imagedata, headangle, prange] = showpencil(fname, settings);
%
% 
% 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
%

% etm mods 4/07 
%  works with data from 3/14/07 tests with known bottom shapes
% assumes .3 degrees/step- should get this from the header if possible
% uses Doug Wilson's method of reading the data and doing the initial
% plot.  To converto to x-space, griddata is used.  Then the max value of
% each angular scan is used to define the bottom.  Simplistic, the the line
% generated agrees with the color image displayed (if you un comment the
% pcolor(Xr, -Yr, D); shading flat line...
%
% Revisions by Chris Sherwood 23 June 2006
% Written by Marinna Martini 8/28/02
% for the U.S. Geological Survey, WOods Hole Field Center

nheaderbytes = 13;

fid = fopen(fname,'r');
if fid == -1,
    disp(['error reading file ',fname]);
    return
else
    disp(['Now reading file ', fname]);
end
[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 ( this is the same computation as below for hp1
 headposhi = bitand(63,sheader(7));
 headposlo = bitand(127,sheader(6));
  headpos = (headposhi.*128)+headposlo;
   hp=bitand(headpos,63);
   stepdir = bitget(sheader(7),6);
% we're doing .3 degree steps and not sure what -600 is for but it works
  headangle = 0.3*(headpos-600);
  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);
% this is how Doug W reads in the data
   A1 = fread(fid,'uchar');
   lA1 = length(A1);
   A1t = A1(62:lA1);
   lA1t = length(A1t);
   npings1 = lA1t/513;
   imagedata = reshape(A1t,513, npings1);
   hp1  = bitand(imagedata(6,:),127) + 128.*bitand(imagedata(7,:), 63);  % high bit is zero, and
% bit 6 is direction of scan, so remove it
   hpangledeg1 = (0.3.*(hp1 - 600));	% Angle in radians for Matlab
   hang = hpangledeg1.*pi./180;
   range_config=imagedata(8);
fclose(fid);
% head angle should be -75 -> 75; 502 points is correct

% remove beginning part nearest the xducer
imagedata=imagedata(14:end,:);

  npoints = 500;
  SampPerMeter = npoints/range_config;

   r  = 1:npoints;    %because the last point is already removed from imagedata
   rr = r./SampPerMeter(1);
	Yr = rr'*cos(hang);
	Xr = rr'*sin(hang);

   colorthreshold = 0;
   %colorthreshold = input('Color threshold for black (1:100)? ');  % This is a noise floor  
   D = imagedata-colorthreshold;  	% Trim off header
   Dz = find(D<0);									% Make all values below
   D(Dz) = 0;											% threshold = 0
   %figure
   %axhan = pcolor(Xr,-Yr,D); shading interp	% you can use surf or mesh if desired here
   %axis ([-3 3 -.9 -.7])
   %title('raw version, plotted against sin and cos')

% interpolate into x,y space
% xx and yy below are arbitrary numbers for getting the right general shape- they can be tweaked!
xx=[-3.16:.0125:3.16];
yy=[.2:.0025:1.4]';
imi=griddata(Xr,Yr,D,xx,yy,'linear');
 %figure
 %pcolor(xx,-yy, imi); shading flat
 %title('interpd version')

  % find the elevations using the max return in each bin
  zz=max(imi);
  for ik=1:length(zz)
   if ~isnan(zz(ik))
    nn=find(imi(:,ik)==zz(ik),1,'first');
   else
    nn=1;
   end
  ly(ik)=nn;
  end
  
  % use the contiguous middle of the plot to select the x range
  % figure(2); plot(diff(yy(ly)))
   jmps=find(abs(diff(yy(ly))) >.2);
      if(jmps(end-1)-jmps(2) < range_config.*100)
        strt_idx=1; end_idx=jmps(end);
      else
        strt_idx=jmps(2)+1;
        end_idx=jmps(end-1)-1;
      end
   x_gd=xx(strt_idx:end_idx);
   y_gd=-yy(ly(strt_idx:end_idx));
   z_gd=zz(ly(strt_idx:end_idx));
    
    %now remove everything greater than 3 std_devs of the mean
    mn_el=mean(y_gd);
    std_el=std(y_gd);
    gdvals=[mn_el-(3*std_el) mn_el+(3*std_el)];
    ng_idx=find(y_gd < gdvals(1) | y_gd > gdvals(2));
    y_gd(ng_idx)=NaN; 
     figure
      plot(x_gd,y_gd,'.')
       title([fname(end-11:end) '; range setting= ' num2str(range_config)])
       xlabel('horizontal distance along seafloor(m)')
       ylabel('depth(m)')
  % now save the values you want 
PencilData.imagedata = D;
PencilData.xmat = Xr;
PencilData.ymat = Yr;
PencilData.range_config = range_config;
PencilData.headangle = hang;
PencilData.filename = fname;
PencilData.intrpx = x_gd;
PencilData.intrpy = y_gd;
PencilData.intrpz = z_gd;