function procsonar_apr07(metaFile, outFileRoot)
% procsonar_apr07.m  A function to process Imagenex fan and pencil 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, E.Montgomery
% USGS Woods Hole Field Center
% csullivan@usgs.gov
%
% Dependencies:
%   USGS NetCDF Toolbox (C. Denham)
%   showfan07.m (M. Martini, C. Sherwood, E. Montgomery)
%   showpen08.m (M. Martini, C. Sherwood, E.Montgomery)

% etm 04/07
% This program just processes the fan and pencil data to images
% procsonar_apr07 must be run with ('penxxxmeta') to create the
% processed pencilbeam data then run with ('fanxxxmeta').  This creates
% fan_data.mat and pen_data.mat.
%
%  Once the fan and pencil .mat files have been created, animtripod can be
%  run to integrat the data from ADCP, ADV and sonars on one plot, then
%  make a movie from all the samples
%
% 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

close all
more off

version = '2.0'; % Version updated by etm 4/07

% 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');
   % Default settings for Eurostrat first deployment 2002-2003
   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


% 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')
            % create the netcdf file      
            if Fanidx==1
              %  initiaize NetCDF attributes
                    ncdims.x=1001; ncdims.y=1001;
               nc = defineSonarNcFile(outFileRoot, settings, ncdims);
               nc.start_time = datestr(gregorian(FanTime(1)));
               hist = nc.history(:);
               hist_new = ['Sonar data written to NetCDF by procsonar_apr07.m V ',version,...
                '; ',hist];
               nc.history = hist_new;
               nc.NOTE = ['To view images in Matlab type the following at the command ',...
                'prompt:  nc=netdcf(''sonarxxx.nc'');',...
                'imagesc(nc{''ycoord''}(:),nc{''xcoord''}(:),flipud(nc{''sonar_image''}(5,:,:)))',...
                 'where n is the time index'];
                nc = write2sonarNcFile(nc, Fanidx, FanTime(Fanidx), Xplot, Yplot, Zi);
            else
                nc = write2sonarNcFile(nc, Fanidx, FanTime(Fanidx), Xplot, Yplot, Zi);
            end

            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);
            %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
  
            % for pencil, X can be different lengths (if range setting changes, and y and z are
            % dependent on x, but we're going to make all the same length
               xx=[-3.16:.0125:3.16];  % from showpen07
                 xdst=ones(1,length(xx))*nan;
                 yelev=ones(length(xx),1)*nan;
                 zss=ones(length(xx),1)*nan;
                 st=find(PencilData(Pencilidx).intrpx(1) == xx); 
                 xdst(st:st+length(PencilData(Pencilidx).intrpx)-1)=PencilData(Pencilidx).intrpx';
                 yelev(st:st+length(PencilData(Pencilidx).intrpy)-1)=PencilData(Pencilidx).intrpy';
                 zss(st:st+length(PencilData(Pencilidx).intrpz)-1)=PencilData(Pencilidx).intrpz;
             % create the netcdf file
             if Pencilidx==1
              %  initiaize NetCDF attributes
               ncdims.x=length(xx); ncdims.y=1;
               nc = defineSonarNcFile(outFileRoot, settings, ncdims);
               nc.start_time = datestr(gregorian(PencilTime(1)));
               hist = nc.history(:);
               hist_new = ['Sonar data written to NetCDF by procsonar_apr07.m V ',version,...
                '; ',hist];
               nc.history = hist_new;
               % I don't know why, but the flipud is needed here
               nc.NOTE = ['Pencil data has the extracted x and y data, with',...
                        'zi containing the strength of the max reflection defining yelev'];
                nc = write2sonarNcFile(nc, Pencilidx, PencilTime(Pencilidx), xdst, yelev, zss);
            else
                nc = write2sonarNcFile(nc, Pencilidx, PencilTime(Pencilidx), xdst, yelev, zss);
            end
            Pencilidx = Pencilidx+1;
            clear outfile imfile
         end
       end
     end
end

% this is where the data is saved
%  writing to netCDF doesn't work, but you get a .mat file for each fan and
%  Pencil run
if strcmpi(settings.SonartoAnimate,'fan'),
    save fan_data.mat FanData FanTime
       nc.stop_time = datestr(gregorian(FanTime(end)))
       nc.DELTA_T = [num2str(gmean(diff(FanTime*1000))*60),' sec'];
        close(nc);
       ncclose;
elseif strcmpi(settings.SonartoAnimate,'pencil'),
   % KLUDGE Alert - should write a netCDF file here
   save pen_data.mat PencilTime PencilData
       nc.stop_time = datestr(gregorian(PencilTime(end)))
       nc.DELTA_T = [num2str(gmean(diff(PencilTime*1000))*60),' sec'];
        close(nc);
       ncclose;

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);

%fan and pencil have different contents and variable names
if strmatch(settings.SonartoAnimate,'fan') 
   nc{'xcoord'} = ncfloat('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);
   
   nc{'ycoord'} = ncfloat('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);
  % 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'} = ncfloat('time','x','y');
   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_ = settings.fillval;
   nc{'sonar_image'}.scale_factor = nclong(1000);
else
   nc{'xdst'} = ncfloat('time','x');
   nc{'xdst'}.long_name = ncchar(['Distance from sonar, m']);
   nc{'xdst'}.units = ncchar('m');
   nc{'xdst'}.valid_range = ncfloat([0 settings.max_range]);
   nc{'xdst'}.FORTRAN_format = ncchar('F10.2');
   nc{'xdst'}.FillValue_ = ncfloat(settings.fillval);

   nc{'yelev'} = ncfloat('time','x');
   nc{'yelev'}.long_name = ncchar(['Distance from Sonar, m']);
   nc{'yelev'}.units = ncchar('m');
   nc{'yelev'}.valid_range = ncfloat([0 settings.max_range]);
   nc{'yelev'}.FORTRAN_format = ncchar('F10.2');
   nc{'yelev'}.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{'s_strn'} = ncfloat('time','x');
   nc{'s_strn'}.long_name = ncchar(['Imagenex Sonar Image']);
   nc{'s_strn'}.units = ncchar('');
   nc{'s_strn'}.valid_range = ncfloat([0 log10(256)*1000]);
   nc{'s_strn'}.sensor_type = nc.INST_TYPE(:);
   nc{'s_strn'}.FORTRAN_format = ncchar('F10.2');
   nc{'s_strn'}.FillValue_ = settings.fillval;
   nc{'s_strn'}.scale_factor = nclong(1000);
end  

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

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

% timeWrd is in decimal days, so have to convert .604 days to msecs
   nc{'time'}(idx) = floor(timeWrd);
   nc{'time2'}(idx) = mod(timeWrd,floor(timeWrd))*1000*3600*24;
   
if strmatch(nc.SonartoAnimate(:),'fan')
 if idx == 1 
   nc{'ycoord'}(1:ny) = Yplot;
   nc{'xcoord'}(1:nx) = Xplot;
 end
   nc{'sonar_image'}(idx,1:ny,1:nx) = round(Zi.*1000); %apply scale factor
else     % pencil
   nc{'xdst'}(idx,1:nx) = Xplot;
   nc{'yelev'}(idx,1:nx) = Yplot;
   nc{'s_strn'}(idx,1:nx) = Zi;
end
  clear alltime timeWrd
% this command works to plot fan beam imaged:
%  pcolor(nc{'xcoord'}(:),nc{'ycoord'}(:),nc{'sonar_image'}(1,:,:));
%  shading flat
% or
%  imagesc(nc{'xcoord'}(:),nc{'ycoord'}(:),rot90(nc{'sonar_image'}(1,:,:)'));
%  shading flat
% not sure why you HAVE to rot90, but you do.  Or when opening the .nc file
% directly, it looks like flipud is necessary.