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,
%                        animexample.txt, is provided in this package of
%                        mfiles.

% Written by Ellyn Montgomery
% USGS Woods Hole Field Center
% emontgomery@usgs.gov
%%This function is heavily based on the script animatesonar.m by M. Martini
%and updatates made by Charlene Sullivan (thanks!). 
%
% Dependencies:
%   USGS NetCDF Toolbox (C. Denham)
%   smartinterp
% This function require users to have run proc_sonar to created fan
% and pencil cdfs, and to have processed the ADCP data to have created
% waves and currents .cdf or .nc 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);'])
    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

if fan
    Fanidx=1;
      nc=netcdf(settings.FanFile)
       tj=nc{'time'}(:)+nc{'time2'}(:)./(3600*1000*24);
       datenum_fan=datenum(gregorian(tj));
     figure
     set(gcf,'Position',[200 200 560 640])
   for jj=1:length(tj)
        fan_ax = axes('pos',[0.2 0.32 0.7 0.611]);
        imagesc(nc{'xcoord'}(:),nc{'ycoord'}(:),squeeze(nc{'sonar_image'}(jj,:,:)))
         set(gca,'Ydir','normal')
         colormap('bone')
        ylabel('range(m)');
        xlabel('range (m)')
        set(gca,'Fontsize',12);
        set(gca,'Xtick',[-5:1:5])
        set(gca,'Ytick',[-5:1:5])
       %axis equal tight
       if(1), % plot range circles
         hold on
        for icirc = 1:5,
	     h = circle(icirc,0,0);
	     set(h,'color',[.7 .7 .6])
        end
       end
         % now set up the labels using time
         tt=nc{'time'}(jj) + nc{'time2'}(jj)/86400000;
         n=gregorian(tt);
         FileName =['S1' num2str(n(2),'%02d') num2str(n(3),'%02d') num2str(n(4),'%02d')];
         FileName = [FileName '.F' num2str(n(5),'%02d') ];
          text(0.01,0.03,FileName,'units','normalized','color','y',...
	      'horizontalalignment','left','fontsize',10)
          text(0.8,0.95,'\uparrow North','units','normalized',...
    	  'color','y','fontsize',10)
          text(0.99,0.03,datestr(n),  'units','normalized','color','y',...
	      'horizontalalignment','right','fontsize',10)
         hold on
         % add the legs of the particle tripod    
         particle_tripod_plot
         % now add currents vector
            if currents,
               curridx1 = find(datenum_curr <= datenum_fan(jj));
               curridx2 = find(datenum_curr >= datenum_fan(jj));
               if ~isempty(curridx1) && ~isempty(curridx2)
                  if datenum_fan(jj)-datenum_curr(curridx1(end))< datenum_curr(curridx2(1))-datenum_fan(jj)
                     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
         % now add waves vector
            if waves,
               xx=find(datenum_wave <= datenum_fan(jj));
               df=datenum_wave(xx(end):xx(end)+1)-datenum_fan(jj);
               idx=find(min(abs(df))==df);
               if idx==1 
                   waveidx=xx(end);
               else
                   waveidx=xx(end)+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
                  wave_dot = plot(datenum_wave(waveidx),Hsig(waveidx),'r.','markersize',14);
               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
    end
            
         if WhichSonar == 'P' && validFile == 1 && strcmpi(settings.SonartoAnimate,'pencil'),
            % first, display the sonar image
                % for after the netcdf of the Pencil is read
                %   tj=nc{'time'}(:)+nc{'time2'}(:)./(3600*1000*24);
                %   datenum_pen=datenum(gregorian(tj));

            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

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