function theMeta = aqamat2cdf(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:  aqa2cdf(settings);
%
% 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: 842
%        Deployment_date: '04-dec-2007 23:54:00 GMT'
%          Recovery_date: '29-feb-2008 00:45:00 GMT'
%              longitude: -118.3542
%               latitude: 33.7100
%     magnetic_variation: 12.9000
%            WATER_DEPTH: 54
%            DATA_ORIGIN: 'USGS WHSC and PSC Sed. Trans. Groups'
%             EXPERIMENT: '2007-2008 PV BBL Experiment'
%            DESCRIPTION: [1x101 char]
%            Conventions: 'PMEL/EPIC'
%                PROJECT: 'USGS/EPA Palos Verdes'
%              INST_TYPE: 'Aquatec Aquascat Acoustic Backscatter Sensor'
%        input_directory: 'E:\PV_2007\842_B3_geoprobe\8427abss'
%       output_directory: 'E:\PV_2007\842_B3_geoprobe\8427abss'
%            outFileRoot: '8427abss'
%          sensor_height: [0.9900 0.9900 0.9900]
%         InterBurstTime: 60
%          BurstDuration: 840
%             goodbursts: [1 Inf]
%              ABSSmodel: 'AQUAscat1000R                   '
%       ABSSSerialNumber: '278-18'
%       ABSSBoardVersion: '3'
%      AbsTransducerName: [3x3 char]
%              frequency: [1035000 2535000 4700000]
%       AbsTxPulseLength: [1.3333e-005 1.3333e-005 1.3333e-005]
%               PingRate: 64
%               NumPings: 53760
%
% 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
%
% Dependencies:
%   readaqa_oleg.m (O. Mouraenko)
%   matpctile.m (C. Sherwood)
%   julian.m
%   gregorian.m

% TODO:  Read frequencies from 'ABS_info' variable if they're not specified
% in the user's metadata file
% TODO:  Add lat/lon variables? For longitude follow ADCP's lead, where
% longitude is negative and its epic key is 502.
% TODO:  Correct the units attribute
% TODO:  Correct the valid range attribute
% TODO: split into three files by channel
% TODO: consider these calculations
%   9.3	RangeCorrectAbsProfiles
%
% 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

more off

version = '0.0';

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

if exist(theMeta,'file'),
    disp(sprintf('Loading metadata from %s',theMeta))
    load(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?');
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

% global attributes which change
gattnames = {'MOORING', 'Deployment_date', 'Recovery_date', 'longitude', 'latitude',...
    'magnetic_variation', 'WATER_DEPTH', 'DATA_ORIGIN', 'EXPERIMENT','DESCRIPTION',...
    'PROJECT', 'WATER_DEPTH_NOTE','INST_TYPE','ABSSmodel'};
for ifield = 1:length(gattnames),
    if ~isfield(theMeta,gattnames{ifield}),
        disp(sprintf('Missing metadata for %s',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

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

% 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 = 500000000;
burstFileNum = 1; %initialize burst netcdf file #
lastBrec = 0; %initialize burst record #
AbsData = []; % this is to quiet mlint
AbsMean = [];

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,
        % define the output files
        %Create burst and statistic netCDF files
        [cdfb, burstFile] = defineAbsBurstCdf(theMeta, burstFileNum);
        [cdfs, statFile] = defineAbsStatCdf(theMeta);
        fprintf('\n')
        fprintf('There are %d bursts to write to netCDF\n',nMatFiles)
        fprintf('Writing burst data to the file %s\n',burstFile)
        fprintf('Writing statistics data to the file %s\n',statFile)

        %Extract dimensions.  These are constant for each
        %Aquascat file, so extract them once at the start.
        binDim = cdfb('bin');
        nBins = size(binDim,1);
        profDim = cdfb('profile');
        nProfiles = size(profDim,1);
        xducerDim = cdfb('xducer');
        nXducers = size(xducerDim,1);
        pctlDim = cdfs('pctile');
        nPctiles = size(pctlDim,1);

        %Some index variables are constant for each Aquascat
        %file, so write these only once at the beginning.
        if nXducers>1
            cdfb{'r'}(1:nBins,1:nXducers) = theMeta.AbsBinRange / 1000;  %Distance from the ABS, m
            cdfs{'r'}(1:nBins,1:nXducers) = theMeta.AbsBinRange / 1000;
        else
            cdfb{'r'}(1:nBins) = theMeta.AbsBinRange / 1000;  %Distance from the ABS, m
            cdfs{'r'}(1:nBins) = theMeta.AbsBinRange / 1000;
        end
        cdfs{'pctile'}(1:nPctiles) = pctiles; %Percentiles

    end

    disp(sprintf('reading burst %d, %s elapsed %f min',...
        BurstNumber, BurstTime, toc/60))

    %Write Aquascat data to the burst netCDF file
    % 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*(nProfiles-1));
    cdfb{'time'}(brec, :) = floor(profTime);
    cdfb{'time2'}(brec, :) = (profTime-floor(profTime)).*(24*3600*1000);
    if ~isequal(BurstNumber,n) %if burst #'s begin to repeat,
        cdfb{'burst'}(brec) = n;     %override them
    else
        cdfb{'burst'}(brec) = BurstNumber;
    end
    % shape of AbsData is [nbins x nprofiles x ntransducers]
    %cdfb{'abs_trans1'}(brec, 1:nProfiles, 1:nBins) = AbsData(:,:,1); %Amplitude at frequency 1
    %cdfb{'abs_trans2'}(brec, 1:nProfiles, 1:nBins) = AbsData(:,:,2); %Ampliude at frequency 2
    %cdfb{'abs_trans3'}(brec, 1:nProfiles, 1:nBins) = AbsData(:,:,3); %Ampliude at frequency 3
    for ixd = 1:3,
        varname = sprintf('abs_trans%d',ixd);
        cdfb{varname}(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
        filechk = dir(burstFile);
        bytesPerAqFile = filechk.bytes;
        nAqFilesPerBurstCdf = floor(maxBytes/bytesPerAqFile);
    end

    %write statistics to the statistic netCDF file
    %cdfs{'time'}(srec) = profTime(1);%use first time stamp in each burst
    cdfs{'time'}(srec) = floor(profStartTime);%use first time stamp in each burst
    %cdfs{'time2'}(srec) = (profTime(1)-floor(profTime(1))).*(24*3600*1000);
    cdfs{'time2'}(srec) = (profStartTime-floor(profStartTime)).*(24*3600*1000);
    if ~isequal(BurstNumber,n) %if burst #'s begin to repeat,
        cdfs{'burst'}(srec) = n;     %override them
    else
        cdfs{'burst'}(srec) = BurstNumber;
    end
    for ixd = 1:3,
        varname = sprintf('abs_trans%d_mean',ixd);
        cdfs{varname}(srec, :) = AbsMean(:,ixd);
        varname = sprintf('abs_trans%d_mean_sq',ixd);
        cdfs{varname}(srec, :) = gmean(AbsData(:,:,ixd)'.^2);
        varname = sprintf('abs_trans%d_std',ixd);
        cdfs{varname}(srec, :) = gstd(AbsData(:,:,ixd)');
        varname = sprintf('abs_trans%d_std_sq',ixd);
        cdfs{varname}(srec, :) = std(AbsData(:,:,ixd)'.^2);
        varname = sprintf('abs_trans%d_pctl',ixd);
        cdfs{varname}(srec, :, :) = ...
            matpctile(AbsData(:,:,ixd)', pctiles(:))';
        varname = sprintf('abs_trans%d_pctl_sq',ixd);
        cdfs{varname}(srec, :, :) = ...
            matpctile(AbsData(:,:,ixd)'.^2, pctiles)';
    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

    %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)==0
        lastBrec = cdfb{'burst'}(end);
        %Update the attributes of the current burst netCDF file
        fprintf('The burst file %s is full\n',burstFile)
        cdfb.burstnum = [];
        start_time = cdfb{'time'}(1,1) + cdfb{'time2'}(1,1)/3600/1000/24;
        cdfb.start_time = ncchar(datestr(gregorian(start_time)));
        stop_time = cdfb{'time'}(end,end) + cdfb{'time2'}(end,end)/3600/1000/24;
        cdfb.stop_time = ncchar(datestr(gregorian(stop_time)));
        hist = cdfb.history(:);
        hist_new = ['Converted to netCDF via Matlab by %s V ', mfilename, version,'; ',hist];
        cdfb.history = ncchar(hist_new);
        close(cdfb);
        %Create a new burst netCDF file for subsequent bursts
        burstFileNum = burstFileNum + 1;
        [cdfb, burstFile] = defineAbsBurstCdf(theMeta, burstFileNum);
        fprintf('Now writing burst data to the file %s\n',burstFile)
        cdfb{'r'}(1:nBins,1:nXducers) = theMeta.AbsBinRange / 1000;  %Distance from the ABS, m
    end


end

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

%update the attributes of the statistics netCDF file and
%calculate min/max values
cdfs.burstnum = [];
start = cdfs{'time'}(1) + cdfs{'time2'}(1)/3600/1000/24;
cdfs.start_time = ncchar(datestr(gregorian(start)));
stop = cdfs{'time'}(end) + cdfs{'time2'}(end)/3600/1000/24;
cdfs.stop_time = ncchar(datestr(gregorian(stop)));
delta_t = mean(diff(cdfs{'time'}(:) + cdfs{'time2'}(:)/3600/1000/24));
cdfs.DELTA_T = ncchar([num2str(delta_t*24*3600),' s']);
hist = cdfs.history(:);
hist_new = ['Converted to netCDF via Matlab by aqa2cdf.m V ',version,'; ',hist];
cdfs.history = ncchar(hist_new);
calc_minmax_vals(cdfs);
close(cdfs);
fprintf('Finished writing all bursts. %6.2f minutes elapsed\n',toc/60)
%end

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

outputDir = theMeta.output_directory;
burstFileNum = num2str(burstFileNum);
burstFile = fullfile(outputDir,[theMeta.outFileRoot,burstFileNum,'b.cdf']);

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'};
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>1
            for ii = 1:nRows
                temp = [temp,' ',theFieldDef(ii,1:nCols)];
            end
            theFieldDef = temp;
        end
        eval(['cdfb.',metaFields{i},' = ncchar(theFieldDef);'])
    else
        eval(['cdfb.',metaFields{i},'= ncfloat(theFieldDef);'])
    end
end

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

% define dimensions
% cdfb('burst') = 0; %unlimited
% cdfb('profile') = theMeta.nProfiles;
% cdfb('bin') = theMeta.nBins;
% cdfb('xducer') = theMeta.nXducers;
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{'time'}.FillValue_ = ncfloat(1.00000004091848e+035);

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);
%cdfb{'time2'}.FillValue_ = ncfloat(1.00000004091848e+035);

% we need depth.
% TODO Not too sure how to handle when transducers are at different depths
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(:,1)+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);

for idx = 1:theMeta.nXducers,
    varname = sprintf('vrange%d',idx);
    cdfb{varname} = ncfloat('depth');
    buf = sprintf('Range from Transducer at %4.3f Mhz',theMeta.frequency(idx)./6);
    cdfb{varname}.long_name = ncchar(buf);
    buf = sprintf('%2.1f Mhz Range',theMeta.frequency(idx)./6);
    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(idx);
    cdfb{varname}.sensor_depth = theMeta.WATER_DEPTH - miscMeta.sensor_height(idx);
    %cdfb{varname}.serial = ncchar(cdfb.serial_number(:));
    cdfb{varname}.frequency = ncchar(sprintf('%2.1f Mhz',theMeta.frequency(idx)./6));
    cdfb{varname}.FORTRAN_format = ncchar('F10.2');
    cdfb{varname}.FillValue_ = ncfloat(1.e+035);
    cdfb{varname}(:) = theMeta.AbsBinRange(:,idx);
end

for idx = 1:theMeta.nXducers,
    varname = sprintf('abs_trans%d',idx);
    cdfb{varname} = ncfloat('time','profile','depth');
    buf = sprintf('Amplitude at %4.3f Mhz',theMeta.frequency(idx)./6);
    cdfb{varname}.long_name = ncchar(buf);
    buf = sprintf('%2.1f Mhz Mean',theMeta.frequency(idx)./6);
    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(idx);
    cdfb{varname}.sensor_depth = theMeta.WATER_DEPTH - miscMeta.sensor_height(idx);
    %cdfb{varname}.serial = ncchar(cdfb.serial_number(:));
    cdfb{varname}.frequency = ncchar(sprintf('%2.1f Mhz',theMeta.frequency(idx)./6));
    cdfb{varname}.FORTRAN_format = ncchar('F10.2');
    cdfb{varname}.FillValue_ = ncfloat(1.e+035);
end

endef(cdfb)

return

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

outputDir = theMeta.output_directory;
statFile = fullfile(outputDir,[theMeta.outFileRoot,'s.cdf']);

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'};
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>1
            for ii = 1:nRows
                temp = [temp,' ',theFieldDef(ii,1:nCols)];
            end
            theFieldDef = temp;
        end
        eval(['cdfs.',metaFields{i},' = ncchar(theFieldDef);'])
    else
        eval(['cdfs.',metaFields{i},'= ncfloat(theFieldDef);'])
    end
end

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

% define dimensions
% cdfs('burst') = 0; %unlimited
% cdfs('bin') = theMeta.nBins;
% cdfs('pctile') = miscMeta.nPctiles;
% cdfs('xducer') = theMeta.nXducers;
cdfs('time') = 0; %unlimited
cdfs('depth') = theMeta.nBins;
cdfs('lat') = 1;
cdfs('lon') = 1;
cdfs('pctile') = miscMeta.nPctiles;
%cdfs('xducer') = theMeta.nXducers;

% 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{'time'}.FillValue_ = ncfloat(1.e+035);

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);
%cdfs{'time2'}.FillValue_ = ncfloat(1.e+035);

% we need depth.
% TODO Not too sure how to handle when transducers are at different depths
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.AbsBinRange(:,1)+theMeta.WATER_DEPTH(1);

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;

for idx = 1:theMeta.nXducers,
    varname = sprintf('vrange%d',idx);
    cdfs{varname} = ncfloat('depth');
    buf = sprintf('Range from Transducer at %4.3f Mhz',theMeta.frequency(idx)./6);
    cdfs{varname}.long_name = ncchar(buf);
    buf = sprintf('%2.1f Mhz Range',theMeta.frequency(idx)./6);
    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(idx);
    cdfs{varname}.sensor_depth = theMeta.WATER_DEPTH - miscMeta.sensor_height(idx);
    %cdfs{varname}.serial = ncchar(cdfs.serial_number(:));
    cdfs{varname}.frequency = ncchar(sprintf('%2.1f Mhz',theMeta.frequency(idx)./6));
    cdfs{varname}.FORTRAN_format = ncchar('F10.2');
    cdfs{varname}.FillValue_ = ncfloat(1.e+035);
    cdfs{varname}(:) = theMeta.AbsBinRange(:,idx);
end

for idx = 1:theMeta.nXducers,
    varname = sprintf('abs_trans%d_mean',idx);
    cdfs{varname} = ncfloat('time','depth');
    buf = sprintf('Mean Amplitude at %4.3f Mhz',theMeta.frequency(idx)./6);
    cdfs{varname}.long_name = ncchar(buf);
    buf = sprintf('%2.1f Mhz Mean',theMeta.frequency(idx)./6);
    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(idx);
    cdfs{varname}.sensor_depth = theMeta.WATER_DEPTH - miscMeta.sensor_height(idx);
    %cdfs{varname}.serial = ncchar(cdfs.serial_number(:));
    cdfs{varname}.frequency = ncchar(sprintf('%2.1f Mhz',theMeta.frequency(idx)./6));
    cdfs{varname}.FORTRAN_format = ncchar('F10.2');
    cdfs{varname}.FillValue_ = ncfloat(1.e+035);
end

for idx = 1:theMeta.nXducers,
    varname = sprintf('abs_trans%d_mean_sq',idx);
    cdfs{varname} = ncfloat('time','depth');
    buf = sprintf('Squared Mean Amplitude at %4.3f Mhz',theMeta.frequency(idx)./6);
    cdfs{varname}.long_name = ncchar(buf);
    buf = sprintf('%2.1f Mhz Mean',theMeta.frequency(idx)./6);
    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(idx);
    cdfs{varname}.sensor_depth = theMeta.WATER_DEPTH - miscMeta.sensor_height(idx);
    %cdfs{varname}.serial = ncchar(cdfs.serial_number(:));
    cdfs{varname}.frequency = ncchar(sprintf('%2.1f Mhz',theMeta.frequency(idx)./6));
    cdfs{varname}.FORTRAN_format = ncchar('F10.2');
    cdfs{varname}.FillValue_ = ncfloat(1.e+035);
end

for idx = 1:theMeta.nXducers,
    varname = sprintf('abs_trans%d_std',idx);
    cdfs{varname} = ncfloat('time','depth');
    buf = sprintf('Amplitude Standard Deviation at %4.3f Mhz',theMeta.frequency(idx)./6);
    cdfs{varname}.long_name = ncchar(buf);
    buf = sprintf('%2.1f Mhz Mean',theMeta.frequency(idx)./6);
    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(idx);
    cdfs{varname}.sensor_depth = theMeta.WATER_DEPTH - miscMeta.sensor_height(idx);
    %cdfs{varname}.serial = ncchar(cdfs.serial_number(:));
    cdfs{varname}.frequency = ncchar(sprintf('%2.1f Mhz',theMeta.frequency(idx)./6));
    cdfs{varname}.FORTRAN_format = ncchar('F10.2');
    cdfs{varname}.FillValue_ = ncfloat(1.e+035);
end

for idx = 1:theMeta.nXducers,
    varname = sprintf('abs_trans%d_std_sq',idx);
    cdfs{varname} = ncfloat('time','depth');
    buf = sprintf('Squared Amplitude Standard Deviation at %4.3f Mhz',theMeta.frequency(idx)./6);
    cdfs{varname}.long_name = ncchar(buf);
    buf = sprintf('%2.1f Mhz Mean',theMeta.frequency(idx)./6);
    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(idx);
    cdfs{varname}.sensor_depth = theMeta.WATER_DEPTH - miscMeta.sensor_height(idx);
    %cdfs{varname}.serial = ncchar(cdfs.serial_number(:));
    cdfs{varname}.frequency = ncchar(sprintf('%2.1f Mhz',theMeta.frequency(idx)./6));
    cdfs{varname}.FORTRAN_format = ncchar('F10.2');
    cdfs{varname}.FillValue_ = ncfloat(1.e+035);
end

for idx = 1:theMeta.nXducers,
    varname = sprintf('abs_trans%d_pctl',idx);
    cdfs{varname} = ncfloat('time','depth','pctile');
    buf = sprintf('Amplitude Percentiles at %4.3f Mhz',theMeta.frequency(idx)./6);
    cdfs{varname}.long_name = ncchar(buf);
    buf = sprintf('%2.1f Mhz Mean',theMeta.frequency(idx)./6);
    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(idx);
    cdfs{varname}.sensor_depth = theMeta.WATER_DEPTH - miscMeta.sensor_height(idx);
    %cdfs{varname}.serial = ncchar(cdfs.serial_number(:));
    cdfs{varname}.frequency = ncchar(sprintf('%2.1f Mhz',theMeta.frequency(idx)./6));
    cdfs{varname}.FORTRAN_format = ncchar('F10.2');
    cdfs{varname}.FillValue_ = ncfloat(1.e+035);
end

for idx = 1:theMeta.nXducers,
    varname = sprintf('abs_trans%d_pctl_sq',idx);
    cdfs{varname} = ncfloat('time','depth','pctile');
    buf = sprintf('Squared Amplitude Percentiles at %4.3f Mhz',theMeta.frequency(idx)./6);
    cdfs{varname}.long_name = ncchar(buf);
    buf = sprintf('%2.1f Mhz Mean',theMeta.frequency(idx)./6);
    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(idx);
    cdfs{varname}.sensor_depth = theMeta.WATER_DEPTH - miscMeta.sensor_height(idx);
    %cdfs{varname}.serial = ncchar(cdfs.serial_number(:));
    cdfs{varname}.frequency = ncchar(sprintf('%2.1f Mhz',theMeta.frequency(idx)./6));
    cdfs{varname}.FORTRAN_format = ncchar('F10.2');
    cdfs{varname}.FillValue_ = ncfloat(1.e+035);
end

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