function plotpencilMM(header, image, plottype)
% Written by Marinna Martini 8/28/02
% for the U.S. Geological Survey, WOods Hole Field Center

settings.pencil_vert_ang = 0.0;
header.Height = 1.5;

imagedata = image.imagedata;
[npoints, nangles] = size(imagedata');
headangle = header.HeadAngle.*pi./180;  % Angle in radians for Matlab
nearfieldcutoff = 1; %52;
farfieldcutoff = header.NPoints;
slantrange = header.Range(1)/header.NPoints:(header.Range(1)/header.NPoints): ...
    header.Range(1);
slantrange = slantrange(nearfieldcutoff:farfieldcutoff);
horizrange = real((slantrange.^2 - header.Height^2).^(1/2));
[pathstr,name,ext] = fileparts(header.FileName);
fname = [name, ext];
gdate = gregorian(header.PencilTime(1));
%sectorswept = 1:header.StepSize:header.SectorWidth+header.StepSize+1;

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

switch plottype
    case 'square'
        % plot linear
        set(gcf,'name','Pencil 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,' Pencil beam image taken ',...
            datestr(datenum(gdate))],...
            'units','normalized','color','w')

    case 'polar'
        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
        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')

    case 'edge'
        colorthreshold = 0;
        %colorthreshold = input('Color threshold for black (1:100)? ');  % This is a noise floor
        D = imagedata-colorthreshold;  	% Trim off header
        % 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]';
        % interpolate into x,y space
        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) < header.Range(1).*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(header.Range(1))])
        xlabel('horizontal distance along seafloor(m)')
        ylabel('depth(m)')

    case 'headangle'
        set(gcf,'name','Pencil Head')
        set(gcf,'numbertitle','off')
        plot(headangle)
        ylabel('Head Angle, degrees')
        text(0.01,0.05,[fname,' Pencil beam image taken ',...
            datestr(datenum(gdate))],'units','normalized')
    otherwise % no plot
end