function animsonar_nc(metaFile, outRoot)
% animsonar_nc.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:  animsonar_nc('metaFile')
%
%       where:  metaFile is the name of your text file containing metadata,
%                        surrounded by single quotes WITHOUT the file
%                        extension .txt. An example metadata file,
%                        anim_metaexample.txt, is provided in this package of
%                        mfiles.
%               outRoot is the root namae for the output movie files

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

%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 netCDF files of fan and pencil data using
%procsonar_07.m, and to have processed the ADCP data to have created
%waves and currents .nc files.
% output 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
%    ** the matlab M variable with the Movie images is what is saved currently
%       use disp_mvmovie.m to show the movie of the M structure

%close all
more off

version = '1.2'; % Version updated by etm, 9/2007

% 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);
outdir='frame/';

% 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 - assumes image already rotated so N is "up"
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);
   clear ncw 
end
%
% now that the data is loaded, see about making a frame for animation
if fan,
   % load the currents during the desired time
   ncf = netcdf(settings.FanFile);
   timeobj = ncf{'time'};
   timeobj = autonan(timeobj, 1);
   time2obj = ncf{'time2'};
   time2obj = autonan(time2obj, 1);
   tj=timeobj(:)+time2obj(:)./(3600*1000*24);
   datenum_fan=datenum(gregorian(tj));
   % find the indeces to the dates of interest
   findidx = find(datenum_fan < datenum(settings.FirstSonarDay));
   if isempty(findidx), dateidx(1) = 1;
   else dateidx(1) = max(findidx)+1;
   end
   if dateidx(1) < 0, dateidx(1) = 1; end
   findidx = find(datenum_fan > datenum(settings.LastSonarDay));
   if isempty(findidx), dateidx(2) = length(datenum_fan);
   else dateidx(2) = min(findidx)-1;
   end
   % trim time
   % tj = tj(dateidx(1):dateidx(2));
   %datenum_fan = datenum_fan(dateidx(1):dateidx(2));
   % now get the data, replace fill values w/ nans,
   ncvarnames = {'x','y','sonar_image',};
   names = {'fx','fy','fsi'};
   for n = 1:length(ncvarnames)-1;
      ncobj = ncf{ncvarnames{n}};
      ncobj = autonan(ncobj,1);
      eval([names{n},' = ncobj(:,1);'])
   end
   % close(ncf)      % keep open for access to large sonar_image matrix
   clear ncobj tg tj timeobj time2obj
end

   wv_rec=0; fcnt=1;
     figure
   % now loop based on the number of fan images
     p=size(ncf{'sonar_image'});
     % input file starts at 8/27- first good data is 8/29, so skip into the
     % fan data using the indices to start
       for ik=dateidx(1):dateidx(2)
         datenum_sonar = datenum_fan(ik);
         disp(['processing ' datestr(datenum_sonar)])
            %establish the two needed axes
            clf
            % currently set up to generate an avi from the matlab movie M
            set(gcf,'Position',[100 50 800 660])
            sonar_ax=axes('pos',[0.16 0.22 0.75 0.75]);  axis square
               wave_ax=axes('pos',[0.228 0.055 0.62 0.11]);
            %then leave movie2avi(M,'tst_mvco.avi') at the end                       
            
            % this wave_ax is for making the .png output images have equal x
            % lengths.  In the Matlab frames, the wave_axis appears
            % shorter, but it ends up near equal in the .png
            %  wave_ax=axes('pos',[0.255 0.055 0.565 0.11]);
            %      also uncomment the print line on 347

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

            FanDatenum=datenum(datenum_fan(ik));
            % fx and fy are stored as integers, so made an approximation
            xx=linspace(min(fx),max(fx),length(fx)); yy=xx;            
           axes(sonar_ax);
            % imagesc(xx,yy,ncf{'sonar_image'}(ik,:,:))
            % this is the same thing as the imagesc, but won't fail
            % if it can't find an ancestor
            if (size(p) == 3)
              image(xx,yy,ncf{'sonar_image'}(ik,:,:),'CDataMapping','scaled');
            else
              image(xx,yy,ncf{'sonar_image'}(ik,1,:,:),'CDataMapping','scaled');
            end
               hold on
                 colormap gray;    axis square;
                set(sonar_ax,'ydir','normal'); colormap gray; 
                set (sonar_ax,'FontSize',14);
          % add some labels
             yl=ylabel('Fan Beam Sonar, range in meters');
             set(yl,'fontsize',14)
             tt=text(0.8,0.955,'\uparrow North',...
    	       'units','normalized','color','y','fontsize',14)    
             %set(tt,'fontsize',14)
             ts = datestr(datenum_fan(ik));
             tt=text(0.99,0.03,ts,...
	          'units','normalized','color','y','horizontalalignment','right','fontsize',12);
             %set(tt,'fontsize',14)

            if currents,    % add current vector and text
               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(sonar_ax,currX, currY);
                  set(curr_arr,'linewidth',2,'color','r')
                  tt=text(4,4.1,[num2str(round(currSpd)),' cm/s'],'color','r','horizontalalignment','center');
                  set(tt,'FontSize',12);
                  tt=text(4,3.8,[num2str(round(currDir)),' deg'],'color','r','horizontalalignment','center');
                  set(tt,'FontSize',12);
                else
                  text(4,4,['--- cm/s'],'color','r','horizontalalignment','center')
                  text(4,3.5,['--- deg '],'color','r','horizontalalignment','center')
               end
            end  
            
            if waves,   % add current vector and text
               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
                  wv_pts(ik)=waveidx;
                  %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(sonar_ax,waveX, waveY);
                  set(wave_arr,'linewidth',2,'color','c')
                  ww=text(4,3.4,[num2str(round(lambda*500)/5),' cm'],'color','c','horizontalalignment','center');
                    set(ww,'FontSize',12)
                  ww=text(4,3.1,[num2str(round(wvdir(waveidx))),' deg'],'color','c','horizontalalignment','center');
                   set(ww,'FontSize',12)
              else
                  text(4,3,['--- cm  '],'color','c','horizontalalignment','center')
                  text(4,2.5,['--- deg '],'color','c','horizontalalignment','center')
               end
                  axes(wave_ax);  %make the wave axis current
                  plot(datenum_wave,Hsig); hold on;
               if ~isempty(waveidx1) && ~isempty(waveidx2) 
                  wave_dot = plot(datenum_wave(waveidx),Hsig(waveidx),'ro','markersize',4,'MarkerFaceColor','r');
                   wv_rec=wv_rec+1;
               end
                 set (wave_ax,'FontSize',14);
                 yl=ylabel({'Significant','wave','height (m)'});
                 set(yl,'fontsize',14)
                 set(yl,'units','normalized')
                 set(yl,'position',[-.07 .49 0])
                 set(wave_ax,'xlim',[datenum(settings.FirstSonarDay) datenum(settings.LastSonarDay)])
                 % datetick('x',1,'keeplimits')
                 set(wave_ax,'Ylim',[0 2.5])
                 set(wave_ax,'xticklabel',[])
                 xt = get(wave_ax,'xtick');
                 for t = [1 length(xt)]
                     label = strvcat([datestr(xt(t),7),'-',datestr(xt(t),3),'-',datestr(xt(t),10)]);
                     text(xt(t),-0.8,label,'horizontalalignment','center','fontsize',12)
                 end
            end
                        
             %eval(['print -dpng ' outfile]); % best for VideoMach use
             % add the frame to the structure containing the "movie"
             h=gcf;
             M(fcnt)=getframe(h);    % add the gcf to get the entire window
             fcnt=fcnt+1;
       end
       disp([num2str(fcnt) ' frames written'])
      close(ncf) 
    % eval(['save ' outRoot '_movie M']);
    movie2avi(M, 'tmp_run.avi')     % this works for sure
    eval(['movie2avi(M, ''' outRoot '_frames.avi'')'])  % this may not

     %hh=figure; 
     %set(hh,'Position',[100 100 800 800])
     %movie(hh,M)
   ncclose;

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