function [pict] = plotpen81a(PenData, HeaderData, plottype);
% PLOTPEN81a - Plot Imagenex imaging sonar data logged by imagenex viewer
%     [pict, x_gd, y_gd, z_gd] = plotpen07(PenData, PenHeader, plottype)
% PenHeader and PenDat are returned by readpencil.m
% plottype = 'square' | 'polar' | 'prange'
% still need to pass 'settings' to get these params: 
% instrument installed height, adcdp compass value, declination
% pict is a structure containing the image from which the elevation and max return are obtained
%  imagesc(pict.x, pict.y, pict.imi) should work

% based on showPen07 by EM
% reading the data is now replaced by readfan()

 SampPerMeter=HeaderData.NPoints/HeaderData.Range(250);  %range settings 
                            % should all be same, but using mid-cval anyhow
 r  = 1:HeaderData.NPoints;   
 rr = r./SampPerMeter(1);
                                       
  % try compensating for transducer angle off the vertical
    settings.Pencil_vert_off=0;
  hd_angle=HeaderData.HeadAngle+settings.Pencil_vert_off;
  %
  % with .81a files, probably have multiple sweeps, so have to separate:
  %  organize into sweeps by when headangle changes
    mx=max(hd_angle);
    mn=min(hd_angle);
    mxind=find(hd_angle==mx);
    mnind=find(hd_angle==mn);    
    % sometimes you get several max and mins before turns around:
      lx=find(diff(mnind)==1);
         mnind(lx)=[];
         clear lx
      lx=find(diff(mxind)==1);
         mxind(lx)=[];
    sweep_strt=sort([mxind mnind]);
     len_swp=sweep_strt(2)-sweep_strt(1);
     
     % scale factor to get the Profile Range correct
     %fctr2=(2/1000)*100;
     rx2=HeaderData.Range(1)*2;
     
     mmpm=0.002;  % this should be .002 2mm/m = 2/1000
    
for ik=1:floor(length(sweep_strt)-1)
    clf
   hang=deg2rad(hd_angle(sweep_strt(ik):sweep_strt(ik+1)-1));  % has to be radians for trig functions                                   
	 Yr = rr'*cos(hang);
	 Xr = rr'*sin(hang);

       raw = PenData.imagedata(1:251,sweep_strt(ik):sweep_strt(ik+1)-1);  
       Dz = find(raw<0);									% Make all values below
       raw(Dz) = 0;											% threshold = 0
  
% Now interpolate into denser x,y space
% xx and yy below are arbitrary numbers for getting the right general shape- they can be tweaked!
xx=[-HeaderData.Range(250):.01:HeaderData.Range(250)];
yy=[.5:.002:2]';
imi=griddata(Xr,Yr,raw,xx,yy,'linear');
 pcolor(xx,-yy, imi); shading flat
 hold on
 %title('interpd version')
 pr=HeaderData.ProfileRange(sweep_strt(ik):sweep_strt(ik+1)-1);
 %plot((pr/SampPerMeter*fctr2).*(sin(hang)),-((pr/SampPerMeter*fctr2).*(cos(hang))),'r.')
 plot((pr/SampPerMeter)/rx2.*(sin(hang)),-((pr/SampPerMeter)/rx2.*(cos(hang))),'y*')
 plot((pr*mmpm).*(sin(hang)),-((pr*mmpm).*(cos(hang))),'r*') % may only work if range == 2 | 4
  xlabel('distance (m)')
  ylabel('depth (m)')
   %gtext('yellow * have (pr/125.5)/4')  % for range == 2
   %gtext('yellow * have (pr/83.667)/6')  % for range == 4
   %gtext('yellow * have (pr/62.75)/8')  % for range == 4
   %gtext('or (pr/SamplesPerMeter)/(2*Range)')
   %gtext('red * are pr*.002 (or pr*(2mm/1000mm/m))')
  t_str=[HeaderData.FileName ': scan ' num2str(ik)];
  title(t_str)
 pict(ik).imi=imi;
 pict(1).x=xx;
 pict(1).y=-yy;
 pause(3)
end