function theMeta = aqamat2cdfbychan(theMeta)
% aqamat2cdf - A function to convert Aquascat ABS data (*.aqa.mat) to netCDF.
%
% Only for Aquascat 1000
%
% function theMeta = aqamat2cdf(theMeta)
% prerequisite:
%       Aquatec's MATLAB mfiles for the ABSS
%       aqa2mat.m to use Aquatec's mfiles to make .mat files and collect
%                   metadata.  This is the easier way to contruct the
%                   theMeta structure and create the .mat files.  Or, one
%                   can construct theMeta (below) by hand and run
%                   ExtractAllAquascat1000Files with no areguments
%
% usage:  aqa2cdfbychan(theMETA);
%
% where:
% if theMeta is a .mat file on disk, theMeta structure is loaded from it
% if theMeta is a structure, it must be as follows
% theMeta =
%                MOORING: 847
%        Deployment_date: '04-dec-2007 22:03:40 GMT'
%          Recovery_date: '10-Apr-2008 18:00:00 GMT'
%              longitude: -118.3088
%               latitude: 33.6838
%     magnetic_variation: 12.9000
%            WATER_DEPTH: 57
%       WATER_DEPTH_NOTE: 'from ships navigation at deployment '
%            DATA_ORIGIN: 'USGS WHSC and PSC Sed Trans Groups'
%             EXPERIMENT: '2007-2008 PV BBL Experiment'
%            DESCRIPTION: [1x101 char]
%                PROJECT: 'USGS/EPA Palos Verdes'
%              Data_Cmnt: 'geoprobe tripod at site B6'
%                  SciPi: 'C. Sherwood, M. Noble, J. Xu, K. Rosenberger'
%              ABSSmodel: 'AQUAscat1000R                   '
%        input_directory: 'E:\PV_2007\847_B6_geoprobe\8476abss\testdata'
%       output_directory: 'E:\PV_2007\847_B6_geoprobe\8476abss\'
%            outFileRoot: '8476abss'
%           metaFileName: '8476abssMeta.mat'
%          sensor_height: [1.0350 1.0350 1.0350]
%       ABSSSerialNumber: '706-30'
%       ABSSBoardVersion: '3'
%      AbsTransducerName: [3x3 char]
%              frequency: [980000 2500000 5000000]
%       AbsTxPulseLength: [1.3333e-005 1.3333e-005 1.3333e-005]
%          BurstDuration: [1 1 1]
%               PingRate: 64
%               NumPings: 53760
%             AbsAverage: [64 64 64]
%         AbsBinLengthMM: [5 5 5]
%            AbsBinRange: [256x3 double]
%             AbsNumBins: [256 256 256]
%         AbsNumProfiles: [840 840 840]
%         AbsProfileRate: [1 1 1]
%        AbsStartingGain: [0 0 0]
%                  nBins: 256
%              nProfiles: 840
%               nXducers: 3
%            Conventions: 'PMEL/EPIC'
%           DATA_SUBTYPE: 'MOORED'
%           COORD_SYSTEM: 'GEOGRAPHIC'
%             WATER_MASS: '?'
%              POS_CONST: 0
%            DEPTH_CONST: 0
%                DRIFTER: 0
%               VAR_FILL: 1.0000e+035
%              FILL_FLAG: 0
%              COMPOSITE: 0
%              INST_TYPE: 'Acoustic Backscatter Sensor'
%                pctiles: [5x1 double]
%               nPctiles: 5
%
% It takes about 1 hour to convert 2200 files.
%

%
% Written by Marinna Martini based on original code by Charlene Sullivan
% USGS Woods Hole Field Center
% csullivan@usgs.gov
%

%
% 5/8/08 MM change over to struct driven script input method, use Aquatec's
% ABSS toolbox to read data so that format changes are no longer USGS'
% problem
% C. Sullivan   08/22/05,   version 1.5
% Add capeability to process USC ABS data.  Read this data with
% read_aquatec_v4.m.  This code is a modified form of the code
% read_aquatec_V3_2.m provided by G. Voulgaris (USC).  Include metadata
% output to 'ABS_info' by readaqa_oleg.m and read_aquatec_v4.m in metadata.
% Change netcdf variables 'Mdata*' to 'abs_trans*'.
% C. Sullivan   07/11/05,   version 1.4
% Add the attributes 'input_directory', 'output_directory', and 'Conventions'
% to the metadata file. Users should place their metadata file in the
% directory specified by the attribute 'output_directory', and run this
% mfile from that directory. Remove the attributes 'SampleRate', 'SerialNumber',
% 'BuildDate', 'DSPVersion', and 'FPGAVersion' from the metadata file because
% these attributes will be read from the .aq(a)(f) files. Removing median
% variable form statistics netCDF file b/c these values are given by the 50th
% percentile.
% C. Sullivan   07/07/05,   version 1.3
% Some netCDF files are not being closed appropriately. Getting rid of
% 'isunix' loop for finding files as we can use the 'dir' command on both
% unix and windows with the same result.
% C. Sullivan   07/05/05,   version 1.2
% This version assumes the user has deleted extra .aqa or .aqf files that
% were collected prior to (after) instrument deployment (recovery). We cannot
% get burst number from 'ABS_info.burstnum', output by readaqa_oleg.m, because
% these are not sequentially numbered from 1:nBursts.  For example, the 300th
% .aqa file has an ABS_info.burstnum=44 and ABS_info.burstnum appears to
% repeats after every 255th burst.
% C. Sullivan   07/01/05,   version 1.1
% Calculate percentiles for the statistics netCDF file. Calculate statistics
% using both the raw amplitudes and the raw amplitudes squared. This version
% will also load ABS data downloaded from the instrument as a series of
% *.aqa and/or *.aqf files.
% C. Sullivan   06/28/05,   version 1.0
% This function converts ABS data, downloaded from the instrument as a series
% of *.aqa files, to raw data burst and statistic netCDF (*.cdf) files. It
% assumes all metadata is defined in a simple text file (see absmetaexample.txt),
% and the user is reading a USGS ABS with the mfile readaqa_oleg.m. A maximum
% size limit of 500 MB is imposed on the burst netCDF files, therefore there
% may be multiple burst netCDF files created. It has been possible to write
% all the statistics to a single statistics netCDF file.

% TODO - ABSS returns nprofiles & nbins for each transducer, make use of
% this so that when it's different for each, it's OK
% TODO:  Correct the units attribute
% TODO:  Correct the valid range attribute
% TODO: consider these calculations
%   9.3	RangeCorrectAbsProfiles

more off

% get the current SVN version- the value is automatically obtained in svn
% is the file's svn.keywords which is set to "Revision"
rev_info = 'SVN $Revision$';
fprintf('%s %s running\n',mfilename,rev_info)

% There is a 500 MB size limit imposed on the raw data burst netCDF (.cdf)
% files. Additional raw data burst netCDF files are created as needed.
maxBytes = 750000000;

% mandatory inputs
if ~exist('theMeta','var'),
    theMeta = [];
end

% ## directory with .aqa and .aqf files (required)
if ~isfield(theMeta,'input_directory') || isempty(theMeta.input_directory),
    theMeta.input_directory = uigetdir(pwd,'Where are the .aqa files?');
    if ~theMeta.input_directory, theMeta.input_directory= []; return; end
end
if ~isfield(theMeta,'output_directory') || isempty(theMeta.output_directory),
    theMeta.output_directory = uigetdir(pwd, 'Where do you want the .cdf file to go?');
end

if ~isfield(theMeta,'metaFileName') ||  ~exist(theMeta.metaFileName,'file'),
    disp('Metadata needed, did you run aqa2mat first?')
    return
else
    if isfield(theMeta,'input_directory'), userdir = theMeta.input_directory; end
    fprintf('Loading metadata from %s\n',theMeta.metaFileName)
    load(theMeta.metaFileName);
    if ~strcmp(userdir,theMeta.input_directory),
        str = {userdir, theMeta.input_directory};
        [s,ok] = listdlg('PromptString','Select an input directory:',...
            'SelectionMode','single',...
            'ListString',str);
        if ok, theMeta.input_directory = str{s}; end
    end
end

% global attributes which change
gattnames = {'MOORING', 'Deployment_date', 'Recovery_date', 'longitude', 'latitude',...
    'magnetic_variation', 'WATER_DEPTH', 'DATA_ORIGIN', 'EXPERIMENT','DESCRIPTION',...
    'PROJECT', 'WATER_DEPTH_NOTE','Data_Cmnt','SciPi','ABSSmodel'};
for ifield = 1:length(gattnames),
    if ~isfield(theMeta,gattnames{ifield}),
        fprintf('Missing metadata for %s\n',gattnames{ifield})
        eval(sprintf('theMeta.%s = [];',gattnames{ifield}))
    end
end
% fixed globals
theMeta.Conventions = 'PMEL/EPIC'; % ## variable conventions, ie: EPIC'
theMeta.DATA_SUBTYPE = 'MOORED'; % ## description of DATA_TYPE
theMeta.COORD_SYSTEM = 'GEOGRAPHIC'; % ## coordinate system of the data   '
theMeta.WATER_MASS = '?'; % ## water mass flag used for EPIC contouring programs
theMeta.POS_CONST = 0; % ## consistent position flag (=1 not consistent)
theMeta.DEPTH_CONST = 0; % ## consistent depth flag (=1 not consistent)
theMeta.DRIFTER = 0; % ## drifter flag (=1 if drifter)
theMeta.VAR_FILL = 1.E35; % ## missing or bad data value identifier
theMeta.FILL_FLAG = 0; % ## data fill flag (=1 if data has fill values)
theMeta.COMPOSITE = 0; % ## number of pieces in composite series
theMeta.INST_TYPE = 'Acoustic Backscatter Sensor';
if ischar(theMeta.longitude), theMeta.longitude = str2double(theMeta.longitude); end
if ischar(theMeta.latitude), theMeta.latitude = str2double(theMeta.latitude); end

% Check for Aquascat MATLAB files
inputDir = theMeta.input_directory;
aqamat = dir(fullfile(inputDir,'*.aqa.mat'));
nMatFiles = size(aqamat,1);
if isempty(aqamat),
    fprintf('\n')
    error(['Aquascat ABS MATLAB data files do not exist in ',inputDir]);
else
    matFiles = cell(nMatFiles);
    for n = 1:nMatFiles
        matFiles{n} = aqamat(n).name;
    end
end

% The 1st, 15th, 50th, 84th, and 99th percentiles will be calculated for
% the statistics file. Include this information in the user's metadata.
pctiles = [1; 16; 50; 84; 99];
theMeta.pctiles = pctiles;
theMeta.nPctiles = length(theMeta.pctiles);

burstFileNum = ones(theMeta.nXducers,1); %initialize burst netcdf file #
lastBrec = 0; %initialize burst record #
AbsData = []; % this is to quiet mlint
AbsMean = [];
cdfb = cell(theMeta.nXducers,1);
cdfs = cell(theMeta.nXducers,1);
burstFile = cell(theMeta.nXducers,1);
statFile = cell(theMeta.nXducers,1);

disp(theMeta)

tic

% The following loop reads the data in each Aquascat .mat data file,
% writes the % burst data to the burst netCDF file,
% performs basic statistics on the burst, and writes
% the statistics to the statistics netCDF file.
for n = 1:nMatFiles,
    
    aqFile = fullfile(inputDir,matFiles{n}); %burst file name
    brec = n - lastBrec; %burst record number
    srec = n; %statistic record number
    
    %Read the data
    load(aqFile)
    
    if n==1,
        fprintf('\n')
        fprintf('There are %d bursts to write to netCDF\n',nMatFiles)
        % define the output files
        %Create burst and statistic netCDF files
        for ixd = 1:theMeta.nXducers,
            [cdfb{ixd}, burstFile{ixd}] = defineAbsBurstCdf(theMeta, burstFileNum(ixd), ixd);
            fprintf('Writing channel %d burst data to the file %s\n',ixd,burstFile{ixd})
            [cdfs{ixd}, statFile{ixd}] = defineAbsStatCdf(theMeta, ixd);
            fprintf('Writing channel %d statistics data to the file %s\n',ixd,statFile{ixd})
        end
    end
    
    fprintf('reading burst %d, %s elapsed %f min\n',BurstNumber, BurstTime, toc/60)
    
    %Write Aquascat data to the burst netCDF files
    % BurstTime = '14-Jan-2008 06:09:54'
    profStartTime = julian(datevec(datenum(BurstTime))); % profile start time, julian days
    sampleRate = theMeta.PingRate(1);        %sample rate
    samples_per_avg = theMeta.AbsAverage(1);    %samples per average
    profDt = (1/(sampleRate/samples_per_avg))/3600/24;
    profTime = profStartTime : profDt : profStartTime+(profDt*(theMeta.nProfiles-1));
    for ixd = 1:theMeta.nXducers,
        cdfb{ixd}{'time'}(brec, :) = floor(profTime);
        cdfb{ixd}{'time2'}(brec, :) = (profTime-floor(profTime)).*(24*3600*1000);
        if ~isequal(BurstNumber,n) %if burst #'s begin to repeat,
            cdfb{ixd}{'burst'}(brec) = n;     %override them
        else
            cdfb{ixd}{'burst'}(brec) = BurstNumber;
        end
        % shape of AbsData is [nbins x nprofiles x ntransducers]
        %cdfb{'abs_trans1'}(brec, 1:theMeta.nProfiles, 1:nBins) = AbsData(:,:,1); %Amplitude at frequency 1
        %cdfb{'abs_trans2'}(brec, 1:theMeta.nProfiles, 1:nBins) = AbsData(:,:,2); %Ampliude at frequency 2
        %cdfb{'abs_trans3'}(brec, 1:theMeta.nProfiles, 1:nBins) = AbsData(:,:,3); %Ampliude at frequency 3
        cdfb{ixd}{'abs_trans'}(brec, :, :) = AbsData(:,:,ixd)'; %Amplitude at frequency 1
    end
    
    %Get the size of the burst netCDF file after the
    %1st Aquascat file is written. Use this size to determine
    %the number of Aquascat files you can write to each
    %burst netCDF file
    if n == 1,
        nAqFilesPerBurstCdf = zeros(theMeta.nXducers,1);
        for ixd = 1:theMeta.nXducers,
            filechk = dir(burstFile{ixd});
            bytesPerAqFile = filechk.bytes;
            nAqFilesPerBurstCdf(ixd) = floor(maxBytes/bytesPerAqFile);
        end
    end
    
    %write statistics to the statistic netCDF file
    for ixd = 1:theMeta.nXducers,
        %cdfs{ixd}{'time'}(srec) = profTime(1);%use first time stamp in each burst
        cdfs{ixd}{'time'}(srec) = floor(profStartTime);%use first time stamp in each burst
        %cdfs{ixd}{'time2'}(srec) = (profTime(1)-floor(profTime(1))).*(24*3600*1000);
        cdfs{ixd}{'time2'}(srec) = (profStartTime-floor(profStartTime)).*(24*3600*1000);
        if ~isequal(BurstNumber,n) %if burst #'s begin to repeat,
            cdfs{ixd}{'burst'}(srec) = n;     %override them
        else
            cdfs{ixd}{'burst'}(srec) = BurstNumber;
        end
        cdfs{ixd}{'abs_trans_mean'}(srec, :) = AbsMean(:,ixd);
        cdfs{ixd}{'abs_trans_mean_sq'}(srec, :) = gmean(AbsData(:,:,ixd)'.^2);
        cdfs{ixd}{'abs_trans_std'}(srec, :) = gstd(AbsData(:,:,ixd)');
        cdfs{ixd}{'abs_trans_std_sq'}(srec, :) = std(AbsData(:,:,ixd)'.^2);
        cdfs{ixd}{'abs_trans_pctl'}(srec, :, :) = ...
            matpctile(AbsData(:,:,ixd)', pctiles(:))';
        cdfs{ixd}{'abs_trans_pctl_sq'}(srec, :, :) = ...
            matpctile(AbsData(:,:,ixd)'.^2, pctiles)';
        % Ken Foote's stats
        xi = sum(AbsData(:,:,ixd),1); % sum along profiles
        s1 = sum(xi);
        s2 = sum(xi.^2);
        [r, n] = size(AbsData(:,:,ixd));
        sav = s1/n;
        ssd = sqrt((s2-n*sav^2)/(n-1));
        
        cdfs{ixd}{'sav'}(srec) = sav;
        cdfs{ixd}{'ssd'}(srec) = ssd;
    end
    
    %Update the user with the code's progress
    if mod(n,50)==0
        fprintf('Finished writing burst %d; %6.2f minutes elapsed\n',...
            n,toc/60)
    end
    
    for ixd = 1:theMeta.nXducers,
        %If n equals the number of bursts you can write to the burst netCDF
        %file, close this burst netCDF file and create another one to which
        %you will write subsequent bursts
        if mod(n,nAqFilesPerBurstCdf(ixd))==0
            lastBrec = cdfb{ixd}{'burst'}(end);
            %Update the attributes of the current burst netCDF file
            fprintf('The burst file %s is full\n',burstFile{ixd})
            cdfb{ixd}.burstnum = [];
            start_time = cdfb{ixd}{'time'}(1,1) + cdfb{ixd}{'time2'}(1,1)/3600/1000/24;
            cdfb{ixd}.start_time = ncchar(datestr(gregorian(start_time)));
            stop_time = cdfb{ixd}{'time'}(end,end) + cdfb{ixd}{'time2'}(end,end)/3600/1000/24;
            cdfb{ixd}.stop_time = ncchar(datestr(gregorian(stop_time)));
            hist = cdfb{ixd}.history(:);
            hist_new = ['Converted to netCDF via Matlab by %s V ', mfilename, rev_info,'; ',hist];
            cdfb{ixd}.history = ncchar(hist_new);
            close(cdfb{ixd});
            %Create a new burst netCDF file for subsequent bursts
            burstFileNum(ixd) = burstFileNum(ixd) + 1;
            [cdfb{ixd}, burstFile{ixd}] = defineAbsBurstCdf(theMeta, burstFileNum(ixd), ixd);
            fprintf('Now writing burst data to the file %s\n',burstFile{ixd})
        end
    end
end

for ixd = 1:theMeta.nXducers,
    %if mod(n,nMatFiles)==0
    %if mod(n,100)==0
    %Update the attributes of the last burst netCDF file
    cdfb{ixd}.burstnum = [];
    start_time = cdfb{ixd}{'time'}(1,1) + cdfb{ixd}{'time2'}(1,1)/3600/1000/24;
    if ~isempty(start_time),
        cdfb{ixd}.start_time = ncchar(datestr(gregorian(start_time)));
        stop_time = cdfb{ixd}{'time'}(end,end) + cdfb{ixd}{'time2'}(end,end)/3600/1000/24;
        cdfb{ixd}.stop_time = ncchar(datestr(gregorian(stop_time)));
    end
    
    hist = cdfb{ixd}.history(:);
    hist_new = ['Converted to netCDF via Matlab by aqa2cdf.m V ',rev_info,'; ',hist];
    cdfb{ixd}.history = ncchar(hist_new);
    close(cdfb{ixd});
end

%update the attributes of the statistics netCDF file and
%calculate min/max values
for ixd = 1:theMeta.nXducers,
    cdfs{ixd}.burstnum = [];
    start = cdfs{ixd}{'time'}(1) + cdfs{ixd}{'time2'}(1)/3600/1000/24;
    cdfs{ixd}.start_time = ncchar(datestr(gregorian(start)));
    stop = cdfs{ixd}{'time'}(end) + cdfs{ixd}{'time2'}(end)/3600/1000/24;
    cdfs{ixd}.stop_time = ncchar(datestr(gregorian(stop)));
    delta_t = mean(diff(cdfs{ixd}{'time'}(:) + cdfs{ixd}{'time2'}(:)/3600/1000/24));
    cdfs{ixd}.DELTA_T = ncchar([num2str(delta_t*24*3600),' s']);
    hist = cdfs{ixd}.history(:);
    hist_new = ['Converted to netCDF via Matlab by aqa2cdf.m V ',rev_info,'; ',hist];
    cdfs{ixd}.history = ncchar(hist_new);
    calc_minmax_vals(cdfs{ixd});
    close(cdfs{ixd});
end

fprintf('Finished writing all bursts. %6.2f minutes elapsed\n',toc/60)
%end

%% ---------------------- Subfunction ------------------------------------ %
function [cdfb, burstFile] = defineAbsBurstCdf(theMeta, burstFileNum, ixdr)

outputDir = theMeta.output_directory;
buf = sprintf('%s%d_%db.cdf',theMeta.outFileRoot,ixdr,burstFileNum);
burstFile = fullfile(outputDir,buf);

cdfb = netcdf(burstFile,'clobber');

% spearate some stuff we know we don't want in the global attributes
omitGatts = {'pctiles','nPctiles','input_directory', 'output_directory',...
    'sensor_height','outFileRoot','metaFileName'};
for i = 1:length(omitGatts),
    if isfield(theMeta, omitGatts{i}),
        eval(sprintf('miscMeta.%s = theMeta.%s;',omitGatts{i},omitGatts{i}))
        theMeta = rmfield(theMeta, omitGatts{i});
    end
end
% this cause multiple files - too many - burst files to be written
%theMeta = check_gmeta(theMeta);

% write metadata
metaFields = fieldnames(theMeta);
for i = 1:length(metaFields)
    %theField = metaFields{i};
    %theFieldDef = getfield(theMeta,theField);
    theFieldDef = theMeta.(metaFields{i});
    nRows = size(theFieldDef,1);
    nCols = size(theFieldDef,2);
    temp = [];
    if ischar(theFieldDef)
        if nRows == theMeta.nXducers,
            theFieldDef = theFieldDef(ixdr,1:nCols);
        elseif nRows > theMeta.nXducers
            for ii = 1:nRows
                temp = [temp,' ',theFieldDef(ii,1:nCols)];
            end
            theFieldDef = temp;
        end
        theCommand = ['cdfb.',metaFields{i},' = ncchar(theFieldDef);'];
    else
        if length(theFieldDef) == theMeta.nXducers,
            theCommand = ['cdfb.',metaFields{i},'= ncfloat(theFieldDef(ixdr));'];
        else
            theCommand = ['cdfb.',metaFields{i},'= ncfloat(theFieldDef);'];
        end
    end
    %disp(theCommand)
    eval(theCommand)
end

% write additional metadata
cdfb.CREATION_DATE = ncchar(datestr(now));
cdfb.DESCRIPT = ncchar('Aquascat ABS raw data burst file');
cdfb.DATA_TYPE = ncchar('ABS');

cdfb('time') = 0; %unlimited
cdfb('depth') = theMeta.nBins;
cdfb('lat') = 1;
cdfb('lon') = 1;
cdfb('profile') = theMeta.nProfiles;

% define burst variables
cdfb{'time'} = nclong('time','profile');
cdfb{'time'}.FORTRAN_format = ncchar('F10.2');
cdfb{'time'}.units = ncchar('True Julian Day');
cdfb{'time'}.type = ncchar('EVEN');
cdfb{'time'}.epic_code = nclong(624);

cdfb{'time2'} = nclong('time','profile') ;
cdfb{'time2'}.FORTRAN_format = ncchar('F10.2');
cdfb{'time2'}.units = ncchar('msec since 0:00 GMT');
cdfb{'time2'}.type = ncchar('EVEN');
cdfb{'time2'}.epic_code = nclong(624);

% we need depth.
cdfb{'depth'} = ncfloat('depth');  ncobj = cdfb{'depth'};
ncobj.FORTRAN_format = ncchar('F10.2');
ncobj.units = ncchar('m');
ncobj.type = ncchar('EVEN');
ncobj.name = ncchar('D');
ncobj.long_name = ncchar('DEPTH (m)');
ncobj.generic_name = ncchar('depth');
ncobj.epic_code = nclong(3);
ncobj.bin_size = ncdouble(theMeta.AbsBinLengthMM/1000);
ncobj.xducer_offset_from_bottom = ncdouble(miscMeta.sensor_height(1));
ncobj(:) = theMeta.AbsBinRange(:,ixdr)+theMeta.WATER_DEPTH(1);

cdfb{'lon'} = ncfloat('lon'); ncobj = cdfb{'lon'};
ncobj.FORTRAN_format = ncchar('f10.4');
ncobj.units = ncchar('degree_east');
ncobj.type = ncchar('EVEN');
ncobj.epic_code = nclong(502);
ncobj.name = ncchar('LON');
ncobj.long_name = ncchar('LONGITUDE');
ncobj.generic_name = ncchar('lon');
ncobj(:) = theMeta.longitude;

cdfb{'lat'} = ncfloat('lat');  ncobj = cdfb{'lat'};
ncobj.FORTRAN_format = ncchar('F10.2');
ncobj.units = ncchar('degree_north');
ncobj.type = ncchar('EVEN');
ncobj.epic_code = nclong(500);
ncobj.name = ncchar('LAT');
ncobj.long_name = ncchar('LATITUDE');
ncobj.generic_name = ncchar('lat');
ncobj(:) = theMeta.latitude;

cdfb{'burst'} = nclong('time');
cdfb{'burst'}.FORTRAN_format = ncchar('F10.2');
cdfb{'burst'}.units = ncchar('counts');
cdfb{'burst'}.type = ncchar('EVEN');
cdfb{'burst'}.FillValue_ = ncfloat(1.e+035);

varname = 'vrange';
cdfb{varname} = ncfloat('depth');
buf = sprintf('Range from Transducer at %4.3f Mhz',theMeta.frequency(ixdr)./1e6);
cdfb{varname}.long_name = ncchar(buf);
buf = sprintf('%2.1f Mhz Range',theMeta.frequency(ixdr)./1e6);
cdfb{varname}.name = ncchar(buf);
cdfb{varname}.generic_name = ncchar(' ');
cdfb{varname}.epic_code = 0;
cdfb{varname}.units = ncchar('m');
cdfb{varname}.valid_range = ncfloat([0 3]);
cdfb{varname}.sensor_type = ncchar(theMeta.INST_TYPE);
cdfb{varname}.initial_sensor_height = miscMeta.sensor_height(ixdr);
cdfb{varname}.sensor_depth = theMeta.WATER_DEPTH - miscMeta.sensor_height(ixdr);
%cdfb{varname}.serial = ncchar(cdfb.serial_number(:));
cdfb{varname}.frequency = ncchar(sprintf('%2.1f Mhz',theMeta.frequency(ixdr)./1e6));
cdfb{varname}.FORTRAN_format = ncchar('F10.2');
cdfb{varname}.FillValue_ = ncfloat(1.e+035);
cdfb{varname}(:) = theMeta.AbsBinRange(:,ixdr);


varname = 'abs_trans';
cdfb{varname} = ncfloat('time','profile','depth');
buf = sprintf('Amplitude at %4.3f Mhz',theMeta.frequency(ixdr)./1e6);
cdfb{varname}.long_name = ncchar(buf);
buf = sprintf('%2.1f Mhz ',theMeta.frequency(ixdr)./1e6);
cdfb{varname}.name = ncchar(buf);
cdfb{varname}.generic_name = ncchar(' ');
cdfb{varname}.epic_code = 0;
cdfb{varname}.units = ncchar('ABS units');
cdfb{varname}.valid_range = ncfloat([0 1e16]);
cdfb{varname}.sensor_type = ncchar(theMeta.INST_TYPE);
cdfb{varname}.initial_sensor_height = miscMeta.sensor_height(ixdr);
cdfb{varname}.sensor_depth = theMeta.WATER_DEPTH - miscMeta.sensor_height(ixdr);
%cdfb{varname}.serial = ncchar(cdfb.serial_number(:));
cdfb{varname}.frequency = ncchar(sprintf('%2.1f Mhz',theMeta.frequency(ixdr)./1e6));
cdfb{varname}.FORTRAN_format = ncchar('F10.2');
cdfb{varname}.FillValue_ = ncfloat(1.e+035);

endef(cdfb)

return

%% ---------------------- Subfunction ------------------------------------
function [cdfs, statFile] = defineAbsStatCdf(theMeta, ixdr)

outputDir = theMeta.output_directory;
%statFile = fullfile(outputDir,[theMeta.outFileRoot,'s.cdf']);
buf = sprintf('%s%ds.cdf',theMeta.outFileRoot,ixdr);
statFile = fullfile(outputDir,buf);

cdfs = netcdf(statFile,'clobber');

% spearate some stuff we know we don't want in the global attributes
omitGatts = {'pctiles','nPctiles','input_directory',...
    'output_directory','sensor_height','outFileRoot','metaFileName'};
for i = 1:length(omitGatts),
    if isfield(theMeta, omitGatts{i}),
        eval(sprintf('miscMeta.%s = theMeta.%s;',omitGatts{i},omitGatts{i}))
        theMeta = rmfield(theMeta, omitGatts{i});
    end
end
% write metadata
metaFields = fieldnames(theMeta);
for i = 1:length(metaFields)
    %theField = metaFields{i};
    %theFieldDef = getfield(theMeta,theField);
    theFieldDef = theMeta.(metaFields{i});
    nRows = size(theFieldDef,1);
    nCols = size(theFieldDef,2);
    temp = [];
    if ischar(theFieldDef)
        if nRows == theMeta.nXducers,
            theFieldDef = theFieldDef(ixdr,1:nCols);
        elseif nRows > theMeta.nXducers
            for ii = 1:nRows
                temp = [temp,' ',theFieldDef(ii,1:nCols)];
            end
            theFieldDef = temp;
        end
        theCommand = ['cdfs.',metaFields{i},' = ncchar(theFieldDef);'];
    else
        if length(theFieldDef) == theMeta.nXducers,
            theCommand = ['cdfs.',metaFields{i},'= ncfloat(theFieldDef(ixdr));'];
        else
            theCommand = ['cdfs.',metaFields{i},'= ncfloat(theFieldDef);'];
        end
    end
    eval(theCommand)
end


% write additional metadata
cdfs.CREATION_DATE = ncchar(datestr(now));
cdfs.DESCRIPT = ncchar('Aquascat ABS raw data statistics file');
cdfs.DATA_TYPE = ncchar('ABS');

cdfs('time') = 0; %unlimited
cdfs('depth') = theMeta.nBins;
cdfs('lat') = 1;
cdfs('lon') = 1;
cdfs('pctile') = miscMeta.nPctiles;

% define stats variables
cdfs{'time'} = nclong('time');
cdfs{'time'}.FORTRAN_format = ncchar('F10.2');
cdfs{'time'}.units = ncchar('True Julian Day');
cdfs{'time'}.type = ncchar('UNEVEN');
cdfs{'time'}.epic_code = nclong(624);

cdfs{'time2'} = nclong('time') ;
cdfs{'time2'}.FORTRAN_format = ncchar('F10.2');
cdfs{'time2'}.units = ncchar('msec since 0:00 GMT');
cdfs{'time2'}.type = ncchar('UNEVEN');
cdfs{'time2'}.epic_code = nclong(624);

% we need depth.
cdfs{'depth'} = ncfloat('depth');  ncobj = cdfs{'depth'};
ncobj.FORTRAN_format = ncchar('F10.2');
ncobj.units = ncchar('m');
ncobj.type = ncchar('EVEN');
ncobj.name = ncchar('D');
ncobj.long_name = ncchar('DEPTH (m)');
ncobj.generic_name = ncchar('depth');
ncobj.epic_code = nclong(3);
ncobj.bin_size = ncdouble(theMeta.AbsBinLengthMM/1000);
ncobj.xducer_offset_from_bottom = ncdouble(miscMeta.sensor_height(1));
ncobj(:) = theMeta.WATER_DEPTH(1)- miscMeta.sensor_height(1)+ theMeta.AbsBinRange(:,ixdr);

cdfs{'lon'} = ncfloat('lon'); ncobj = cdfs{'lon'};
ncobj.FORTRAN_format = ncchar('f10.4');
ncobj.units = ncchar('degree_east');
ncobj.type = ncchar('EVEN');
ncobj.epic_code = nclong(502);
ncobj.name = ncchar('LON');
ncobj.long_name = ncchar('LONGITUDE');
ncobj.generic_name = ncchar('lon');
ncobj(:) = theMeta.longitude;

cdfs{'lat'} = ncfloat('lat');  ncobj = cdfs{'lat'};
ncobj.FORTRAN_format = ncchar('F10.2');
ncobj.units = ncchar('degree_north');
ncobj.type = ncchar('EVEN');
ncobj.epic_code = nclong(500);
ncobj.name = ncchar('LAT');
ncobj.long_name = ncchar('LATITUDE');
ncobj.generic_name = ncchar('lat');
ncobj(:) = theMeta.latitude;

cdfs{'burst'} = nclong('time');
cdfs{'burst'}.FORTRAN_format = ncchar('F10.2');
cdfs{'burst'}.units = ncchar('counts');
cdfs{'burst'}.type = ncchar('EVEN');
cdfs{'burst'}.FillValue_ = ncfloat(1.e+035);

cdfs{'pctile'} = nclong('pctile');
cdfs{'pctile'}.FORTRAN_format = ncchar('F10.2');
cdfs{'pctile'}.units = ncchar('%');
cdfs{'pctile'}.type = ncchar('EVEN');
cdfs{'pctile'}.FillValue_ = ncfloat(1.e+035);
cdfs{'pctile'}(:) = miscMeta.pctiles;

    varname = 'vrange';
    cdfs{varname} = ncfloat('depth');
    buf = sprintf('Range from Transducer at %4.3f Mhz',theMeta.frequency(ixdr)./1e6);
    cdfs{varname}.long_name = ncchar(buf);
    buf = sprintf('%2.1f Mhz Range',theMeta.frequency(ixdr)./1e6);
    cdfs{varname}.name = ncchar(buf);
    cdfs{varname}.generic_name = ncchar(' ');
    cdfs{varname}.epic_code = 0;
    cdfs{varname}.units = ncchar('m');
    cdfs{varname}.valid_range = ncfloat([0 3]);
    cdfs{varname}.sensor_type = ncchar(theMeta.INST_TYPE);
    cdfs{varname}.initial_sensor_height = miscMeta.sensor_height(ixdr);
    cdfs{varname}.sensor_depth = theMeta.WATER_DEPTH - miscMeta.sensor_height(ixdr);
    %cdfs{varname}.serial = ncchar(cdfs.serial_number(:));
    cdfs{varname}.frequency = ncchar(sprintf('%2.1f Mhz',theMeta.frequency(ixdr)./1e6));
    cdfs{varname}.FORTRAN_format = ncchar('F10.2');
    cdfs{varname}.FillValue_ = ncfloat(1.e+035);
    cdfs{varname}(:) = theMeta.AbsBinRange(:,ixdr);

    varname = 'abs_trans_mean';
    cdfs{varname} = ncfloat('time','depth');
    buf = sprintf('Mean Amplitude at %4.3f Mhz',theMeta.frequency(ixdr)./1e6);
    cdfs{varname}.long_name = ncchar(buf);
    buf = sprintf('%2.1f Mhz Mean',theMeta.frequency(ixdr)./1e6);
    cdfs{varname}.name = ncchar(buf);
    cdfs{varname}.generic_name = ncchar(' ');
    cdfs{varname}.epic_code = 0;
    cdfs{varname}.units = ncchar('ABS units');
    cdfs{varname}.valid_range = ncfloat([0 1e16]);
    cdfs{varname}.sensor_type = ncchar(theMeta.INST_TYPE);
    cdfs{varname}.initial_sensor_height = miscMeta.sensor_height(ixdr);
    cdfs{varname}.sensor_depth = theMeta.WATER_DEPTH - miscMeta.sensor_height(ixdr);
    %cdfs{varname}.serial = ncchar(cdfs.serial_number(:));
    cdfs{varname}.frequency = ncchar(sprintf('%2.1f Mhz',theMeta.frequency(ixdr)./1e6));
    cdfs{varname}.FORTRAN_format = ncchar('F10.2');
    cdfs{varname}.FillValue_ = ncfloat(1.e+035);

    varname = 'abs_trans_mean_sq';
    cdfs{varname} = ncfloat('time','depth');
    buf = sprintf('Squared Mean Amplitude at %4.3f Mhz',theMeta.frequency(ixdr)./1e6);
    cdfs{varname}.long_name = ncchar(buf);
    buf = sprintf('%2.1f Mhz Mean',theMeta.frequency(ixdr)./1e6);
    cdfs{varname}.name = ncchar(buf);
    cdfs{varname}.generic_name = ncchar(' ');
    cdfs{varname}.epic_code = 0;
    cdfs{varname}.units = ncchar('ABS units^2');
    cdfs{varname}.valid_range = ncfloat([0 1e16]);
    cdfs{varname}.sensor_type = ncchar(theMeta.INST_TYPE);
    cdfs{varname}.initial_sensor_height = miscMeta.sensor_height(ixdr);
    cdfs{varname}.sensor_depth = theMeta.WATER_DEPTH - miscMeta.sensor_height(ixdr);
    %cdfs{varname}.serial = ncchar(cdfs.serial_number(:));
    cdfs{varname}.frequency = ncchar(sprintf('%2.1f Mhz',theMeta.frequency(ixdr)./1e6));
    cdfs{varname}.FORTRAN_format = ncchar('F10.2');
    cdfs{varname}.FillValue_ = ncfloat(1.e+035);

    varname = 'abs_trans_std';
    cdfs{varname} = ncfloat('time','depth');
    buf = sprintf('Amplitude Standard Deviation at %4.3f Mhz',theMeta.frequency(ixdr)./1e6);
    cdfs{varname}.long_name = ncchar(buf);
    buf = sprintf('%2.1f Mhz Mean',theMeta.frequency(ixdr)./1e6);
    cdfs{varname}.name = ncchar(buf);
    cdfs{varname}.generic_name = ncchar(' ');
    cdfs{varname}.epic_code = 0;
    cdfs{varname}.units = ncchar('ABS units');
    cdfs{varname}.valid_range = ncfloat([0 1e16]);
    cdfs{varname}.sensor_type = ncchar(theMeta.INST_TYPE);
    cdfs{varname}.initial_sensor_height = miscMeta.sensor_height(ixdr);
    cdfs{varname}.sensor_depth = theMeta.WATER_DEPTH - miscMeta.sensor_height(ixdr);
    %cdfs{varname}.serial = ncchar(cdfs.serial_number(:));
    cdfs{varname}.frequency = ncchar(sprintf('%2.1f Mhz',theMeta.frequency(ixdr)./1e6));
    cdfs{varname}.FORTRAN_format = ncchar('F10.2');
    cdfs{varname}.FillValue_ = ncfloat(1.e+035);

    varname = 'abs_trans_std_sq';
    cdfs{varname} = ncfloat('time','depth');
    buf = sprintf('Squared Amplitude Standard Deviation at %4.3f Mhz',theMeta.frequency(ixdr)./1e6);
    cdfs{varname}.long_name = ncchar(buf);
    buf = sprintf('%2.1f Mhz Mean',theMeta.frequency(ixdr)./1e6);
    cdfs{varname}.name = ncchar(buf);
    cdfs{varname}.generic_name = ncchar(' ');
    cdfs{varname}.epic_code = 0;
    cdfs{varname}.units = ncchar('ABS units^2');
    cdfs{varname}.valid_range = ncfloat([0 1e16]);
    cdfs{varname}.sensor_type = ncchar(theMeta.INST_TYPE);
    cdfs{varname}.initial_sensor_height = miscMeta.sensor_height(ixdr);
    cdfs{varname}.sensor_depth = theMeta.WATER_DEPTH - miscMeta.sensor_height(ixdr);
    %cdfs{varname}.serial = ncchar(cdfs.serial_number(:));
    cdfs{varname}.frequency = ncchar(sprintf('%2.1f Mhz',theMeta.frequency(ixdr)./1e6));
    cdfs{varname}.FORTRAN_format = ncchar('F10.2');
    cdfs{varname}.FillValue_ = ncfloat(1.e+035);

    varname = 'abs_trans_pctl';
    cdfs{varname} = ncfloat('time','depth','pctile');
    buf = sprintf('Amplitude Percentiles at %4.3f Mhz',theMeta.frequency(ixdr)./1e6);
    cdfs{varname}.long_name = ncchar(buf);
    buf = sprintf('%2.1f Mhz Mean',theMeta.frequency(ixdr)./1e6);
    cdfs{varname}.name = ncchar(buf);
    cdfs{varname}.generic_name = ncchar(' ');
    cdfs{varname}.epic_code = 0;
    cdfs{varname}.units = ncchar('ABS units');
    cdfs{varname}.valid_range = ncfloat([0 1e16]);
    cdfs{varname}.sensor_type = ncchar(theMeta.INST_TYPE);
    cdfs{varname}.initial_sensor_height = miscMeta.sensor_height(ixdr);
    cdfs{varname}.sensor_depth = theMeta.WATER_DEPTH - miscMeta.sensor_height(ixdr);
    %cdfs{varname}.serial = ncchar(cdfs.serial_number(:));
    cdfs{varname}.frequency = ncchar(sprintf('%2.1f Mhz',theMeta.frequency(ixdr)./1e6));
    cdfs{varname}.FORTRAN_format = ncchar('F10.2');
    cdfs{varname}.FillValue_ = ncfloat(1.e+035);

    varname = 'abs_trans_pctl_sq';
    cdfs{varname} = ncfloat('time','depth','pctile');
    buf = sprintf('Squared Amplitude Percentiles at %4.3f Mhz',theMeta.frequency(ixdr)./1e6);
    cdfs{varname}.long_name = ncchar(buf);
    buf = sprintf('%2.1f Mhz Mean',theMeta.frequency(ixdr)./1e6);
    cdfs{varname}.name = ncchar(buf);
    cdfs{varname}.generic_name = ncchar(' ');
    cdfs{varname}.epic_code = 0;
    cdfs{varname}.units = ncchar('ABS units^2');
    cdfs{varname}.valid_range = ncfloat([0 1e16]);
    cdfs{varname}.sensor_type = ncchar(theMeta.INST_TYPE);
    cdfs{varname}.initial_sensor_height = miscMeta.sensor_height(ixdr);
    cdfs{varname}.sensor_depth = theMeta.WATER_DEPTH - miscMeta.sensor_height(ixdr);
    %cdfs{varname}.serial = ncchar(cdfs.serial_number(:));
    cdfs{varname}.frequency = ncchar(sprintf('%2.1f Mhz',theMeta.frequency(ixdr)./1e6));
    cdfs{varname}.FORTRAN_format = ncchar('F10.2');
    cdfs{varname}.FillValue_ = ncfloat(1.e+035);

varname = 'sav';
cdfs{varname} = ncfloat('time');
buf = sprintf('sav at %4.3f Mhz',theMeta.frequency(ixdr)./1e6);
cdfs{varname}.long_name = ncchar(buf);
buf = sprintf('sav %2.1f Mhz',theMeta.frequency(ixdr)./1e6);
cdfs{varname}.name = ncchar(buf);
cdfs{varname}.generic_name = ncchar(' ');
cdfs{varname}.epic_code = 0;
cdfs{varname}.units = ncchar(' ');
cdfs{varname}.valid_range = ncfloat([0 1e16]);
cdfs{varname}.sensor_type = ncchar(theMeta.INST_TYPE);
cdfs{varname}.initial_sensor_height = miscMeta.sensor_height(ixdr);
cdfs{varname}.sensor_depth = theMeta.WATER_DEPTH - miscMeta.sensor_height(ixdr);
%cdfs{varname}.serial = ncchar(cdfs.serial_number(:));
cdfs{varname}.frequency = ncchar(sprintf('%2.1f Mhz',theMeta.frequency(ixdr)./1e6));
cdfs{varname}.FORTRAN_format = ncchar('F10.2');
cdfs{varname}.FillValue_ = ncfloat(1.e+035);

varname = 'ssd';
cdfs{varname} = ncfloat('time');
buf = sprintf('ssd at %4.3f Mhz',theMeta.frequency(ixdr)./1e6);
cdfs{varname}.long_name = ncchar(buf);
buf = sprintf('ssd %2.1f Mhz',theMeta.frequency(ixdr)./1e6);
cdfs{varname}.name = ncchar(buf);
cdfs{varname}.generic_name = ncchar(' ');
cdfs{varname}.epic_code = 0;
cdfs{varname}.units = ncchar(' ');
cdfs{varname}.valid_range = ncfloat([0 1e16]);
cdfs{varname}.sensor_type = ncchar(theMeta.INST_TYPE);
cdfs{varname}.initial_sensor_height = miscMeta.sensor_height(ixdr);
cdfs{varname}.sensor_depth = theMeta.WATER_DEPTH - miscMeta.sensor_height(ixdr);
%cdfs{varname}.serial = ncchar(cdfs.serial_number(:));
cdfs{varname}.frequency = ncchar(sprintf('%2.1f Mhz',theMeta.frequency(ixdr)./1e6));
cdfs{varname}.FORTRAN_format = ncchar('F10.2');
cdfs{varname}.FillValue_ = ncfloat(1.e+035);

endef(cdfs)

return

% ---------------------- Subfunction ------------------------------------ %
function calc_minmax_vals(cdf)

theVars = var(cdf);
for i = 1:length(theVars),
    if ~strcmp(ncnames(theVars{i}),'time') && ...
            ~strcmp(ncnames(theVars{i}),'time2')
        data = theVars{i}(:);
        N = ndims(data);
        if N==1   %1D data
            theVars{i}.minimum = ncfloat(min(data));
            theVars{i}.maximum = ncfloat(max(data));
        elseif N==2 %2D data
            [row, col] = size(data);
            for icol = 1:col
                mins(icol) = min(data(:,icol));
                maxs(icol) = max(data(:,icol));
            end
            theVars{i}.minimum = ncfloat(mins);
            theVars{i}.maximum = ncfloat(maxs);
            clear mins maxs
        elseif N==3 %3D data
            [row, col, lvl] = size(data);
            for jlvl = 1:lvl
                mins(:,jlvl) = min(data(:,:,jlvl))';
                maxs(:,jlvl) = max(data(:,:,jlvl))';
            end
            theVars{i}.minimum = ncfloat(mins);
            theVars{i}.maximum = ncfloat(maxs);
            clear mins maxs
        end
        clear data
    end
end
return