function animatesonar2nc(metaFile, outFileRoot)
% animatesonar2nc.m  A function to write Imagenex sonar data to netCDF and
%                    to create a series of images which can be used to 
%                    animate the 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
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Use of this program is described in:
%
% Sullivan, C.M., Warner, J.C., Martini, M.A., Voulgaris, G., 
% Work, P.A., Haas, K.A., and Hanes, D.H. (2006) 
% South Carolina Coastal Erosion Study Data Report for Observations
% October 2003 - April 2004., USGS Open-File Report 2005-1429.
%
% Program written in Matlab v7.1.0 SP3
% Program ran on PC with Windows XP Professional OS.
%
% "Although this program has been used by the USGS, no warranty, 
% expressed or implied, is made by the USGS or the United States 
% Government as to the accuracy and functioning of the program 
% and related program material nor shall the fact of distribution 
% constitute any such warranty, and no responsibility is assumed 
% by the USGS in connection therewith."
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
% Dependencies:
%   USGS NetCDF Toolbox (C. Denham)
%   smart_interp.m (P. Alexander)
%   showfan.m (M. Martini, C. Sherwood)

% C. Sullivan   02/15/06,   version 1.1
% Change handling of metadata information to read this information from an
% m-file. Using an mfile allows metadata fields and information to have
% an unlimited number of characters.
% 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.0';

% Check for metadata file
metaPath = pwd;
meta = dir([metaFile,'.m']);
if isempty(meta)
    fprintf('\n')
    fprintf('The metadata file %s.m 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,'.m']));
    if isempty(meta)
          error('Still cannot find the metadata file ',fullfile(metaPath,[metaFile,'.m']))
    end
end
%metaFile = fullfile(metaPath,meta.name);

% Get user's metadata structure that is contained in
% the variable 'settings'
eval(metaFile)%metaFile;

% Check that the metadata contains required fields.  If a required field
% is missing, ask the user for it.
requiredFields = {'FirstSonarDay','LastSonarDay','WavesFile','RootDataDir',...
             'ADCPbin','CurrentsFile','SectorSweep','datayear',...
             'AngleSweepAround','StepSize','DataPoints','Height',...
             'rot','SonartoAnimate'};
for f = 1:length(requiredFields)
    if ~isfield(settings,requiredFields{f})
        disp(['The field ''',requiredFields{f},''' is not specified in ',metaFile,'.m'])
        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';
    requiredFields = {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   = 'Fields required for the sonar setup';
    lineNo  = 1;
    dlgresult  = inputdlg(prompt,dlgtitle,lineNo,requiredFields);
    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 requiredFields missingFields

% Load currents and wave data
currents = 1;
if exist(settings.CurrentsFile) ~= 2,
    currents = 0;
end

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

if currents, 
    % load the currents during the desired time
    nc = netcdf(settings.CurrentsFile);
    timeobj = nc{'time'};
    time2obj = nc{'time2'};
    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);
        data = ncobj(dateidx(1):dateidx(2),settings.ADCPbin);
        data_int = smart_interp(tj, data, tj, 1);
        eval([names{n},' = data_int;'])
    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_peak','wvdir'};
    names = {'Hsig','Tp','wvdir_from'};
    for n = 1:length(ncvarnames);
        ncobj = ncw{ncvarnames{n}};
        ncobj = autonan(ncobj,1);
        data = ncobj(wavedateidx(1):wavedateidx(2),:);
        data_int = smart_interp(julian_wave, data, julian_wave, 3);
        eval([names{n},' = data_int;'])

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



% 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
            outfile = datestr(datenum_sonar,30);
            outfile = fullfile(outdir,outfile);

            % 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,
                % first, display the sonar image
                FanTime(Fanidx) = julian([settings.datayear mon day hr mn 0]);
                FanDatenum=datenum([settings.datayear mon day hr mn 0]);
                
                imfig = 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 direction from direction FROM to
                        %direction TOWARDS
                        wvdir_to = wvdir_from(waveidx);
                        if wvdir_to > 180
                            wvdir_to = wvdir_from(waveidx) - 180;
                        elseif wvdir_to <= 180
                            wvdir_to = wvdir_from(waveidx) + 180;
                        end
                        %calculate ub (wave orbital velocity) with Sherwood's ubqf.m
                        ub = ubqf( Tp(waveidx), Hsig(waveidx), waterdepth );
                        %amplitude of do (wave orbital diameter) is ub/omega:
                        omega = 2*pi/Tp(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_to);
                        wave_arr = compass(waveX, waveY);
                        set(wave_arr,'linewidth',2,'color','c')
                        if isnan(lambda) | isnan(wvdir_to)
                            text(4,3,['--- cm  '],'color','c','horizontalalignment','center')
                            text(4,2.5,['--- deg '],'color','c','horizontalalignment','center')
                        else
                            text(4,3,[num2str(round(lambda*500)/5),' cm'],'color','c','horizontalalignment','center')
                            text(4,2.5,[num2str(round(wvdir_to)),' deg'],'color','c','horizontalalignment','center')
                        end
                    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])
                    ylim([0 2])
                    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,
                % first, display the sonar image
                PencilTime(Pencilidx) = julian([settings.datayear mon day hr mn 0]);
                if strcmp(settings.SonartoAnimate,'Pencil'),%a double check that yes you want to animate the pencil
                    subplot(2,1,1)
                    showpencil(fullfile(RootDataDir,SubDataDir(d).name,DataFiles(f).name),...
                        'polar',settings);
                    %PencilMovie(Pencilidx) = getframe(gcf);
                    imfile = sprintf('%s%03d.png',outfile,Pencilidx);
                    %print(imfig,'-djpeg90',imfile); % image was a bit muddy, text looked OK
                    print(imfig,'-dpng',imfile); % best for videomach
                end
                Pencilidx = Pencilidx+1;
            end
        end
    end
end

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

% ---------------- 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;
else
    WhichSonar = 'P';
    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


