function [retnx,retny,retnz]=plotrange_cdf(ncname, tidx, settings)
%plotrange_cdf reads azimuth _raw.cdf, and plots one of 4 types of image
%
% usage [retnx,retny,retnz]=plotrange_cdf('az2009-02-02_raw.cdf, 1, settings)
%  first argument is filename to open  (expects raw azimuth .cdf)
%  tidx can be any time index present in the file
%  third argument is a structure with these fields
%    Pencil_tilt- angle from vertical of pencil beam (-2.248 is right for
%        Hatteras09)
%    plottype- can be : '2d' | '3draw' | '3d' | '3d_frm_img' | 'scat_frm_img'
%    2d - plots profile_range from each individual scan in the series
%    3draw - uses profile_range as source of depth, and plots with no interpolation
%    3d - uses profile_range as source of depth, and griddats to interpolate (sparse result)
%    3d_frm_img - uses data from the raw_image and does griddata to fill space
%    scat_frm_img - uses the raw image data and uses scatter to plot
%     (this is my preference for really seeing what's there)
%    dxy is needed to size the output from griddata- .05 is good for testing
%     (only needed for 3d and 3d_frm_img)
%    rot2compass - the value needed to get image plotted so that +y is north
%     this value should be the pen_adcp_offset -90(since scan 1 is perpencicular
%     to pencil 0)-90 (for math to cartesian) +declination,
%     so for Hatteras09 should be 182-90+10
%   returns x,y,z points used to make the plot displayed.
%
%  original by MM July07
%  EM mods to read az*cdf June 09
%
% note- autonan MUST be on for 3d_frm_image and scat_frm_img to work
%
if nargin < 3;
    help (mfilename); return
end

if ~exist(ncname,'file')
    [filename, pathname, filterindex] = uigetfile('*.*','Open an az*CDF file');
    if ~filename, return; end
    ncname = fullfile(pathname, filename);
end

%Open the data file
nc=netcdf(ncname);

if isempty(nc.Height(:))
    nc.Height = 1;   % should be read from header, but isn't in there now.
end

if ~isfield(settings, 'plottype'),
    plottype = '3draw'; % default
else
    plottype=settings.plottype;
end
if ~isfield(settings, 'Pencil_tilt'),
    tilt = 0; % default
else
    tilt=settings.Pencil_tilt;
end
if ~isfield(settings, 'rot2compass'),
    settings.rot2compass=0;
end
%clf
switch plottype
    %     case '2draw'  %doesn't work as expected
    %         % trim
    %         ntrim = 57;
    %         imagedata = squeeze(nc{'raw_image'}(idx,:,:,:));
    %         A = nc.Height(:); % height of pencil sweep vertex above nominal sea bed
    %         beta = squeeze(nc{'headangle'}(idx,1,1:end-ntrim)); % 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, nScans] = size(imagedata);
    %         for iAz = 1:nAz,
    %             nidx = find(imagedata(iAz,:,:) == 0);
    %             imagedata(iAz,nidx) = ones(length(nidx),1).*NaN;
    %             plot(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',nc{'azangle'}(iAz)))
    %             pause(2)
    %             cla
    %         end
    case '2d'           % this works well
        % trim
        ntrim = 57;
        ztrim = 5; % in cm
        if nc.Range(:) < 5,
            factor = 0.002;
        else
            factor = 1;
        end
        PR = squeeze(nc{'profile_range'}(tidx,:,1:end-ntrim)).*factor;
        A = nc.Height(:); % height of pencil sweep vertex above nominal sea bed
        beta = squeeze(nc{'headangle'}(tidx,:,1:end-ntrim))+tilt; % the Pencil Head Angle (in degrees)
        Ro=0.1;         % 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',[nc.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 .6 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, nc{'azangle'}(iAz), zmed))
            pause(.2)
            cla
        end
        close(nc);
        retnx=[];retny=[];retnz=[];
        
    case '3draw'
        % trim
        ntrim = 57;
        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.1; %meters, pencil sweep apex offset from azimuth centerline of rotation
        if nc.Range(:) < 5,
            factor = 0.002;
        else
            factor = 1;
        end
        PR = squeeze(nc{'profile_range'}(tidx,:,1:end-ntrim)).*factor;
        x = ones(size(PR)).*NaN;
        y = x;
        z = x;
        A = nc.Height(:); % height of pencil sweep vertex above nominal sea bed
        beta = squeeze(nc{'headangle'}(tidx,:,1:end-ntrim))+tilt; % the Pencil Head Angle (in degrees)
        beta = beta.*(pi/180); % convert to radians
        m = A.*tan(beta); % horizontal distance from sweep apex to measurement M
        alpha = nc{'azangle'}(tidx,:)+settings.rot2compass; % 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));
            %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;
        axis([-2.5 2.5 -2.5 2.5 -2 0])
        caxis([-1.6 0])
        view (0,90)
        colorbar
        close(nc);
        retnx=x; retny=y; retnz=z;
        %surfc(x,y,z);
        
    case '3d'
        % trim
        ntrim = 57;
        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 nc.Range(:) < 5,
            factor = 0.002;
        else
            factor = 1;
        end
        PR = squeeze(nc{'profile_range'}(tidx,:,1:end-ntrim)).*factor;
        x = ones(size(PR)).*NaN;
        y = x;
        z = x;
        A = nc.Height(:); % nominal height of pencil sweep vertex above sea bed
        beta = squeeze(nc{'headangle'}(tidx,:,1:end-ntrim))+tilt; % the Pencil Head Angle (in degrees)
        beta = beta.*(pi/180); % convert to radians
        m = A.*tan(beta); % horizontal distance from sweep apex to measurement M
        alpha = nc{'azangle'}(tidx,:)+settings.rot2compass; % 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
        if (isfield(settings,'dxy'))
            dxy = settings.dxy;
        else
            dxy=.02;
        end
        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')
        close(nc);
        retnx=Xplot;retny= Yplot; retnz=Zi;
        
        
        % uses the line on the raw image to approximate profile_range
    case 'scat_frm_img'
        % trim
        ntrim=57;
        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
        % if Ro is 0.0, there's a nan in the middle of the computed grids
        % which messes up griddata- use a small positive value
        Ro = 0.05; %meters, pencil sweep apex offset from azimuth centerline of rotation
        beta = squeeze(nc{'headangle'}(tidx,:,:)+tilt); % the Pencil Head Angle (in degrees)
        [nt naz nang]=size(beta);
        if find(isnan(beta(1,:)))
            nang=find(isnan(beta(1,:)),1,'first');
        end
        beta=beta(:,1:nang-1);
        if isempty(beta)
            disp('AutoNan needs to be on for this to work- please set it and try again')
            close(nc);
            return
        end
        
        x = ones(size(beta)).*NaN;
        y = x;
        z = x;
        A = nc.Height(:); % nominal height of pencil sweep vertex above sea bed
        beta = beta.*(pi/180); % convert to radians
        m = A.*tan(beta); % horizontal distance from sweep apex to measurement M
        alpha = nc{'azangle'}(tidx,:)+settings.rot2compass; % 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(beta);
        
        for iAz = 1:nAz,
            npoints = nc.DataPoints(:);
            range_config=nc.Range(:);
            SampPerMeter = npoints/range_config;
            raw_img=squeeze(nc{'raw_image'}(tidx,iAz,:,1:nang-1));
            %Find the equivalent of profile_range
            mx=max(raw_img(50:end-1,:));
            for jj=1:length(mx)
                mtch=find(mx(jj)==raw_img(50:end,jj));
                if isempty(mtch)
                    locs(jj)=0;
                else
                    locs(jj)=mtch(1);
                end
            end
            newPR=(locs+50)*nc.Range(:);
            % 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
            idx = find(newPR(1:end-ntrim) == 0); % remove zeros
            nzero = length(idx);
            newPR(idx) = ones(length(idx),1).*NaN;
            z(iAz,1:nang-1) = (newPR(1:nang-1)*2/1000).*cos(beta(iAz,1:nang-1)); % measured distance from apex A to bed
            % z is the elevation based on strongest return of each bin
        end
        %trim all to remove extraneous last line
        xn=x(:,1:nang-1);  %minus x needed to get the correct right/left
        yn=y(:,1:nang-1);
        % z has opposite + and - direction for each scan, so have to flip
        zn=z(:,1:nang-1);
        scatter(xn(:),yn(:),3,-(zn(:)),'filled')
        axis ([-2.5 2.5 -2.5 2.5])
        caxis([-1.5 -.3])
        colorbar
        xlabel('meters')
        ylabel('meters')
        title('scatter plot of individual sweeps from images')
        close(nc);
        retnx=xn; retny=yn; retnz=zn;
        
    case '3d_frm_img'
        % trim
        ntrim=57;
        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
        % if Ro is 0.0, there's a nan in the middle of the computed grids
        % which messes up griddata- use a small positive value
        Ro = 0.05; %meters, pencil sweep apex offset from azimuth centerline of rotation
        beta = squeeze(nc{'headangle'}(tidx,:,:)+tilt); % the Pencil Head Angle (in degrees)
        [nt naz nang]=size(beta);
        if find(isnan(beta(1,:)))
            nang=find(isnan(beta(1,:)),1,'first');
        end
        beta=beta(:,1:nang-1);
        if isempty(beta)
            disp('AutoNan needs to be on for this to work- please set it and try again')
            close(nc);
            return
        end
        
        x = ones(size(beta)).*NaN;
        y = x;
        z = x;
        A = nc.Height(:); % nominal height of pencil sweep vertex above sea bed
        beta = beta.*(pi/180); % convert to radians
        m = A.*tan(beta); % horizontal distance from sweep apex to measurement M
        alpha = nc{'azangle'}(tidx,:)+settings.rot2compass; % 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(beta);
        
        for iAz = 1:nAz,
            npoints = nc.DataPoints(:);
            range_config=nc.Range(:);
            SampPerMeter = npoints/range_config;
            raw_img=squeeze(nc{'raw_image'}(tidx,iAz,:,1:nang));
            %Find the equivalent of profile_range
            mx=max(raw_img(50:end,:));
            for jj=1:length(mx)
                mtch=find(mx(jj)==raw_img(50:end,jj));
                if isempty(mtch)
                    locs(jj)=0;
                else
                    locs(jj)=mtch(1);
                end
            end
            newPR=(locs+50)*nc.Range(:);
            % 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));
            % used to debug rotation question
            %if iAz < 5
            %    [alpha(1:10) x(iAz,1:10)' y(iAz,1:10)']
            %end
            % z = distance of sea bed surface from point A
            z(iAz,1:nang-1) = (newPR(1:nang-1)*2/1000).*cos(beta(iAz,1:nang-1)); % measured distance from apex A to bed
        end
        if (isfield(settings,'dxy'))
            dxy = settings.dxy;
        else
            dxy=.02;
        end
        sr = max(max(x));
        %Xplot = min(min(x)):dxy:max(max(x));
        Xplot = -sr:dxy:sr;
        Yplot = Xplot;
        [Xi, Yi] = meshgrid(Xplot,Yplot);
        % then compute Z, which interpolates in some areas and fills others
        Zi = griddata(x,y,z, Xi,Yi);
        %imagesc(Xplot,Yplot,Zi);
        pcolor(Xplot,Yplot,-Zi); shading flat;
        axis square
        colorbar
        axis ([-2.5 2.5 -2.5 2.5])
        xlabel('meters')
        ylabel('meters')
        caxis([-1.5 -.3])
        title('image interpolated onto a regular grid: nan for median +/- 1m')
        close(nc);
        retnx=Xplot; retny=Yplot; retnz=Zi;
        
    otherwise % no plot
        disp('The plottype entered didn''t match the options')
        disp('please use one of: 2d 3d 3draw 3d_frm_img scat_frm_img')
        close(nc)
        hold off
        ret_mat=[];
end
hold off