function animatesonar2nc(metaFile, outFileRoot)
% animsonar1023.m  A function to evaluate Imagenex sonar data %
%
% usage:
%
%       where:  metaFile is the name of your text file containing metadata,
%                        surrounded by single quotes WITHOUT the file
%                        extension .txt. An example metadata file,
%                        sonarmetaexample.txt, is provided in this package of
%                        mfiles.
%               outFileRoot is the name given to the output netCDF file,
%                           surrounded by single quotes WITHOUT the file
%                           extension .cdf.  This file be named:
%                             'outFileRoot'im-cal.nc
%
% Written by Charlene Sullivan
% USGS Woods Hole Field Center
% csullivan@usgs.gov
%
% Dependencies:
%   USGS NetCDF Toolbox (C. Denham)
%   smart_interp.m (P. Alexander)
%   showfan.m (M. Martini, C. Sherwood)

% C. Sherwood 23 June 2006

% C. Sullivan   10/03/05,   version 1.0
%This function is heavily based on the script animatesonar.m by M. Martini
%(thanks!). Use of this function does not require users to manually specify
%metadata information, copy sonar data directories to a root directory, and
%delete sonar data files that are less than 5 kb (as required by
%animatesonar.m). To use this function, users must include required metadata
%in a text file with specific formatting (see sonarmetaexample.txt in this
%package of m-files). This text file must reside in the current directory.
%Users can show currents and/or wave data in conjunction with sonar images
%by specifying the complete paths to the netcdf files containing this data
%in their metadata file. Sonar images are saved as a series of .png's in a
%subdirectory named 'frame' of the current directory. The naming convention
%used w/ the .pngs's is yyyymmddTHHMMSS.png, where 'yyyy','mm','dd','HH',
%'MM','SS', are the year, month, day, hour, minute, and second at which the
%image was colleccted.

% Notes from M. Martini on animating images from animatesonar.m:
% -- frame images are written to the local directory and also saved to
% Fanmovie or Pencilmovie in the workspace, so if you like MATLAB movies,
% save the workspace when this script is done.
% -- the best results so far have been by digitizing frames using VideoMach
% Microsoft Video 1 codec, 75% compression quality
% turn keep duration off
% turn automatic off
% 5 frames per second

%TODO: The pencil - how do we determine the files are not duds?

close all
more off

version = '1.1'; % Version updated by CRS June 2006

% Check for metadata file
metaPath = pwd;
meta = dir([metaFile,'.txt']);
if isempty(meta)
   fprintf('\n')
   fprintf('The metadata file %s.txt does not exist in this directory\n',metaFile)
   metaPath = input('Please enter the full path to the directory with your metadata file:  ','s');
   meta = dir(fullfile(metaPath,[metaFile,'.txt']));
   if isempty(meta)
      error('Still cannot find the metadata file ',fullfile(metaPath,[metaFile,'.txt']))
   end
end
metaFile = fullfile(metaPath,meta.name);

% Get user's metadata structure
settings = readSonarMeta(metaFile);

% Check that the metadata contains required fields.  If a required field
% is missing, ask the user for it.
reqFields = {'FirstSonarDay','LastSonarDay','WavesFile','RootDataDir',...
   'ADCPbin','CurrentsFile','SectorSweep','datayear',...
   'AngleSweepAround','StepSize','DataPoints','Height',...
   'SonartoAnimate'};
for f = 1:length(reqFields)
   if ~isfield(settings,reqFields{f})
      disp(['The field ''',reqFields{f},''' is not specified in ',metaFile,'.txt'])
      missingFields(f) = 1;
   else
      missingFields(f) = 0;
   end
end
if any(missingFields)
   disp('Required fields missing from the metadata');
   % ask user to check/change defaults
   % Default settings for Myrtle Beach fall deployment 2003-2004
   %settings.CurrentsFile = 'None specified';%'E:\CSULLIVAN\SOUTH_CAROLINA\DATAFILES\7321adc-a1h.nc';
   %settings.ADCPbin = 'None specified';%'7';
   %settings.WavesFile = 'E:\CSULLIVAN\SOUTH_CAROLINA\DATAFILES\7221adw-a.nc';
   %settings.FirstSonarDay = '01-Nov-2003';
   %settings.LastSonarDay = '21-Jan-2004';
   %settings.RootDataDir = 'E:\CSULLIVAN\SOUTH_CAROLINA\DATAFILES\732\7323im\movieTest';
   %settings.SonartoAnimate = 'Fan'; %'Pencil'
   %settings.SectorSweep = '348'; %'132'
   %settings.AngleSweepAround = '0';
   %settings.StepSize = '1';
   %settings.DataPoints = '500';  % number of 'ranges' recorded in each ping
   %settings.Height = '0.55'; %1.06
   %settings.datayear = '2003';
   %settings.rot = '112';
   % Default settings for Eurostrat first deployment 2002-2003
   settings.CurrentsFile = 'c:\home\data\MVCO'
    %pick a bin according to 0.42x total depth
    %mooring 701 in 9.2 m = 3.864 m depth
    %ADCPbin = 6; % at 3.5448 m
   settings.ADCPbin = '6';
   settings.FirstSonarDay = julian(2002,11,7,0);
   settings.LastSonarDay = julian(2002,12,8,0);
   settings.RootDataDir = 'C:\home\data\eutostrat\original_images\';
   settings.SonartoAnimate = 'Pencil';
   settings.SectorSweep = '132';
   settings.AngleSweepAround = '0';
   settings.StepSize = '1';
   settings.DataPoints = '500';  % number of 'ranges' recorded in each ping
   settings.Height = '1.06'; 
   settings.datayear = '2002';
   settings.rot = '112';
   defaultSonarSettings = {settings.CurrentsFile,...
      settings.ADCPbin,...
      settings.WavesFile,...
      settings.FirstSonarDay,...
      settings.LastSonarDay,...
      settings.RootDataDir,...
      settings.SonartoAnimate,...
      settings.SectorSweep,...
      settings.AngleSweepAround,...
      settings.StepSize,...
      settings.DataPoints,...
      settings.Height,...
      settings.datayear,...
      settings.rot};
   prompt  = {'Enter the complete path to the currents netcdf file:',...
      'Enter the ADCP bin to use:',...
      'Enter the complete path to the waves netcdf file:',...
      'Enter the date of the first sonar image:',...
      'Enter the date of the last sonar image:',...
      'Enter the complete path to the sonar data directories:',...
      'Enter the sonar you would like to animate:',...
      'Enter the sector that was swept:',...
      'Enter the center to sweep around:',...
      'Enter the Step size used:',...
      'Enter the data points setting (points/10):',...
      'Enter the sonar''s height off the bottom, m:',...
      'Enter the year the instrument was deployed:',...
      'Enter the rotation angle to rotate the image to true N, deg:'};
   dlgtitle   = 'Sonar input metadata for the setup';
   lineNo  = 1;
   dlgresult  = inputdlg(prompt,dlgtitle,lineNo,defaultSonarSettings);
   settings.CurrentsFile = str2num(dlgresult{1});
   settings.ADCPbin = str2num(dlgresult{2});
   settings.WavesFile = str2num(dlgresult{3});
   settings.FirstSonarDay = str2num(dlgresult{4});  % number of 'ranges' recorded in each ping
   settings.LastSonarDay = str2num(dlgresult{5});
   settings.RootDataDir = str2num(dlgresult{6});
   settings.SonartoAnimate = str2num(dlgresult{7});
   settings.SectorSweep = str2num(dlgresult{8});
   settings.AngleSweepAround = str2num(dlgresult{9});
   settings.StepSize = str2num(dlgresult{10});
   settings.DataPoints = str2num(dlgresult{11});
   settings.Height = str2num(dlgresult{12});
   settings.datayear = str2num(dlgresult{13});
   settings.rot = str2num(dlgresult{14});
end
clear reqFields missingFields

% Load currents and wave data
%  etm 08/06 turned off both currents and waves for playing with pencil
currents = 0;
if exist(settings.CurrentsFile) ~= 2,
   currents = 0;
end

waves = 0;
if exist(settings.WavesFile) ~= 2,
   waves = 0;
end

if currents,
   % load the currents during the desired time
   nc = netcdf(settings.CurrentsFile);
   timeobj = nc{'time'};
   timeobj = autonan(timeobj, 1);
   time2obj = nc{'time2'};
   time2obj = autonan(time2obj, 1);
   tj=timeobj(:)+time2obj(:)./(3600*1000*24);
   datenum_curr=datenum(gregorian(tj));
   % find the indeces to the dates of interest
   findidx = find(datenum_curr < datenum(settings.FirstSonarDay));
   if isempty(findidx), dateidx(1) = 1;
   else dateidx(1) = max(findidx);
   end
   if dateidx(1) < 0, dateidx(1) = 1; end
   findidx = find(datenum_curr > datenum(settings.LastSonarDay));
   if isempty(findidx), dateidx(2) = length(datenum_curr);
   else dateidx(2) = min(findidx);
   end
   % trim time
   tj = tj(dateidx(1):dateidx(2));
   datenum_curr = datenum_curr(dateidx(1):dateidx(2));
   % now get the data, replace fill values w/ nans, and linearly
   % interpolate over the nans
   ncvarnames = {'u_1205','v_1206','w_1204',};
   names = {'u','v','w'};
   for n = 1:length(ncvarnames);
      ncobj = nc{ncvarnames{n}};
      ncobj = autonan(ncobj,1);
      eval([names{n},' = ncobj(dateidx(1):dateidx(2),settings.ADCPbin);'])
      %         eval(['bads = find(isnan(',names{n},'));'])
      %         eval([names{n},'(bads) = nan;'])
      %         if ~isempty(bads)
      %             eval([names{n},' = smart_interp(tj, ',names{n},', tj, 1);'])
      %         end
   end
   bindepth = nc{'depth'}(:);
   waterdepth = nc.WATER_DEPTH(:);
   close(nc)
   clear nc ncobj tg tj timeobj time2obj
end

if waves,
   % load the waves during the desired time
   ncw = netcdf(settings.WavesFile);
   julian_wave = ncw{'time'}(:) + ncw{'time2'}(:)/3600/1000/24;
   datenum_wave = datenum(gregorian(julian_wave));
   % find the indeces to the dates of interest
   findidx = find(datenum_wave < datenum(settings.FirstSonarDay));
   if isempty(findidx), wavedateidx(1) = 1;
   else wavedateidx(1) = max(findidx);
   end
   if wavedateidx(1) < 0, wavedateidx(1) = 1; end
   findidx = find(datenum_wave > datenum(settings.LastSonarDay));
   if isempty(findidx), wavedateidx(2) = length(datenum_wave);
   else wavedateidx(2) = min(findidx);
   end
   % trim time
   julian_wave = julian_wave(wavedateidx(1):wavedateidx(2));
   datenum_wave = datenum_wave(wavedateidx(1):wavedateidx(2));
   % now get the data, replace fill values w/ nans,
   % linearly interpolate over the nans, remove leftover
   % nans at beginning and end of the timeseries
   ncvarnames = {'wh_4061','wp_4060','wvdir'};
   names = {'Hsig','Tm','wvdir'};
   for n = 1:length(ncvarnames);
      ncobj = ncw{ncvarnames{n}};
      eval([names{n},' = ncobj(wavedateidx(1):wavedateidx(2));'])
      eval(['bads = find(',names{n},'>=1e35);'])
      eval([names{n},'(bads) = nan;'])
      if ~isempty(bads)
         eval([names{n},' = smart_interp(julian_wave, ',names{n},', julian_wave, 3);'])
      end
      eval(['goods = find(~isnan(',names{n},'));'])
      eval([names{n},' = ',names{n},'(goods);'])
   end
   datenum_wave = datenum_wave(goods);
   julian_wave = julian_wave(goods);
   close(ncw);
end

% Get a listing of sonar data directories
SubDataDir = dir(settings.RootDataDir);

% Individual movie frames saved to a subdirectory
% of the current directory named 'frame'
outdir = 'frame';
[success,message,messageid] = mkdir(outdir);

clear FanTime PencilTime
Fanidx = 1;
Pencilidx = 1;
save settings settings
% Go through each sonar data directory and parse. The
% first two directories are always . and ..
for d=3:length(SubDataDir),

   disp(SubDataDir(d).name)

   % Get the year in which the data was collected
   [yyyy] = getSonarYear(SubDataDir(d).name, settings.FirstSonarDay,...
      settings.LastSonarDay);
   yyyy = str2num(yyyy);
   mm = str2num(SubDataDir(d).name(2:3));
   dd = str2num(SubDataDir(d).name(4:5));
   dirDay = julian(yyyy,mm,dd,0);
   fileDay = gregorian(dirDay - 1); %files are 1 day behind directory day
   settings.datayear = fileDay(1);

   % Act only on directories
   if exist(fullfile(settings.RootDataDir,SubDataDir(d).name)) == 7,
      % Get a listing of data files
      DataFiles = dir(fullfile(settings.RootDataDir,SubDataDir(d).name));
      for f = 3:length(DataFiles),

         sonarDataFile = fullfile(settings.RootDataDir,...
            SubDataDir(d).name,...
            DataFiles(f).name);
         disp(sonarDataFile)

         % The file name convention for sonar data files is
         % S1111020.F50 = S1MMDDHH.?MM where ? = F for fan
         % or P for pencil
         a = sscanf(DataFiles(f).name,'S1%c%c%c%c%c%c.%c%c%c')
         mon = str2num(a(1:2));
         day = str2num(a(3:4));
         hr = str2num(a(5:6));
         mn = str2num(a(8:9));
         datenum_sonar = datenum(settings.datayear,mon,day,hr,mn,0); %used to interpolate currents?

         % The file name convention for individual
         % sonar frames is yyyymmddTHHMMSS.png
            outf = ['B',datestr(datenum_sonar,30)];
         outfile = fullfile(outdir,outf);

         % Fan files less than 5 kb are duds and
         % should not be used in the animation
         [WhichSonar, validFile] = checkSonarFile(sonarDataFile);

         % For valid (>5 kb) fan files, display the sonar image
         % and write the image data to NetCDF
         if WhichSonar == 'F' && validFile == 1 && strcmpi(settings.SonartoAnimate,'fan'),
            % first, display the sonar image
            FanTime(Fanidx) = julian([settings.datayear mon day hr mn 0]);
            FanDatenum=datenum([settings.datayear mon day hr mn 0]);

            clearimfig = figure;
            set(gcf,'Position',[9 37 827 659])

            sonar_ax = axes('pos',[0.2 0.25 0.73 0.73]);
            [imagedata, Xplot, Yplot, Zi, settings] = ...
               showfan(fullfile(settings.RootDataDir,SubDataDir(d).name,...
               DataFiles(f).name),'polar',settings);
            ylabel('Fan Beam Sonar, range in meters')
            if currents,
               curridx1 = find(datenum_curr <= datenum_sonar);
               curridx2 = find(datenum_curr >= datenum_sonar);
               if ~isempty(curridx1) && ~isempty(curridx2)
                  if datenum_sonar-datenum_curr(curridx1(end))< datenum_curr(curridx2(1))-datenum_sonar
                     curridx = curridx1(end);
                  else
                     curridx = curridx2(1);
                  end
                  %convert u,v to polar coordinates
                  [currDir, currSpd] = uv2polar(u(curridx), v(curridx));
                  %convert current speed, direction to x, y. dividing
                  %speed by 10 to scale arrows for for 10 cm/s range rings
                  [currX, currY] = xycoord(currSpd/10,currDir);
                  %place the arrow on the sonar image
                  curr_arr = compass(currX, currY);
                  set(curr_arr,'linewidth',2,'color','r')
                  text(4,4,[num2str(round(currSpd)),' cm/s'],'color','r','horizontalalignment','center')
                  text(4,3.5,[num2str(round(currDir)),' deg'],'color','r','horizontalalignment','center')
               else
                  text(4,4,['--- cm/s'],'color','r','horizontalalignment','center')
                  text(4,3.5,['--- deg '],'color','r','horizontalalignment','center')
               end
            end
            if waves,
               waveidx1 = find(datenum_wave <= datenum_sonar);
               waveidx2 = find(datenum_wave >= datenum_sonar);
               if ~isempty(waveidx1) && ~isempty(waveidx2)
                  if abs(datenum_sonar-datenum_wave(waveidx1(end))) < abs(datenum_wave(waveidx2(1))-datenum_sonar)
                     waveidx = waveidx1(end);
                  else
                     waveidx = waveidx2(1);
                  end
                  %convert wave height (magnitude) and direction to x, y.
                  %(direction is direction FROM.  we want to see
                  %direction TO)
                  if wvdir(waveidx) <= 180
                     wvdir(waveidx) = wvdir(waveidx)+180;
                  else
                     wvdir(waveidx) = wvdir(waveidx)-180;
                  end
                  %calculate ub (wave orbital velocity) with Sherwood's ubqf.m
                  ub = ubqf( Tm(waveidx), Hsig(waveidx), waterdepth );
                  %amplitude of do (wave orbital diameter) is ub/omega:
                  omega = 2*pi/Tm(waveidx);
                  amp = ub / omega;
                  %ripple wavelength is 1.33*A
                  lambda = 1.33 * amp;
                  %place the arrow on the sonar image. scale lambda by 5
                  %so you see 5 wavelengths on the image
                  [waveX, waveY] = xycoord(lambda*5, wvdir(waveidx));
                  wave_arr = compass(waveX, waveY);
                  set(wave_arr,'linewidth',2,'color','c')
                  text(4,3,[num2str(round(lambda*500)/5),' cm'],'color','c','horizontalalignment','center')
                  text(4,2.5,[num2str(round(wvdir(waveidx))),' deg'],'color','c','horizontalalignment','center')
               else
                  text(4,3,['--- cm  '],'color','c','horizontalalignment','center')
                  text(4,2.5,['--- deg '],'color','c','horizontalalignment','center')
               end
               wave_ax = axes('pos',[0.2 0.055 0.73 0.13]);
               plot(datenum_wave,Hsig);
               hold on
               if ~isempty(waveidx1) && ~isempty(waveidx2)
                  wave_dot = plot(datenum_wave(waveidx),Hsig(waveidx),'r.','markersize',14);
               end
               yl=ylabel({'Significant','wave','height (m)'});
               set(yl,'fontsize',12)
               set(yl,'units','normalized')
               set(yl,'position',[-.07 .49 0])
               set(wave_ax,'fontsize',12)
               set(wave_ax,'xlim',[datenum(settings.FirstSonarDay) datenum(settings.LastSonarDay)])
               datetick('x',1,'keeplimits')
               set(wave_ax,'xticklabel',[])
               xt = get(wave_ax,'xtick');
               for t = 1:length(xt)
                  if xt(t)>=datenum(settings.FirstSonarDay) && xt(t)<=datenum(settings.LastSonarDay)
                     label = strvcat([datestr(xt(t),7),'-',datestr(xt(t),3),'-',datestr(xt(t),10)]);
                     text(xt(t),-0.6,label,'horizontalalignment','center','fontsize',12)
                  end
               end
            end
             
            %FanMovie(Fanidx) = getframe(gcf);
            %imfile = sprintf('%s%03d.png',outfile,Fanidx);
          %%%% imfile = sprintf('%s.png',outfile);
            %imfile = sprintf('%s.eps',outfile);
            %print(imfig,'-djpeg90',imfile); % image was a bit muddy, text looked OK
            %set(gcf,'papersize',[8.5 12])
            %set(kids(3),'PlotBoxAspectRatioMode','manual') % don't do this
            %set(kids(2),'PlotBoxAspectRatioMode','manual') % it crunches the time series plots
            %set(kids(1),'PlotBoxAspectRatioMode','manual')
            %set(gcf,'papersize',[8 5.5])
            %set(gcf,'paperposition',[0 0 8 5.5])
            %%set(gcf,'papersize',[4.75 4.75])
            %%set(gcf,'paperposition',[0 0 4.75 4.75])
            %%set(gcf,'paperunits','normalized','renderer','zbuffer')
          %%%%print(imfig,'-dpng',imfile); % best for videomach
            %print(imfig,'-depsc','-r150',imfile); % best for Adobe?
            % set(gcf,'paperunits','inches','renderer','painters')
            % Define the netcdf file using dimensions returned from the
            % first sonar data file
           % if Fanidx == 1
           %    ncdims.time = 0;
           %    ncdims.x = length(Xplot);
           %    ncdims.y = length(Yplot);
           %    nc = defineSonarNcFile(outFileRoot, settings, ncdims);
           % end
            % Write the sonar data to the netcdf file
           % nc = write2sonarNcFile(nc, Fanidx, datenum_sonar, Xplot, Yplot, Zi);
           % Fanidx = Fanidx+1;
           % close(imfig)
            clear imagedata Xplot Yplot Zi
         elseif WhichSonar == 'P' && validFile == 1 && strcmpi(settings.SonartoAnimate,'pencil'),
            % first, display the sonar image
            PencilTime(Pencilidx) = julian([settings.datayear mon day hr mn 0]);
            settings.PencilTime = PencilTime(Pencilidx);
           [PencilData(Pencilidx),settings]=...
               showpencil(fullfile(settings.RootDataDir,SubDataDir(d).name,DataFiles(f).name),...
               settings);
           % save intermediates
            %  save a_pencil_data.mat PencilTime PencilData
            %PencilMovie(Pencilidx) = getframe(gcf);
            figure(2);
            imfig = 2;
            imfile = sprintf('%s.png',outfile);
            %imfile = sprintf('p%03d.png',Pencilidx);
            set(gcf,'papersize',[4.75 4.75])
            set(gcf,'paperposition',[0 0 4.75 4.75])
            set(gcf,'paperunits','normalized','renderer','zbuffer')
            set(gcf,'paperunits','inches','renderer','painters')
            drawnow
            % print(imfig,'-dpng',imfile); % best for videomach
            % pause(0.1)
            %print(imfig,'-depsc','-r150',imfile); % best for Adobe?

            Pencilidx = Pencilidx+1;
            clear outfile imfile
         end
       end
     end
   end

if strcmpi(settings.SonartoAnimate,'fan'),
   % Update/Add NetCDF attributes
   nc.start_time = datestr(gregorian(min(nc{'time'}(1)+nc{'time2'}(1)/3600/1000/24)),0);
   nc.stop_time = datestr(gregorian(max(nc{'time'}(end)+nc{'time2'}(end)/3600/1000/24)),0);
   nc.DELTA_T = [num2str(gmean(diff(sort(nc{'time'}(:)+nc{'time2'}(:)/3600/1000/24)))*24*3600),' sec'];
   hist = nc.history(:);
   hist_new = ['Sonar data displayed as individual frames for an animation',...
      ' and written to NetCDF by animatesonar2nc.m V ',version,...
      '; ',hist];
   nc.history = hist_new;
   nc.NOTE = ['To view images in Matlab type the following at the command ',...
      'prompt:  imagesc(ycoord,xcoord,flipud(squeeze(sonar_image(n,:,:))''));',...
      'where n is the time index'];
   nc = close(nc);
   ncclose;
elseif strcmpi(settings.SonartoAnimate,'pencil'),
   % KLUDGE Alert - should write a netCDF file here
   save pen_data.mat PencilTime PencilData
end

% ---------------- Subfunction: readSonarMeta.m ------------------------- %
function userMeta = readSonarMeta(metaFile);
[atts, defs] = textread(metaFile,'%s %63c','commentstyle','shell');
defs = cellstr(defs);
for i = 1:length(atts)
   theAtt = atts{i}(:)';
   theDef = defs{i}(:)';
   % deblank removes trailing whitespace
   theAtt = deblank(theAtt);
   theDef = deblank(theDef);
   % check for and replace spaces in
   % the attributes with underscores
   f1 = find(isspace(theAtt));
   f2 = strfind(theAtt,'-');
   f = union(f1,f2);
   if ~isempty(f)
      theAtt(f) = '_';
   end
   % attribute definitions read in as characters; convert to
   % numbers where appropriate
   theDefNum = str2double(theDef);
   if ~isnan(theDefNum)
      theDef = theDefNum;
   end
   eval(['userMeta.',theAtt,'= theDef;'])
end

% ---------------- Subfunction:  getSonarYear.m ------------------------- %
function [yyyy] = getSonarYear(DirName, FirstSonarDay, LastSonarDay);

first = datevec(FirstSonarDay);
first_yyyy = first(1);
first_mm = first(2);
first_dd = first(3);

last = datevec(LastSonarDay);
last_yyyy = last(1);
last_mm = last(2);
last_dd = last(3);

mm = str2num(DirName(2:3));
dd = str2num(DirName(4:5));

if (mm >= first_mm)
   yyyy = num2str(first_yyyy);
else
   yyyy = num2str(last_yyyy);
end

% ---------------- Subfunction:  checkSonarFile.m ----------------- %
function [WhichSonar, validFile] = checkSonarFile(sonarFileName);

F = dir(sonarFileName); % FAN files only!
filesize = F.bytes;
if strcmp(F.name(10),'F') && filesize <= 5000
   disp(['### The Fan file ',sonarFileName, ' is a dud file and will not be used in animations'])
   WhichSonar = 'F';
   validFile = 0;
elseif strcmp(F.name(10),'F') && filesize > 5000
   WhichSonar = 'F';
   validFile = 1;
elseif strcmp(F.name(10),'P') && filesize > 100
   WhichSonar = 'P';
   validFile = 1;
else
   validFile = 0;
end

% ---------------- Subfunction:  defineSonarNcFile.m -------------------- %
function nc = defineSonarNcFile(outFileRoot, settings, ncdims);

nc = netcdf([outFileRoot,'.nc'],'clobber');

% dimensions
nc('time') = 0; %unlimited dimension
nc('x') = ncdims.x;
nc('y') = ncdims.y;

% global attributes
% write metadata
metaFields = fieldnames(settings);
for i = 1:length(metaFields)
   theField = metaFields{i};
   theFieldDef = getfield(settings,theField);
   nRows = size(theFieldDef,1);
   nCols = size(theFieldDef,2);
   temp = [];
   if ischar(theFieldDef)
      if nRows>1
         for ii = 1:nRows
            temp = [temp,' ',theFieldDef(ii,1:nCols)];
         end
         theFieldDef = temp;
      end
      eval(['nc.',theField,' = ncchar(theFieldDef);'])
   else
      eval(['nc.',theField,'= ncfloat(theFieldDef);'])
   end
end

% write additional metadata
nc.CREATION_DATE = ncchar(datestr(now));

% variables
nc{'time'} = nclong('time');
nc{'time'}.FORTRAN_format = ncchar('F10.2');
nc{'time'}.units = ncchar('True Julian Day');
nc{'time'}.type = ncchar('UNEVEN');
nc{'time'}.epic_code = nclong(624);
nc{'time'}.FillValue_ = ncfloat(1.00000004091848e+035);

nc{'time2'} = nclong('time') ;
nc{'time2'}.FORTRAN_format = ncchar('F10.2');
nc{'time2'}.units = ncchar('msec since 0:00 GMT');
nc{'time2'}.type = ncchar('UNEVEN');
nc{'time2'}.epic_code = nclong(624);
nc{'time2'}.FillValue_ = ncfloat(1.00000004091848e+035);

nc{'ycoord'} = ncshort('y');
nc{'ycoord'}.long_name = ncchar(['Distance from Sonar, m']);
nc{'ycoord'}.units = ncchar('m');
nc{'ycoord'}.valid_range = ncfloat([0 settings.max_range]);
nc{'ycoord'}.FORTRAN_format = ncchar('F10.2');
nc{'ycoord'}.FillValue_ = ncfloat(settings.fillval);

nc{'xcoord'} = ncshort('x');
nc{'xcoord'}.long_name = ncchar(['Distance from sonar, m']);
nc{'xcoord'}.units = ncchar('m');
nc{'xcoord'}.valid_range = ncfloat([0 settings.max_range]);
nc{'xcoord'}.FORTRAN_format = ncchar('F10.2');
nc{'xcoord'}.FillValue_ = ncfloat(settings.fillval);

% a scale factor is applied to the variable 'sonar_image' in order to store
% the image data as an integer rather than a float.
nc{'sonar_image'} = ncshort('time','y','x');
nc{'sonar_image'}.long_name = ncchar(['Imagenex Sonar Image']);
nc{'sonar_image'}.units = ncchar('');
nc{'sonar_image'}.valid_range = ncfloat([0 log10(256)*1000]);
nc{'sonar_image'}.sensor_type = nc.INST_TYPE(:);
nc{'sonar_image'}.FORTRAN_format = ncchar('F10.2');
nc{'sonar_image'}.FillValue_ = nclong(round(log10(settings.fillval)*1000));
nc{'sonar_image'}.scale_factor = nclong(1000);

% ---------------- Subfunction:  write2sonarNcFile.m -------------------- %
function nc = write2sonarNcFile(nc, fanidx, sonar_datenum, Xplot, Yplot, Zi);

y = nc('y');
ny = size(y,1);
x = nc('x');
nx = size(x,1);

nc{'time'}(fanidx) = floor(julian(datevec(sonar_datenum)));
nc{'time2'}(fanidx) = (julian(datevec(sonar_datenum))-floor(julian(datevec(sonar_datenum)))) * (3600*1000*24);
if fanidx == 1
   nc{'ycoord'}(1:ny) = Yplot;
   nc{'xcoord'}(1:ny) = Xplot;
end
nc{'sonar_image'}(fanidx,1:ny,1:nx) = round(Zi'.*1000); %apply scale factor