function plotrange(header, image, switches, plottype)
%
%  makes various styles of plots from a single azimuth drive record
%  inputs are readrangeall outputs
%
% usage plotrange(PencilHeader, PencilData, switches, Plottype)
%  plottype can be : 2draw | 2d | 3draw | 3d | polar
%  used with the new readrangeall that has .ProfileRange under PenHeader
%  MM and EM July 07
%
% this program ONLY used the .ProfileRange returned by Imagenex for the
% bottom- no analysis of the image is done here.
%
% em 7/24/07 - 2draw and polar don't work as I expect


header.Height = 1;   % should be read from header, but isn't in there now.

if exist('plottype','var') ~= 1,
    plottype = 'polar'; % default
end

clf
switch plottype
    case '2draw'  %doesn't work as expected
        % trim
        ntrim = 0;
        imagedata = image.imagedata(:,1:end-ntrim);
        A = header.Height; % height of pencil sweep vertex above nominal sea bed
        beta = header.HeadAngle; % the Pencil Head Angle
        Ro = 0.08; %meters, pencil sweep apex offset from azimuth centerline of rotation
        beta = beta.*(pi/180); % convert to radians
        m = A.*tan(beta); % horizontal distance from sweep apex
        [nAz, nPoints] = size(imagedata);
        for iAz = 1:nAz,
            idx = find(imagedata(iAz,:) == 0);
            imagedata(iAz,idx) = ones(length(idx),1).*NaN;
            plot(m(iAz,:), imagedata(iAz,:))
            %set(gca,'ylim',[120 200]);
            %set(gca,'xlim',[-switches.Range switches.Range]);
            ylabel('Profile Range, Slant, cm')
            xlabel('Distance from Pencil Apex along sweep plane')
            title(sprintf('Azimuth Angle = %d',image.AAngle(iAz)))
            pause(2)
            cla
        end
    case '2d'           % this works well
        % trim
        ntrim = 0;
        ztrim = 5; % in cm
        if switches.Range < 5,
            factor = 0.002;
        else
            factor = 1;
        end
        PR = header.ProfileRange(:,1:end-ntrim).*factor;
        A = header.Height; % height of pencil sweep vertex above nominal sea bed
        beta = header.HeadAngle; % the Pencil Head Angle (in degrees)
        Ro=0.0;         % nothing done with this at the moment
        %Ro = 0.08; %meters, pencil sweep apex offset from azimuth centerline of rotation
        beta = beta.*(pi/180); % convert to radians
        m = A.*tan(beta); % horizontal distance from sweep apex
        [nAz, nPoints] = size(PR);
        for iAz = 1:nAz,
%        for iAz = 18   % 18 had a lot of junk in the water
            idx = find(PR(iAz,:) == 0); % remove zeros                
            nzero = length(idx);
            PR(iAz,idx) = ones(length(idx),1).*NaN;
            z(iAz,:) = PR(iAz,:).*cos(beta(iAz,:));
            zmed = gmedian(z(iAz,:));
            zcut = zmed + ztrim;
            idx = find(z(iAz,:) > zcut); % remove outliers
            ncut = length(idx);
            z(iAz,idx) = ones(length(idx),1).*NaN;
            plot(m(iAz,:), z(iAz,:),'*',m(iAz,:),zmed,'-.')

            if iAz == 1,
                ylims = get(gca,'ylim');
                xlims = get(gca,'xlim');
            else
                %set(gca,'ylim',[header.Height switches.Range]);
                %set(gca,'xlim',[-switches.Range switches.Range]);
                set(gca,'ylim',ylims);
                set(gca,'xlim',xlims);
            end
            set(gca,'ydir','reverse');
            axis([-2.5 2.5 .8 1.4])
            ylabel('Profile Range, cm')
            xlabel('Distance from Pencil Apex along sweep plane')
            text(0.1,0.1,sprintf('zeros = %4.1f percent',(nzero/nPoints)*100),'units','normalized')
            text(0.1,0.2,sprintf('beyond median + %2d cm = %4.1f percent',ztrim, (ncut/nPoints)*100),'units','normalized')
            title(sprintf('scan %d; Az Angle = %d, median Range = %5.1f',iAz, image.AAngle(iAz),zmed))
            pause(.2)
            cla
        end

    case '3draw'
        % trim
        ntrim = 0;
        ztrim = 5; % in cm
        % A is the Apex point of the pencil beam sweep
        % M is the point of measurement where pencil beam intersect sea bed
        Ro = 0.08; %meters, pencil sweep apex offset from azimuth centerline of rotation
        if switches.Range < 5,
            factor = 0.002;
        else
            factor = 1;
        end
        PR = header.ProfileRange(:,1:end-ntrim).*factor;
        x = ones(size(PR)).*NaN;
        y = x;
        z = x;
        A = header.Height; % nominal height of pencil sweep vertex above sea bed
        beta = header.HeadAngle; % the Pencil Head Angle, [nAz, nPen] positions
        beta = beta.*(pi/180); % convert to radians
        m = A.*tan(beta); % horizontal distance from sweep apex to measurement M
        alpha = header.AAngle; % the azimuth rotation angle, [nAz] positions, on x-y plane
        alpha = alpha.*(pi/180); % convert to radians
        gamma = atan(m./Ro); % angle between points A and M
        [nAz, nPoints] = size(PR);
        for iAz = 1:nAz,
            % x,y = the cartesian location of point A on the plane of the sea bed
            x(iAz,:) = sqrt(m(iAz,:).^2+Ro.^2).*sin(gamma(iAz,:)+alpha(iAz));
       %% !! y is wacked       
            y(iAz,:) = sqrt(m(iAz,:).^2+Ro.^2).*cos(gamma(iAz,:)+alpha(iAz));
            %plot(x(iAz,:),y(iAz,:),'.')
            %hold on
            % z = distance of sea bed surface from point A
            % remove noise
            idx = find(PR(iAz,:) == 0); % remove zeros
            nzero = length(idx);
            PR(iAz,idx) = ones(length(idx),1).*NaN;
            z(iAz,:) = PR(iAz,:).*cos(beta(iAz,:)); % measured distance from apex A to bed
            zmed = gmedian(z(iAz,:));
            zcut = zmed + ztrim;
            idx = find(z(iAz,:) > zcut); % remove outliers
            ncut = length(idx);
            z(iAz,idx) = ones(length(idx),1).*NaN;           
        end
        % z is OK
        surf(x,y,-z); shading flat;
        view (0,90)
        colorbar
        
        %surfc(x,y,z); 

    case '3d'
        % trim
        ntrim = 0;
        ztrim = 1; % in m
        % A is the Apex point of the pencil beam sweep
        % M is the point of measurement where pencil beam intersect sea bed
        %Ro = 0.08; %meters, pencil sweep apex offset from azimuth centerline of rotation
        Ro = 0.1;
        if switches.Range < 5,
            factor = 0.002;
        else
            factor = 1;
        end
        PR = header.ProfileRange(:,1:end-ntrim).*factor;
        x = ones(size(PR)).*NaN;
        y = x;
        z = x;
        A = header.Height; % nominal height of pencil sweep vertex above sea bed
        beta = header.HeadAngle; % the Pencil Head Angle, [nAz, nPen] positions
        beta = beta.*(pi/180); % convert to radians
        m = A.*tan(beta); % horizontal distance from sweep apex to measurement M
        alpha = header.AAngle; % the azimuth rotation angle, [nAz] positions, on x-y plane
        alpha = alpha.*(pi/180); % convert to radians
        gamma = atan(m./Ro); % angle between points A and M
        [nAz, nPoints] = size(PR);
        for iAz = 1:nAz,
            % x,y = the cartesian location of point A on the plane of the sea bed
            x(iAz,:) = sqrt(m(iAz,:).^2+Ro.^2).*sin(gamma(iAz,:)+alpha(iAz));
            y(iAz,:) = sqrt(m(iAz,:).^2+Ro.^2).*cos(gamma(iAz,:)+alpha(iAz));
            % z = distance of sea bed surface from point A
            % remove noise
            idx = find(PR(iAz,:) == 0); % remove zeros
            nzero = length(idx);
            PR(iAz,idx) = ones(length(idx),1).*NaN;
            z(iAz,:) = PR(iAz,:).*cos(beta(iAz,:)); % measured distance from apex A to bed
            zn=z;
            zmed = gmedian(z(iAz,:));
            zdp = zmed + ztrim;
            zsh = zmed - ztrim;
            idx = find(z(iAz,:) > zdp | z(iAz,:) < zsh); % remove outliers
            ncut = length(idx);
            % for now, this cuts out too much
            % z(iAz,idx) = ones(length(idx),1).*NaN;  
            % if iAz == 18; keyboard: end   % sweep 18 had early returns 
        end
        
      % extent of plot is equal to the width of the pencil sweep
        %dxy = 0.005;
        dxy = 0.05;
        sr = max(max(x));
        %Xplot = min(min(x)):dxy:max(max(x));
        Xplot = -sr:dxy:sr;
        Yplot = Xplot;
        [Xi, Yi] = meshgrid(Xplot,Yplot);
        % then try Z, which has more fancy trimming and filling
        Zi = griddata(x,y,z, Xi,Yi);        
        %imagesc(Xplot,Yplot,Zi);
        pcolor(Xplot,Yplot,-Zi); shading flat;
        axis square
        colorbar
     xlabel('meters')
     ylabel('meters')     
     title('image interpolated onto a regular grid: nan for median +/- 1m')

    % does't work - under development    
    case 'polar'
        % trim
        ntrim = 0;
        ztrim = 20; % in cm
        % A is the Apex point of the pencil beam sweep
        % M is the point of measurement where pencil beam intersect sea bed
        Ro = 0.08; %meters, pencil sweep apex offset from azimuth centerline of rotation
        PR = header.ProfileRange(:,1:end-ntrim);
        x = ones(size(PR)).*NaN;
        y = x;
        z = x;
        A = header.Height; % nominal height of pencil sweep vertex above sea bed
        beta = header.HeadAngle; % the Pencil Head Angle, [nAz, nPen] positions
        beta = beta.*(pi/180); % convert to radians
        m = A.*tan(beta); % horizontal distance from sweep apex to measurement M
        alpha = image.AAngle; % the azimuth rotation angle, [nAz] positions, on x-y plane
        alpha = alpha.*(pi/180); % convert to radians
        gamma = atan(m./Ro); % angle between points A and M
        [nAz, nPoints] = size(PR);
        for iAz = 1:nAz,
            % x,y = the cartesian location of point A on the plane of the sea bed
            x(iAz,:) = sqrt(m(iAz,:).^2+Ro.^2).*sin(gamma(iAz,:)+alpha(iAz));
            y(iAz,:) = sqrt(m(iAz,:).^2+Ro.^2).*cos(gamma(iAz,:)+alpha(iAz));
            % z = distance of sea bed surface from point A
            % remove noise
            idx = find(PR(iAz,:) == 0); % remove zeros
            nzero = length(idx);
            PR(iAz,idx) = ones(length(idx),1).*NaN;
            z(iAz,:) = (PR(iAz,:)*2/1000).*cos(beta(iAz,:)); % measured distance from apex A to bed
            zmed = gmedian(z(iAz,:));
            %zcut = zmed + ztrim;
            %idx = find(z(iAz,:) > zcut); % remove outliers
            %ncut = length(idx);
            %z(iAz,idx) = ones(length(idx),1).*NaN;           
        end

        % extent of plot is equal to the width of the pencil sweep
        % The radius of the swept area = tan(SectorWidth/2)*head height
        %sr = tan((header.SectorWidth/2)*(pi/180))*header.Height;         
        dxy = 0.01;    % Key setting...determines image resolution at cost of
        rot = 0; % rotation, ifneeded
        sr = max(max(x));
        Xplot = -sr:dxy:sr;
        Yplot = Xplot;
        [Xi, Yi] = meshgrid(Xplot,Yplot);
        [thi, ri] = cart2pol(Xi, Yi);
        % convert radians to degrees, then add the rotation offset
        thi = (thi/pi)*180;  %thi = thi*180/pi;
        thi = thi+rot;
        thi(thi>180)=thi(thi>180)-360;
        thi(thi< -180)=thi(thi< -180)+360;
        
        %[theta, rho] = cart2pol(x,y);
        %Zi = interp2(rho, theta, z, ri, thi);
        %imagesc(Xplot,Yplot,Zi);
        

    otherwise % no plot
end