function f = pfan( settings )
% pfan - Process Imagenex fan sonar data from the raw netCDF file
%
%  ** This is a dead-end experiment!!  use procfan instead!!
%
% Input:
%  settings.fname - name of netcdf file with raw sonar data, e.g. '
%  settings.oname - name of output netcdf file (if desired)
%  settings.metaFile: string with name of metadata file [default MMMfan_raw.nc]
%  settings.img_nums: vector with image indices to process [default: all]
%  settings.make_nc
%  settings.make_plot
%  settings.sweep
%  settings.dxy - determines image resolution at cost of speed and storage
%          (reasonable range 0.02 to 0.005)
%  settings.fangle =
%      [ jtime angle;...
%        jtime angle;...
%        ...
%        jtime angle]; % array of times and associated azimuths
%  OR
%  settings.fangle = nnn; % constant azimuth to be used entire time
%
%      Note: If not present or has any NaNs, assume constant angle = 0.
%        If only one time and one angle (dimension [1 2]), use that angle
%        for all images.
%        Otherwise, interpolate angle from settings.fangle for image time.
%        User should ensure this results in reasonable values, esp. when
%        tripod is bumped or angle crosses north.
%
%
%           fname is the netcdf file containing the raw data.  The
%                     rootname will be used to create the name of the
%                     processed file
%           img_nums is the array of image indices to process
%                     can use [1 10 25] or [132:181],
%                     default is to process all
%              ** nb: if you choose discontinuous elements, they will be
%                     put into sequential elements in the output file, so
%                     the timebase is likely to be irregular.
%           ncp is the processed output netcdf object
%
% csherwood@usgs.gov
% emontgomery@usgs.gov
% Based on DO_FAN_ROTS

% Dependencies:
%   USGS NetCDF Toolbox (C. Denham)
%   -dap enabled mexnc doesn't work! ==> mexnc_win_2006a\mexnc.mexw32 works
%   p2xyfan
%   definefanprocncfile.m

% SVN revision number
rev_info= 'SVN $Revision: 982 $';

% Check for raw sonar netcdf file
if(exist(settings.fname,'file')~=2),
   error('Cannot find sonar file %s\n',settings.fname);
end

% open existing cdf file of raw fan data
ncr=netcdf(settings.fname);
% set up how many images to process
if(isfield(settings,'img_nums')),
   nimg_nums = length(settings.img_nums);
else
   nimg_nums=[1:1:length(ncr{'time'})];
end

% decide how big the xy matrices need to be
x=[-ncr.Range(:):settings.dxy(:):ncr.Range(:)];
dim_nc.x=length(x);
dim_nc.y=length(x);
dim_nc.sweep=settings.sweep;

if(settings.make_nc),
   % set up netcdf output file
   if isfield(settings,'oname')
      ofproc=settings.oname;
   else
      ofproc='sonar_proc.nc'
   end
   % instantiate the output ncfile
   ncp = definefanprocnc(ofproc, settings, dim_nc);
   % copy attributes from raw file
   rawAtts=ncnames(att(ncr));
   for ik=1:length(rawAtts)
      eval(['ncp.' char(rawAtts(ik)) '= ncr.' char(rawAtts(ik)) '(:);'])
   end
   % if there's information in settings, replace the ncp attributes
   % with the values in settings
   nn=fieldnames(settings);
   for ik = 1:length(nn)
      eval(['ncp.' nn{ik} '(:)=settings.' nn{ik} ';'])
   end
   % since StepSize is wrong in header, add degreesPerStep here
   ncp.DegPerStep=ncr.StepSize(:);
   ncp.DegPerStep(:)= ncr{'headangle'}(5)-ncr{'headangle'}(4);
   %reset creation date to now
   ncp.CREATION_DATE = ncchar(datestr(now));

   % do the right number of time elements
   ncp{'time'}(1:length(nimg_nums))=ncr{'time'}(nimg_nums);
   ncp{'time2'}(1:length(nimg_nums))=ncr{'time2'}(nimg_nums);
   % put the outputs into processed netcdf file
   for kj=1:settings.sweep
      ncp{'sweep'}(kj)=kj;
   end
   ncp{'sonar_image'}.scale_factor(:)=10000;
end

Fanidx = 1;
fsettings.ncr = ncr;
fsettings.plottype = 'polar';
fsettings.dxy = settings.dxy;
for jj=(nimg_nums(1):nimg_nums(end))
   % process the images and put into output netcdf file
   fsettings.ik = jj;
   jt= ncr{'time'}(jj)+ncr{'time2'}(jj)./86400000;
   fprintf(1,'Processing %s\n',datestr(gregorian(jt),'HHMM dd-mmm-yy'));
   fsettings.rot = interp1(settings.fangle(:,1),settings.fangle(:,2),jt);
   fprintf(1,'Rot = %f\n',fsettings.rot);
   fxy=p2xyfan(fsettings);
   Xplot = fxy.x;
   Yplot = fxy.y;
   Zs = fxy.zs;
   % and put what's returned in the output file and object
   if(settings.make_nc)
      if Fanidx==1
         ncp{'x'}(1:length(Xplot))=Xplot;
         ncp{'y'}(1:length(Yplot))=Yplot;
      end
      % Zs is float- needs to be multiplied by 10000 to store as short
      for kk=1:settings.sweep
         tmp1=floor(Zs(kk,:,:)*10000);
         ncp{'sonar_image'}(Fanidx,kk,1:length(Xplot),1:length(Yplot))=squeeze(tmp1);
         clear tmp1
      end
   end
   Fanidx=Fanidx+1;
end
% done with raw sonar data
close(ncr);

if(settings.make_nc)
   % add to history & make some notes to the netcdf file
   hist = ncp.history(:);
   hist_new = ['Sonar processed with ' ,mfilename, ', ', rev_info, ', using Matlab ' ,...
      version, '; ',hist];
   ncp.history = hist_new;

   ncp.NOTE =['radial data interpolated onto x-y grid to make image;',...
      'image rotated so that +y (up) is N'];
   ncp.NOTE1 = ['To view images in Matlab type the following at the command ',...
      'prompt:  nc=netdcf(''sonarxxx.nc'');',...
      'imagesc(nc{''x''}(:),nc{''y''}(:),squeeze(nc{''sonar_image''}(n,p,:,:)));',...
      'set(gca,''ydir'',''normal''); **where n & p are the time and sweep indexes'];

   % this is where the data is saved
   t1= ncp{'time'}(1)+ncp{'time2'}(1)./86400000;
   ncp.start_time = datestr(gregorian(t1));
   tt= ncp{'time'}(end)+ncp{'time2'}(end)./86400000;
   ncp.stop_time = datestr(gregorian(tt));
   t_all= ncp{'time'}(:)+ncp{'time2'}(:)./86400000;
   % if time data is evenly spaced
   if length(t_all) > 1 & isempty(find(diff(diff(t_all))) ~= 0)
      ncp.DELTA_T = [num2str(gmean(diff(t_all))*24*60),' sec'];
      % time and time2 are EVEN by default
   else
      ncp{'time'}.type(:)='UNEVEN';
      ncp{'time2'}.type(:)='UNEVEN';
      ncp.DELTA_T = ['? sec'];
   end
   % close the writeable version
   close(ncp);
   ncclose;
   % re-open it read-only to return to matlab (why? - CRS)
   % eval(['ncp=netcdf(''' ofproc '.cdf'');'])
end
f = 'done';