function animtripod(metaFile)
% animtripod.m  A function to make a movie of all the frames of data from
% ADCP, and and pencil sonar and other sensors on a tripod
%
% 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.

% Written by Ellyn Montgomery
% USGS Woods Hole Field Center
% emontgomery@usgs.gov
%
% Dependencies:
%   USGS NetCDF Toolbox (C. Denham)

%This function is heavily based on the script animatesonar.m by M. Martini
%and updatates made by Charlene Sullivan (thanks!). Use of this function 
%require users to have created fan_data.mat and pencildata.mat using
%procsonar_apr07.m, and to have processed the ADCP data to have created
%waves and currents files.
%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


close all
more off

version = '1.2'; % Version updated by ETM June 2007
% reads fan data from .cdf file creaded by procsonar

% Check for metadata file
% this one should contain 4 things: file names for currents, waves, pencil
% and fan data.  If there is no data for any of the four, use a fake name,
% thae doesn't exist.
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.
currents=1; waves=1; pencil=1; fan=1;

% See what data we're combining
if exist(settings.CurrentsFile) ~= 2,
   currents = 0;
end
if exist(settings.WavesFile) ~= 2,
   waves = 0;
end   
if exist(settings.PenFile) ~= 2,
   pencil = 0;
end
if exist(settings.FanFile) ~= 2,
  fan = 0;
end


%get the adcp data (could be adv)
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

% get the waves data
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, FanData(Fanidx), settings] = ...
               showfan07(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)
            Fanidx=Fanidx+1;
            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]=...
               showpen07(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

% ---------------- 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