function [Xplot, Yplot, thi, ri, Zs, rot] = procfan(ncr, ik, plottype, settings);
% PROCFAN - rotate Imagenex fan sonar data from _raw.cdf and output images
%  called by do_fan_rots.m - follows execution of mk_rawcdf_08
% [x,y,thi, ri, zs] = procfan(nc, ik, plottype, settings)
%
% reads netcdf file created by mkrawcdf
% nc{'imagedata'}(time,points,scan) image data with header removed
% nc{'head_pos'}(time,scan) head position in counts
% nc{'headangle'}(scan) head rotation angle in degrees
% nc{'profile_range'}(time,scan) Imagenex computed range to surface
% nc{'nDataBytes'}(scan) = number of byte in the ping
% nc.HeadType. = IGX/IMX etc.
% nc.HeadID. = ID number
% nc.Range = range setting
% nc.NReturnBytes = number of bytes in the ping + header
% nc.DegPerStep = from FanSwitches.StepSize;
% nc.SectorWidth = from FanSwitches.SectorWidth;
%
% ik holds the indices of the scans to process
% plottype = 'square' | 'polar' | 'prange'
% still need to pass 'settings' to get or change these params:
%    Height = instrument installed height,
%    fanadcp_off = adcp compass value,
%    magnetic_variation= declination
%    dxy= size of grid= .01 is smallest resolution
%    mkplt= 'y' or 'n', whether to make plots
%   these should come from the metafile, but if supplied, will be used
%
%  originally, this program did the whole raw file together.  With high
%  resolution, it was impossible, so now does one image at a time
%    ik = time index to process
%
% outputs:
% Single matrices of Xplot and Yplot that data is interpolated into
% Single matrices of thi and rho, after rotation
% Zs matrix of rotated image (images, if more than 1 sweep)
% rot is the total angle the image was rotated by 
%   rot=fanadcp_off+adcp3val+(-magnetic_variation)
%
%  emontgomery@usgs.gov
%  April 2, 2008

if (nargin < 4); help mfilename; return; end

tmp = size(ncr{'raw_image'});
ntimes=tmp(1); NPoints=tmp(2); nscans=tmp(3);
% instead of using user supplied settings,  use what's read in the header
% fanbeam sometimes returns .6 for StepSize: we know it should be .15, thus the /4

% replace values in ncr used in plotting if there's a settings value to use
if (isfield(settings,'fanadcp_off') & isfield(settings,'adcp3val'))
    rval=settings.fanadcp_off+settings.adcp3val;
elseif isfield(settings,'fanadcp_off')
    rval=settings.fanadcp_off(:);
elseif isfield(settings,'adcp3val')
    rval=settings.adcp3val(:);
elseif ~isempty(ncr.fanadcp_off(:))
    rval=ncr.fanadcp_off(:);
else
    rval = 0;
end
if isfield(settings,'Height')
    hgt=settings.Height;
else
    hgt=ncr.Height(:);
end
if isfield(settings,'dxy')
    dxy=settings.dxy;
else
    dxy=ncr.dxy(:);
end
if isfield(settings,'magnetic_variation')
    mag_var=settings.magnetic_variation;
else
    if isstr(ncr.magnetic_variation(:))
     mag_var=str2num(ncr.magnetic_variation(:));
    else
     mag_var=ncr.magnetic_variation(:);
   end
end
if ~isfield(settings,'mkplt')
    mkplt='n';
else 
    mkplt='y';
end

% when more than one sweep, the first and middle scans are redundant (based
% on headangle- have to remove them later- the last scan should be good.
% trim the redundant first scan
imagedata = ncr{'raw_image'}(ik,:,2:nscans);
%  headangle=HeaderData(:).HeadAngle(2:nscans-1);
prange = ncr{'profile_range'}(ik,2:nscans);
headangle=ncr{'headangle'}(2:nscans); %remove duplicate first point
mnStepSize=mean(diff(headangle(10:20)));
[npoints,nangles] = size(imagedata);  % necessary to
slantrange = ncr.Range(:)/NPoints:(ncr.Range(:)/NPoints): ncr.Range(:);

sectorswept = 1:mnStepSize:ncr.SectorWidth(:)+(mnStepSize)+1;
sectorswept = sectorswept(2:length(sectorswept)-1);
% trim range to real area of data
nearfieldcutoff = 52;
farfieldcutoff = NPoints;
imagedata = imagedata(nearfieldcutoff:farfieldcutoff,:);
slantrange = slantrange(nearfieldcutoff:farfieldcutoff);

% other calculations
%horizrange = slantrange.*acos(asin((settings.Height./slantrange)));
%horizrange = settings.Height.*cos(asin(slantrange./settings.Height));
horizrange = real((slantrange.^2 - hgt^2).^(1/2));

% set up the interplolation grid
Xplot = (- ncr.Range(:):dxy:ncr.Range(:));
Yplot = Xplot;
% pre-allocate Zs
Zs=ones(ncr.sweep(:),length(Xplot), length(Yplot));

%etm fix 3/27/07
%   ** simply flipping the image puts the tray on the correct side.
%    imagedata=flipud(imagedata);
% but making [ripad and thipad] from positive to negative does the same
% thing. I guess it depends on what sweep direction we believe happened
%
n=0;
for kk=1:ncr.sweeps(:)
    % may want to only do the first sweep because second doesn't seem to
    % overlay correctly
    % etm 11/20/09 --for two sweeps, indexing for the two sweeps should be
    %  sweep 1 = 2:length(headangle)/2
    %  sweep 2 =  length(headangle)/2+2:end
    % starting at 2 and adding 2 in the middle is IMPORTANT

    if ncr.sweeps(:) > 1
        if kk==1  % this only does 2 sweeps now
           %imagedata already has first scan cut off, so imgadata is indexed +1 from
           % raw_data here.  Still have to avoid the middle point
           imdata=imagedata(:,1:floor(length(imagedata)/2));
           ha=ncr{'headangle'}(2:length(sectorswept)+2);
        elseif kk == 2
            % imagedata(2323:4643) is same as ncr{'raw_image'}(2424:4644)
            imdata=imagedata(:,floor(length(imagedata)/2)+2:end);
            ha=ncr{'headangle'}(length(sectorswept)+4:end);
        end
    else
        imdata=imagedata(:,2:end);
        ha=headangle;
    end
    % have to flip the even sweeps because interp2 requires monotonic increasing ha
    if (mod(kk,2)==0)
        imdata=fliplr(imdata);
        %ha must be flipped to get monotonic increase for the interp
        % first sweep is always - to +, 2nd, 4th, etc are + to -
        ha=flipud(ha);
        % in back sweep ha never gets to -174.  To get sweep 1 & 2 images to
        % overlay after interp, have to decrease all by .15
        % ha=ha+1.05;   % I think there's about a 7 step lag
        %ha=ha-mnStepSize;
    end
    if (ha(1) == ha(2))
        ha(1)=ha(2)-mnStepSize;
    end
    if (ha(end)== ha(end-1))
        ha(end)=ha(end)+mnStepSize;
    end

    % now get the new size for plotting
    [npoints,nangles] = size(imdata);

    if exist('plottype') == 1,
        switch plottype
            case 'square'
                % plot linear
                set(gcf,'name','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,' Fan beam image taken ',...
                    datestr(datenum(gdate))],...
                    'units','normalized','color','w')
            case 'polar'
                clear thi ri Xi Yi sr fakeheadangle

                settings.fillval = 255; % sets value of no-data region to make
                % a white background

                %  new way : have this program compute rot rather than have the user!!!
                %    -90 is needed to go from math convention where 0 is horizontal +
                %    to map or geographic wher 0 is UP.
                %    rot = -90 + ((settings.adcp3val-360) + settings.magnetic_variation + settings.fanadcp_off);
                %    rot=0       % apply NO rotation to display raw image.
                rot= rval+(-mag_var);  % should be opposite sign since matlab + angle is cartisian !
                %    convention is ccw is + from 0 at x
                if rot > 360; rot=rot-360; end
                if rot < -360; rot=rot+360; end


                [Xi, Yi] = meshgrid(Xplot,Yplot);
                % [thi, ri]=cart2pol(Xi,Yi);
                % with outputs cartesian coords, so don't have to fuss with math coords
                [ri, thi]=pcoord(Xi, Yi);   %  ri is -180-> 180 cw, with 0 up

                % convert radians to degrees, then add the rotation offset
                % which needs to be negative to be done in math coords
                %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;

                % pad ranges and angle of sweep
                % headangle(1)=-174; headangle(nangles)=161.25
                % fakeheadangle = linspace(headangle(1), headangle(nangles), ...
                %		     nangles)';
                sr = [0,slantrange(1,1)-.1,...
                    slantrange(1,:), slantrange(1,npoints)+.1, 9]';
                % pad angles to get complete circle
                % first 3 of first image are always the same
                %if kk == 1
                %xx=[ha(3):mnStepSize:ha(nangles-1)+mnStepSize];
                %hfill=[ha(3)-(2*mnStepSize) ha(3)-mnStepSize];
                fh=[ha(1:nangles)];
                % I don't know how the diff(diff(fh)) condition came from, but mostly we
                % don't want to use it- else is better.  This was needed
                % for one of the older datasets, so is removed for the
                % moment.  Don't use that confition, if the section is
                % replaced in the code
                %if (diff(diff(fh)) ~= median(diff(fh)))
                %    fakeheadangle=[ha(1):mnStepSize:ha(nangles-1)+mnStepSize]';
                %    % there was also often duplication of the first 3
                %    % headangles, which may account for problems with the size
                %    % of the computed fakeheadangle, so add angles to the end if
                %    % needed
                %    if length(fakeheadangle)< nangles
                %        nmiss=nangles-length(fakeheadangle);
                %        for ik=1:nmiss
                %            fakeheadangle(end+ik)=fakeheadangle(end)+(ik+1)*mnStepSize;
                %        end
                %    end
                %    lopad = [-179.999 min(fakeheadangle)-mnStepSize];
                %    hipad = [max(fakeheadangle)+mnStepSize 179.99];  % was 180
                %    nlo = length(lopad);
                %    nhi = length(hipad);
                %    fakeheadangle = [lopad fakeheadangle' hipad];
                %else
                    fakeheadangle=[ha(1:nangles)];
                    % pad for the dead zone area
                    lopad = [-179.999 min(ha)-mnStepSize];
                    hipad = [max(ha)+mnStepSize 179.99];  %was 180
                    nlo = length(lopad);
                    nhi = length(hipad);
                    fakeheadangle = [lopad fakeheadangle' hipad];
               % end
                                             
                % pad data
                Zpad = (settings.fillval)*ones(npoints+nlo+nhi,nangles+nlo+nhi);
                %Zpad = (settings.fillval)*ones(npoints+4,nangles);
                % if we let the condition by <0, it does a better job in terms of look
                % but the log10 that follows makes those 0's -inf, which makes writing
                % into the short later barf.
                imdata(imdata<(0))=(settings.fillval);
                % set the 0's == 1 to allow the log10 to operate correctly and still
                % get the effect of having missing values replaced by a dark pixel
                % instead of a light one
                imdata(imdata==0)=1;

                % kludge...last row (449) of range values always seems to be 252
                imdata(449,:)=(settings.fillval);
                Zpad(nlo+1:npoints+nhi,nlo+1:nangles+nhi)=imdata;
                % the log10 is needed to more evenly distribute the data in bins
                % omitting it leads to a very dark image
                Zpad = log10(Zpad);

                % interpolate onto polar shape
                Zi = interp2(fakeheadangle,sr,Zpad, thi,ri);
                %    Zi(Zi<1)=(8); % produces weird contrasty plot
                if ~isempty(strfind(mkplt,'y'))
                    imagesc(Xplot,Yplot,Zi)
                    set(gca,'ydir','normal') % because imagesc is sidedownup
                    colormap gray;
                    set(gca,'Fontsize',12);
                    set(gca,'Xtick',[-5:1:5])
                    set(gca,'Ytick',[-5:1:5])
                    
                    if(1), % plot range circles
                        hold on
                        for icirc = 1:5,
                            h = circle(icirc,0,0);
                            set(h,'color',[.7 .7 .6])
                        end
                    end
                    % add something to indicate which sweep
                    if kk==1
                        swtext='sweep 1-';
                        fh=ones(length(fakeheadangle),2);
                    else
                        swtext='sweep 2-';
                    end
                    axis square
                    text(0.01,0.03,[swtext 'rotation used = ' num2str(rot)],...
                        'units','normalized','color','r','horizontalalignment','left','fontsize',10)
                    % text(0.8,0.95,'\uparrow North',...
                    %     'units','normalized','color','y','fontsize',10)
                    % until the time thing in 2009a gets fixed, forced
                    % conversion to double
                    tt= double(ncr{'time'}(ik))+double(ncr{'time2'}(ik))./86400000;
                    ts = datestr(gregorian(tt));
                    text(0.99,0.03,ts,...
                        'units','normalized','color','r','horizontalalignment','right','fontsize',10)
                    fhd(:,kk)=fakeheadangle';  % sometimes needs transpose
                    hold off
                    title ('rotation =(initial adcp compass + fan\_adcp\_offset + -magnetic variation)')
                    %                 gtext([num2str(settings.adcp3val) ' + ' num2str(settings.fanadcp_off) ' + -(' ncr.magnetic_variation(:) ')'])
                end
                clear imdata sr ha fakeheadangle xx hfill Zpad

            case 'glacial' % this is really slow
                set(gcf,'name','Polar')
                set(gcf,'numbertitle','off')
                [th, r] = meshgrid(sectorswept.*(pi/180), horizrange);
                [X, Y] = pol2cart(th, r);
                %Z = imdata';
                Z = log10(imdata+eps)';
                surf(X, Y, Z)
                view([0 90]); % az, el
                shading flat
                colormap jet
                text(0.01,0.05,[fname,' Fan beam image taken ',datestr(datenum(gdate))],...
                    'units','normalized','color','k','fontsize',9)
                %axis equal tight
                %set(gca,'xlim',[0 settings.Height+0.5]);

            case 'prange'
                % plot headposition and prange
                [ax, h1, h2] = plotyy(1:nangles,headangle,1:nangles,prange);
                title([fname,' Fan beam image taken ',datestr(datenum(gdate))])
                xlabel('Head position and profile range by ping')
                ylabel('Deg')
                %set(ax(2),'ylabel',text(0,0,'string','Prange'))
            otherwise % no plot
                %return values
                Xplot = slantrange;
                Yplot = headangle;
                Zi = imdata;
                % end  %end for kk=
        end
        % put Zi into a structure, so both images may be accessed
        %  Zs{ik,kk}=Zi;
        Zs(kk,:,:)=Zi;
    end
    % let the user know where we are
    if mod(ik,10) == 0
        disp(['processing record ' num2str(ik) ' of ' num2str(length(ncr{'time'}))])
    end
    clear Xi Yi sr fakeheadangle
end
end
% keep3 imdata Xplot Yplot Zi FanData settings