function [pict, btm_depth, val_max, rawvals] = plotpen07(PenData, HeaderData, plottype, settings);
% PLOTPEN07 - Plot Imagenex imaging sonar data
%  called by procsonarMM: [readfan & plotfan07] replace showfan07
% [pict, btm_depth, val_max, rawvals] = plotpen07(PenData, HeaderData, plottype, settings)
%
% readfan outputs structures containing the image, file and scan headers
% PenData.imagedata: [ping x bin double], image data with header removed
% HeaderData(:).ReturnDataHeaderType = IGX/IMX etc.
% HeaderData(:).HeadID = ID number
% HeaderData(:).SerialStatus 
% HeaderData(:).HeadPosition = head position in counts
% HeaderData(:).HeadAngle = head rotation angle in degrees
% HeaderData(:).StepDirection = direction of rotation, -1 = CCW, 1 = CW
% HeaderData(:).Range = range setting
% HeaderData(:).ProfileRange
% HeaderData(:).NDataBytes = number of byte in the ping
% HeaderData(:).NReturnBytes = number of bytes in the ping + header
% HeaderData.FileName = fname;
% HeaderData.NPoints = PenSwitches.DataPoints;
% HeaderData.StepSize = PenSwitches.StepSize;
% HeaderData.SectorWidth = PenSwitches.SectorWidth;
% plottype = 'square' | 'polar' | 'prange'
% still need to pass 'settings' to get these params: 
%  settings.Pencil_tilt = angle Pencil 0 is from vertical
%     May also want to incorporate these, but not used now.
%     instrument installed height, adcdp compass value, declination
% outputs:
%  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
%   pict.x and pict.y are the values used for interpolating the image
%  btm_depth is the elevation (position of maximum return for each scan)
%  val_max is the value of the maximum return for each scan

% 23 jul 07 updated by MM to remove some hardwired stuff and see if it will
%           work with readrangeall data
% based on showPen07 by EM
% reading the data is now replaced by readfan()
mm=find(HeaderData.HeadAngle==max(HeaderData.HeadAngle));
midpoint = mm(1)/2; %range settings should all be same
SampPerMeter=HeaderData.NPoints/HeaderData.Range(midpoint);  

  % see if we have one or more sweeps in the raw file
  imgsPerSweep = floor(length(PenData.imagedata)/mm(1));
  
% Pencil_tilt must be set prior to execution
  hd_angle=HeaderData.HeadAngle+settings.Pencil_tilt;
   hang=(hd_angle/180*pi);  % has to be radians for trig functions                                   
    r  = 1:HeaderData.NPoints;    %because the last point is already removed from imagedata
    rr = r./SampPerMeter(1);
	 Yr = rr'*cos(hang);
	 Xr = rr'*sin(hang);

       raw = PenData.imagedata;  
       Dz = find(raw<0);									% Make all values below
       raw(Dz) = 0;											% threshold = 0
  
   % find the elevations using the max return in each bin
   % this method works great on MVCO when signal is huge  and oversaturated 
   %  doesn't like tst0314
%   [fx,fy]=gradient(raw);
%   [p,q]=size(fy);
%    for(ik=1:q);
%     nn=find(fy(30:end,ik)==0);
%      if length(nn) > 1
%          if length(nn) > 100   %for low signal data
%            loc=find(abs(fy(30:end-10,ik)) >5);
%            ly(ik)=round(mean(loc))+30;
%          else
%           % now test for contiguous-ness
%             loc=find(diff(nn)==1);  
%             ly(ik)=ceil(mean(nn(loc(1:end-1))))+30;
%          end
%      else
%        ly(ik)=nn+30;
%      end
%    end
  %plot(ly(200),Yr(ly(200),200),'rx')
  % hold on
  % plot(Yr(:,200))
 
  %again tst0314 had too large a nan gaps for this to work
%  nanloc=find(isnan(ly));
%   if length(nanloc) < 50
%    if ~isempty(nanloc); 
%     for jj=1:length(nanloc)  %average over the adjacent points
%      vv(jj)=ceil(nanmean(ly(nanloc(jj)-3:nanloc(jj)+3)));
%     end
%      ly(nanloc)=vv;      %replace nan's with averages
%    end
%   end
%  % now find the distances using the raw data not the interp'd:
%  okdat=find(~isnan(ly));
%   for kk=okdat
%      r_elev(kk)=Yr(ly(kk),kk);
%      r_dist(kk)=Xr(ly(kk),kk);
%   end  
%  rawvals.elev=r_elev;
%  rawvals.dist=r_dist;
  rawvals.elev=0;
  rawvals.dist=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(midpoint):.0125:HeaderData.Range(midpoint)];
yy=[.2:.0025:1.4]';

n=2;
for ik = 1:imgsPerSweep
    strt=n; stp=mm(1)*ik;
     if imgsPerSweep > 1
      imdata=raw(:,strt:stp);
    else
      imdata=raw;
     end
  imi=griddata(Xr(:,strt:stp),Yr(:,strt:stp),imdata(:,strt:stp),xx,yy,'linear');
  pcolor(xx,-yy, imi); shading flat
  %title('interpd version')
  pict.imi(ik,:,:)=imi;
  pict.x=xx;
  pict.y=-yy;

 clear ly nn 
 % find the locatioin of the max point
  zz=max(imi);
  for kk=1:length(zz)
   if ~isnan(zz(kk))
    nn=find(imi(:,kk)==zz(kk),1,'first');
   else
    nn=1;
   end
  ly(kk)=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) < HeaderData.Range(midpoint).*100)
        strt_idx=1; end_idx=jmps(end);
      else
        strt_idx=jmps(2)+1;
        end_idx=jmps(end-1)-1;
      end
      [strt_idx end_idx end_idx-strt_idx]
     % this fails if the indices don't match between two scans
     %  x_gd(ik,:)=xx(strt_idx:end_idx);
     %  y_gd(ik,:)=-yy(ly(strt_idx:end_idx));
     %  z_gd(ik,:)=zz(strt_idx:end_idx);
  % pass everything through for the moment
   y_gd(ik,:)=-yy(ly);  % bottom depth
   z_gd(ik,:)=zz;       % value of max return
    
    %now remove everything greater than 4 std_devs of the mean in mid-sweep
    mn_el=nanmean(y_gd(50:stp-50));
    std_el=std(y_gd(50:stp-50));
    gdvals=[mn_el-(4*std_el) mn_el+(4*std_el)];
    ng_idx=find(y_gd(ik,:) < gdvals(1) | y_gd(ik,:) > gdvals(2));
    z_gd(ik,ng_idx)=1e35; 
    y_gd(ik,ng_idx)=1e35;
  n=stp+1;
end
  btm_depth=y_gd;       
  val_max=z_gd;

%      clf
%       plot(x_gd,y_gd,'.')
%        title([HeaderData.FileName(end-11:end) '; range setting= ' num2str(HeaderData.Range(midpoint))])
%        xlabel('horizontal distance along seafloor(m)')
%        ylabel('depth(m)')
%