% sg2nc   SeaGuage to nc converion tool
%
% function sg2nc(wbFile, tidFile, wasFile, outFileRoot, metadata)
% sgRoot = root name of the seabird data files, *.tid or *.wb
%
% *.wb file structure
% SBE 26plus
% *    0   197513221    1.00  512
%          14.260035     14.260035     14.258990     14.259408
%          14.258990     14.257946     14.258363     14.257946
%          14.258990     14.259408     14.258990     14.258990
%          14.258990     14.257946     14.258990     14.258990
%
% *.tid file structure
%                                Meters      Deg C
%      1   04/04/2006 23:58:00    -0.2318     13.577
%      2   04/05/2006 00:13:00    -0.2316     13.407
%      3   04/05/2006 00:28:00    -0.2359     13.229
%      4   04/05/2006 00:43:00    -0.2404     13.044
% ------------- Step II what is the metadata?
%
% % here is the metadata structure for this deployment
% % if you use an empty char string denoted by ''
% metadata.mooring_number = '8183'; %four digit NNNL, NNN mooring + L logger
% metadata.deployment_start = '09-apr-2006 20:41:00'; % on bottom time
% metadata.deployment_end = '21-jun-2006 11:50:31'; % the out of water time
% metadata.lon = -73.6888; % decimal degrees, West = negative
% metadata.lat = 40.1615; % decimal degrees, South = negative
% % 14� 58' W changing by 0� 3' E/year from
% % http://www.ngdc.noaa.gov/seg/geomag/jsp/Declination.jsp
% metadata.declination = -13.22; % degrees west is negative
% metadata.water_depth = 64; % m or 109.6 ft (DO NOT OMIT!)
% metadata.press.serial = '1099'; % serial number of the SeaGuage
% metadata.press.height = 0.51; % height of pressure port above bottom, m
% % SBE37's record data in real units and do not need applied cals
% metadata.metafile = mfilename('fullpath');
% metadata.metafile_date = datestr(now);
% metadata.metafile_version = '0.0';
% metadata.metafile_author = 'MM';
% % if the following fields are missing, UNKOWN will be used
% metadata.origin = 'USGS/WHFC'; % with collaborator's data, could be USC, etc.
% metadata.experiment = 'NY Bight';
% metadata.project = 'NY Bight'; % might also use OFA funding agency, such as MWRA, EPA, WCMG
% metadata.comment = '';
% metadata.description = 'Site F wave and tide guage';

function sg2nc(wbFile, tidFile, wasFile, wtsFile, outFileRoot, metadata)


% do the tide file
if ~isempty(tidFile),

    disp(sprintf('Reading the tide data in %s', tidFile))
    fid = fopen(tidFile,'r');
    if fid < 0,
        disp(sprintf('Unable to open %s',tidFile))
        return
    end
    thePath = fileparts(tidFile);

    % skip header line
    data = fgetl(fid);
    % see if Salinity is present
    if findstr(data,'PSU'),
        disp('Salinity data found')
        %1   11/13/92 10:27:16   14.8125    22.102  3.55682    23.909
        data = fscanf(fid,'%d %d/%d/%d %d:%d:%d %f %f %f %f\n',[11,Inf]);
    else
        disp('Salinity data not found')
        %1   11/13/92 10:27:16   14.8125    22.102   
        data = fscanf(fid,'%d %d/%d/%d %d:%d:%d %f %f\n',[9,Inf]);
    end
    fclose(fid);

    ncFile = fullfile(thePath, [outFileRoot,'tid.nc']);
    nc = definetid(ncFile, length(data), metadata);

    disp(sprintf('Writing tide data to the netCDF %s',ncFile))
    metadata.instrument = ''; 
    setup_metadata(nc, metadata)

    % TODO set up the variable attributes

    if ~isempty(data);
        tj = julian([data(4,:)' data(2,:)' data(3,:)' data(5,:)' data(6,:)' data(7,:)']);
        nc{'time'}(:) = floor(tj);
        nc{'time2'}(:) = (tj-floor(tj)).*(1000*3600*24);
        nc{'burst'}(:,1,1,1) = data(1,:); % tide measurement number
        nc{'T_20'}(:,1,1,1) = data(9,:); % temperature in C
        nc{'P_1'}(:,1,1,1) = data(8,:)*0.0689746*10; % pressure, psia converted to db
        if ~isempty(metadata.cond.serial),
            nc{'C_51'}(:,1,1,1) = data(10,:);
            nc{'S_41'}(:,1,1,1) = data(11,:);
        end
    else
        disp('Did not read any data in tide file')
    end

    tj = nc{'time'}(:)+nc{'time2'}(:)./(1000*3600*24);
    dt = mean(diff(tj));
    nc.DELTA_T = ncchar(sprintf('%d',dt));
    nc.history = ncchar('converted by sg2nc.m');
    add_vardesc(nc);
    add_minmaxvalues(nc);
    add_fillvalues(nc,1E35);
    add_depthvalues(nc,metadata);

    close(nc);
end

%% do the burst file
if ~isempty(wbFile),

    disp(sprintf('Reading the wave burst data in %s', wbFile))

    fid = fopen(wbFile,'r');
    if fid < 0,
        disp(sprintf('Unable to open %s',wbFile))
        return
    end
    thePath = fileparts(wbFile);

    % get the number of samples per burst
    header1 = fgetl(fid);
    if strcmp(header1, 'SBE 26plus'),
        metadata.instrument = 'SBE 26plus';
        header1 = fscanf(fid,'%c    %d %d %f %d',5);
    else
        metadata.instrument = 'SBE 26';
        header1 = sscanf(header1,'%c    %d %d %f %d',5);
    end
    metadata.nsamples = header1(5);
    metadata.dt = header1(4);
    frewind(fid);

    % now define the file
    ncFile = fullfile(thePath, [outFileRoot,'wb.nc']);
    nc = definewb(ncFile, metadata);

    disp(sprintf('Writing wave burst data to the netCDF %s',ncFile))
    setup_metadata(nc, metadata)

    % TODO set up the variable attributes

    % skip header line only in 26plus
    if strcmp(metadata.instrument, 'SBE 26plus'); fgetl(fid); end
    ncidx = 1;
    tic

    while ftell(fid) > 0,
        % find the beginning of the next wave burst by hunting for the *
        findstar(fid);
        % read a burst header
        % *    0   197513221    1.00  512
        %          14.260035     14.260035     14.258990     14.259408
        [header1, count] = fscanf(fid,'%d %d %f %d',4);
        if count ~= 4, break; end
        tj = julian(2000,1,1)+header1(3)./(3600*24);
        tg = gregorian(tj);
        % expand time so there is a time stamp for every sample assume that a burst
        % is never going to be longer than a day
        time1 = ones(1,metadata.nsamples).*julian(tg(1),tg(2),tg(3)); % julian day
        time2 = (0:metadata.dt:((metadata.nsamples*metadata.dt)-...
            metadata.dt)).*1000; % in msec from start of burst
        % time2 is the start of burst in milliseconds from the start of the julian day
        time2 = time2 + tg(4).*(3600*1000)+... % hours to milliseconds
            tg(5).*(60*1000)+... % minutes to milliseconds
            tg(6).*(1000); % seconds to milliseconds
        % but account for those bursts that extend over the day boundary
        rolloveridx = find(time2 >= 1000*60*60*24); % find amounts greater than the # of milliseconds in the day
        if ~isempty(rolloveridx),
            time2(rolloveridx) = time2(rolloveridx) - 1000*60*60*24; % we will add a day to time 1 for these
            time1(rolloveridx) = time1(rolloveridx) + 1; % these are technically in the next day
        end
        nc{'time'}(ncidx,:) = time1;
        nc{'time2'}(ncidx,:) = time2;
        nc{'burst'}(ncidx) = header1(1);
        data = fscanf(fid,'%f',metadata.nsamples);
        nc{'P_4022'}(ncidx,:) = data.*0.0689746.*10; % measured pressure, psia converted to db
        fgetl(fid); % need this to skip a \n that the fscanf doesn't read
        ncidx = ncidx+1;
        if ~rem(ncidx,100),
            disp(sprintf('Finished burst #%d at %s %4.2f min elapsed',...
                ncidx,datestr(datenum(gregorian(tj))),toc./60)),
        end

    end

    tj = nc{'time'}(:)+nc{'time2'}(:)./(1000*3600*24);
    dt = mean(diff(tj));
    nc.DELTA_T = ncchar(sprintf('%d',dt));
    nc.history = ncchar('converted by sg2nc.m');
    add_vardesc(nc);
    add_minmaxvalues(nc);
    add_fillvalues(nc,1E35);
    add_depthvalues(nc,metadata);
    
    fclose(fid);
    close(nc);
end

%% do the burst auto-spectrum statistics file
if ~isempty(wasFile),

    disp(sprintf('Reading the burst auto-spectrum statistics data in %s', wasFile))

    fid = fopen(wasFile,'r');
    if fid < 0,
        disp(sprintf('Unable to open %s',wasFile))
        return
    end
    thePath = fileparts(wasFile);

    % get the number of samples per burst
    header1 = fgetl(fid);
    if strcmp(header1, 'SBE 26plus'),
        metadata.instrument = 'SBE 26plus';
        header1 = fscanf(fid,'%c %d %d %f %d %d %f %f %f %d %f %f',12);
    else
        metadata.instrument = 'SBE 26';
        header1 = fscanf(fid,'%c %d %d %f %d %d %f %f %f %d %f %f',12);
    end
    metadata.dt = header1(4); % wave integration time
    metadata.nsamples = header1(5); % samples in burst
    metadata.estperband = header1(6); % spectral estimates per band
    header2 = fscanf(fid,'%d %f %f %f %f %f %f',7);
    % number of frequency bands calculated - this can vary!
    metadata.maxfbands = header2(1); % I hope!
    metadata.f1 = header2(2);
    metadata.df = header2(3);
    frewind(fid);

    % now define the file
    ncFile = fullfile(thePath, [outFileRoot,'was.nc']);
 
    nc = definewas(ncFile, metadata);

    disp(sprintf('Writing burst auto-spectrum data to the netCDF %s',ncFile))
    setup_metadata(nc, metadata)

    % define the frequencies
    nc{'frequency'}(:) = metadata.f1:metadata.df:(metadata.f1+metadata.df*(metadata.maxfbands-1));

    % TODO set up the variable attributes

    ncidx = 1;
    tic

    while ~feof(fid),
        % find the beginning of the next wave burst by hunting for the *
        findstar(fid);
        %         % read a burst header
        % SBE 26plus
        % *    0   197513221    1.00  512    5     0.510    0.000  1025.627    90    0.546    2.538
        %     51    5.859375e-003    9.765625e-003    1.5395e-007    1.5484e-003    1.7067e+002    1.5695e-003
        %      2.214949e-006   1.257661e-006   1.496020e-006   5.784202e-007
        %      1.444130e-006   4.891811e-007   7.384434e-007   6.027650e-007
        [header1, count] = fscanf(fid,'%d %d %f %d %d %f %f %f %d %f %f',11);
        if count ~= 11,  disp('stopped at header1'); break; end
        [header2, count] = fscanf(fid,'%d %f %f %f %f %f %f',7);
        if count ~= 7,  disp('stopped at header2'); break; end
        tj = julian(2000,1,1)+header1(2)./(3600*24);
        % make the time stamp the middle of the burst
        tjend = tj+(metadata.nsamples*metadata.dt/(3600*24));
        tjmid = tj+(tjend-tj)/2;
        nc{'time'}(ncidx) = floor(tjmid);
        nc{'time2'}(ncidx) = (tjmid-floor(tjmid))*(1000*3600*24);
        nc{'burst'}(ncidx) = header1(1);
        nc{'D_3'}(ncidx) = header1(6); % water depth
        nc{'hght_18'}(ncidx) = header1(7); % pressure sensor depth
        nc{'RHO_84'}(ncidx) = header1(8); % density kg/m3
        nc{'chi2'}(ncidx) = header1(9); % chi squared confidence interval %
        nc{'chi2upper'}(ncidx) = header1(10); % multiplier for chi squared upper bound
        nc{'chi2lower'}(ncidx) = header1(11); % multiplier for chi squared lower bound
        nc{'nfreq'}(ncidx) = header2(1); % number of frequency bands calculated - this can vary!
        nc{'tvar'}(ncidx) = header2(4); % Total Variance, m^2
        nc{'tjoul'}(ncidx) = header2(5); % Total Energy (joules)
        nc{'wp_sig'}(ncidx) = header2(6); % Significant Wave Period (s)
        nc{'wh_4061'}(ncidx) = header2(7); % Significant wave period (m) = 4xsqrt(total variance)
        [data, count] = fscanf(fid,'%f',header2(1));
        if count ~= header2(1), disp('stopped at data'); break; end
        nc{'pspec'}(ncidx,1:header2(1)) = data;
        ncidx = ncidx+1;
        if ~rem(ncidx,100),
            disp(sprintf('Finished burst #%d at %s %4.2f min elapsed',...
                ncidx,datestr(datenum(gregorian(tj))),toc./60)),
        end

    end

    tj = nc{'time'}(:)+nc{'time2'}(:)./(1000*3600*24);
    dt = mean(diff(tj));
    nc.DELTA_T = ncchar(sprintf('%d',dt));
    nc.history = ncchar('converted by sg2nc.m');
    add_vardesc(nc);
    add_minmaxvalues(nc);
    add_fillvalues(nc,1E35);
    add_depthvalues(nc,metadata);
    
    fclose(fid);
    close(nc);
end

%% do the surface wave time series statistics file
if ~isempty(wtsFile),

    disp(sprintf('Reading the burst auto-spectrum statistics data in %s', wtsFile))

    fid = fopen(wtsFile,'r');
    if fid < 0,
        disp(sprintf('Unable to open %s',wtsFile))
        return
    end
    thePath = fileparts(wtsFile);

    % get the number of samples per burst
    header = fgetl(fid);
    if strcmp(header, 'SBE 26plus'),
        metadata.instrument = 'SBE 26plus';
        header = fscanf(fid,'%c %d %d %f %d %d %f %f %f',9);
    else
        metadata.instrument = 'SBE 26';
        header = fscanf(fid,'%c %d %d %f %d %d %f %f %f',9);
    end
    metadata.dt = header(4); % wave integration time
    metadata.nsamples = header(5); % samples in burst
    frewind(fid);

    % now define the file
    ncFile = fullfile(thePath, [outFileRoot,'wts.nc']);
 
    nc = definewts(ncFile, metadata);

    disp(sprintf('Writing surface wave time series data to the netCDF %s',ncFile))
    setup_metadata(nc, metadata)

    % TODO set up the wts variable attributes

    ncidx = 1;
    tic

    while ~feof(fid),
        % find the beginning of the next wave burst by hunting for the *
        findstar(fid);
        %         % read a burst header
        % SBE 26plus
        % *    0  197513221   1.00  512    0      0.510      0.000   1025.627
        %      1.937564e-007    1.948795e-003    0.000000e+000    0.000000e+000
        %      2.431568e-003    0.000000e+000    0.000000e+000    0.000000e+000    0.000000e+000        
        [header, count] = fscanf(fid,'%d %d %f %d %d %f %f %f',8);
        if count ~= 8,  disp('stopped at header'); break; end
        tj = julian(2000,1,1)+header(2)./(3600*24);
        % make the time stamp the middle of the burst
        tjend = tj+(metadata.nsamples*metadata.dt/(3600*24));
        tjmid = tj+(tjend-tj)/2;
        nc{'time'}(ncidx) = floor(tjmid);
        nc{'time2'}(ncidx) = (tjmid-floor(tjmid))*(1000*3600*24);
        nc{'burst'}(ncidx) = header(1);
        nc{'nw'}(ncidx) = header(5); % number of waves
        nc{'D_3'}(ncidx) = header(6); % water depth
        nc{'hght_18'}(ncidx) = header(7); % pressure sensor depth
        %nc{'RHO_84'}(ncidx) = header(8); % density kg/m3
        [data, count] = fscanf(fid,'%f',9);
        if count ~= 9, disp('stopped at data'); break; end
        nc{'awh'}(ncidx) = data(3); % average wave height
        nc{'wp_4060'}(ncidx) = data(4); % average wave period
        nc{'mwh_4064'}(ncidx) = data(5); % maximum wave height
        nc{'wh_4061'}(ncidx) = data(6); % significant wave height
        nc{'dwp_4063'}(ncidx) = data(7); % dominant wave period
        nc{'awh10'}(ncidx) = data(8); % average wave height of largest 1/10 waves
        nc{'awh100'}(ncidx) = data(9); % average wave height of largest 1/100 waves
        ncidx = ncidx+1;
        if ~rem(ncidx,100),
            disp(sprintf('Finished burst #%d at %s %4.2f min elapsed',...
                ncidx,datestr(datenum(gregorian(tj))),toc./60)),
        end

    end

    tj = nc{'time'}(:)+nc{'time2'}(:)./(1000*3600*24);
    dt = mean(diff(tj));
    nc.DELTA_T = ncchar(sprintf('%d',dt));
    nc.history = ncchar('converted by sg2nc.m');
    add_vardesc(nc);
    add_minmaxvalues(nc);
    add_fillvalues(nc,1E35);
    add_depthvalues(nc,metadata);
    
    fclose(fid);
    close(nc);
end

return

%% --------------------------------
% set up the metadata
function setup_metadata(nc, metadata)
nc.CREATION_DATE = ncchar(datestr(now));
nc.DATA_TYPE = ncchar('TIME');
nc.DATA_SUBTYPE = ncchar(' ');
nc.DATA_ORIGIN = ncchar('USGS/WHFC');
nc.COORD_SYSTEM = ncchar('GEOGRAPHICAL');
nc.WATER_MASS = ncchar('?');
nc.EXPERIMENT = ncchar(metadata.experiment);
nc.PROJECT = ncchar(metadata.project);
nc.MOORING = ncchar(metadata.mooring_number);
nc.COMPOSITE = nclong(0);
nc.DATA_CMNT = ncchar(metadata.comment);
nc.POS_CONST = nclong(0);
nc.DEPTH_CONST = nclong(0);
nc.WATER_DEPTH = nclong(metadata.water_depth);
nc.DESCRIPT = ncchar(metadata.description);
nc.DRIFTER = nclong(0);
nc.FILL_FLAG = nclong(1E35);
nc.start_time = ncchar(metadata.deployment_start);
nc.stop_time = ncchar(metadata.deployment_end);
nc.latitude = ncfloat(metadata.lat);
nc.longitude = ncfloat(metadata.lon);
nc.magnetic_variation = ncfloat(metadata.declination);
nc.water_depth = ncfloat(metadata.water_depth);
nc.CREATION_DATE = ncchar(datestr(now));
nc.INST_TYPE = ncchar(metadata.instrument);
return

%% ------------------------
% ADD_VARDESC construct the EPIC attribute var_desc on the fly
function add_vardesc(nc)
theVars = var(nc);
theDims = dim(nc);
var_desc = '';
for ivar = 1:length(theVars),
    adim = 0;
    for idim = 1:length(theDims),
        if strcmp(ncnames(theVars{ivar}),ncnames(theDims{idim})) ...
                || strcmp(ncnames(theVars{ivar}),'time2'),
            adim = 1;
        end
    end
    if ~adim,
        var_desc = sprintf('%s:%s',var_desc,theVars{ivar}.name(:));
    end
end
var_desc = var_desc(2:end);
nc.VAR_DESC = ncchar(var_desc);
disp(sprintf('Writing variables %s',var_desc));
return

%% ADD_MINMAXVALUES calculate min and max for variables in a netCDF file and
% add as an attribute.

function add_minmaxvalues(nc)
theVars = var(nc);
for i = 1:length(theVars),
    if ~strcmp(ncnames(theVars{i}),'time') && ~strcmp(ncnames(theVars{i}),'time2'),
        data = theVars{i}(:);
        [row, col] = size(data);
        %mins = zeros(1,col);
        %maxs = zeros(1,col);
        if col == 1,
            theVars{i}.minimum = ncfloat(gmin(data));
            theVars{i}.maximum = ncfloat(gmax(data));
        elseif col == 2,
            theVars{i}.minimum = ncfloat(gmin(gmin(data)));
            theVars{i}.maximum = ncfloat(gmax(gmax(data)));
        elseif col > 2, % then these are amp and cor or some other 3d var
            for icol = 1:col,
                mins(icol) = gmin(data(:,icol));
                maxs(icol) = gmax(data(:,icol));
            end
            theVars{i}.minimum = ncfloat(mins);
            theVars{i}.maximum = ncfloat(maxs);
        end
    end
end
return

%% ADD_FILLVALUES add the requisite fill value to all netCDF var EXCEPT time & time 2

function add_fillvalues(nc,fill)
theVars = var(nc);  % variable objects
for iVar = 1:length(theVars),
    % time does not get a fill value!
    if ~(strcmp(char(ncnames(theVars{iVar})),'time') || ...
            strcmp(char(ncnames(theVars{iVar})),'time2')),
        switch datatype(theVars{iVar}),
            case 'float'
                theVars{iVar}.FillValue_ = ncfloat(fill);
            case 'double'
                theVars{iVar}.FillValue_ = ncdouble(fill);
            case 'long'
                theVars{iVar}.FillValue_ = nclong(fill);
        end
    end
end
return

%% ADD_DEPTHVALUES add the requisite fill value to all netCDF var EXCEPT time & time 2

function add_depthvalues(nc,metadata)
theVars = var(nc);  % variable objects
for iVar = 1:length(theVars),
    % time does not get a fill value!
    if ~( strcmp(char(ncnames(theVars{iVar})),'time') || ...
            strcmp(char(ncnames(theVars{iVar})),'time2') || ...
            strcmp(char(ncnames(theVars{iVar})),'depth') || ...
            strcmp(char(ncnames(theVars{iVar})),'frequency') || ...
            strcmp(char(ncnames(theVars{iVar})),'burst') || ...
            strcmp(char(ncnames(theVars{iVar})),'lat') || ...
            strcmp(char(ncnames(theVars{iVar})),'lon')) ,
        theVars{iVar}.inst_depth = ncfloat(metadata.water_depth-metadata.press.height);
        theVars{iVar}.inst_height = ncfloat(metadata.press.height);
        theVars{iVar}.water_depth = ncfloat(metadata.water_depth);       
    end
end
return


%% -----------------------
function findstar(fid)
while ~feof(fid),
    c = fread(fid,1,'char');
    if c == '*', return; end
end
return

%% ----------------------------
function nc = definetid(ncFile, nrecs, metadata)

nc = netcdf(ncFile, 'clobber');
if isempty(nc), nc = []; return, end

% Dimensions:

nc('time') = nrecs;
nc('depth') = 1;
nc('lat') = 1;
nc('lon') = 1;

% Variables and attributes:
%   coordinate variables:
nc{'time'} = nclong('time');
nc{'time'}.FORTRAN_format = ncchar('F10.2');
nc{'time'}.units = ncchar('True Julian Day');
nc{'time'}.type = ncchar('EVEN');
nc{'time'}.name = ncchar('time');
nc{'time'}.generic_name = ncchar('time');
nc{'time'}.epic_code = nclong(624);
nc{'time'}.comment = ncchar('UT Julian days that begin at midnight; 1968-05-23 = 2440000');
nc{'time'}.comment1 = ncchar('time, taken from burst data, fixed at center of burst');

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

nc{'depth'} = ncfloat('depth') ;
nc{'depth'}.FORTRAN_format = ncchar('F10.2');
nc{'depth'}.units = ncchar('m');
nc{'depth'}.type = ncchar('EVEN');
nc{'depth'}.epic_code = nclong(3);

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

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

%1:P  :PRESSURE (DB)            :depth:dbar:f10.1:
nc{'P_1'} = ncfloat('time', 'depth', 'lat', 'lon'); ncobj = nc{'P_1'}; 
ncobj.name = ncchar('P');
ncobj.long_name = ncchar('PRESSURE (DB)          ');
ncobj.generic_name = ncchar('depth');
ncobj.FORTRAN_format = ncchar('f10.1');
ncobj.units = ncchar('dbar');
ncobj.epic_code = nclong(1);
ncobj.serial_number = ncchar(metadata.press.serial);
ncobj.valid_range = ncfloat([0 10000]);
ncobj.comment = ncchar('Absolute pressure reduced by 10.1325 to remove atmospheric pressure.');

%  20:T  :TEMPERATURE (C)          :temp:C:f10.2:IPTS-1968 standard
nc{'T_20'} = ncfloat('time', 'depth', 'lat', 'lon'); ncobj = nc{'T_20'};
ncobj.name = ncchar('T');
ncobj.long_name = ncchar('TEMPERATURE (C)          ');
ncobj.generic_name = ncchar('temp');
ncobj.FORTRAN_format = ncchar('f10.2');
ncobj.units = ncchar('C');
ncobj.epic_code = nclong(20);
ncobj.serial_number = ncchar(metadata.press.serial);
ncobj.valid_range = ncfloat([0 100]);
ncobj.comment = ncchar('IPTS-1968 Standard');

if isfield(metadata,'cond') && isfield(metadata.cond,'serial') && ...
        ~isempty(metadata.cond.serial),
    
    %    51:C  :CONDUCTIVITY             :con:S/m:f10.3:
    nc{'C_51'} = ncfloat('time', 'depth', 'lat', 'lon'); ncobj = nc{'C_51'};
    ncobj.name = ncchar('C');
    ncobj.long_name = ncchar('CONDUCTIVITY  (S/m)');
    ncobj.generic_name = ncchar('con');
    ncobj.FORTRAN_format = ncchar('f10.2');
    ncobj.units = ncchar('S/m');
    ncobj.epic_code = nclong(51);
    ncobj.serial_number = ncchar(metadata.cond.serial);
    ncobj.sensor_type = ncchar('SBE04');
    ncobj.valid_range = ncfloat([0 100]);

    %  41:S  :SALINITY (PSU)           :sal:PSU: :Practical Salinity Units
    nc{'S_41'} = ncfloat('time', 'depth', 'lat', 'lon'); ncobj = nc{'S_41'};
    ncobj.name = ncchar('S');
    ncobj.long_name = ncchar('SALINITY (PSU)           ');
    ncobj.generic_name = ncchar('sal');
    ncobj.units = ncchar('PSU');
    ncobj.epic_code = nclong(41);
    ncobj.serial_number = ncchar(metadata.cond.serial);
    ncobj.valid_range = ncfloat([0 100]);

end

% basic indeces
nc{'burst'} = ncfloat('time');
nc{'burst'}.units = ncchar('count');
nc{'burst'}.name = ncchar('burst');
nc{'burst'}.generic_name = ncchar('record');
nc{'burst'}.epic_code = nclong(1207);
nc{'burst'}.long_name = ncchar('Burst Number');

return

%% ----------------------------
function nc = definewb(ncFile, metadata)

nc = netcdf(ncFile, 'clobber');

nc('time') = 0; % the record dimension, set fixed as we are in define mode
nc('depth') = 1;
nc('lon') = 1;
nc('lat') = 1;
nc('sample') = metadata.nsamples;
% define first, then load, as we must do statistics, etc.

% basic indeces
nc{'time'} = nclong('time','sample');
nc{'time'}.FORTRAN_format = ncchar('F10.2');
nc{'time'}.units = ncchar('True Julian Day');
nc{'time'}.type = ncchar('EVEN');
nc{'time'}.name = ncchar('time');
nc{'time'}.generic_name = ncchar('time');
nc{'time'}.epic_code = nclong(624);
nc{'time'}.comment = ncchar('UT Julian days that begin at midnight; 1968-05-23 = 2440000');
nc{'time'}.comment1 = ncchar('time, taken from burst data, fixed at center of burst');

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

nc{'depth'} = ncfloat('depth');
nc{'depth'}.FORTRAN_format = ncchar('F10.2');
nc{'depth'}.units = ncchar('m');
nc{'depth'}.name = ncchar('depth');
nc{'depth'}.type = ncchar('EVEN');
nc{'depth'}.generic_name = ncchar('depth');
nc{'depth'}.epic_code = nclong(3);
nc{'depth'}.long_name = ncchar('DEPTH (m)');

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

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

nc{'burst'} = ncfloat('time');
nc{'burst'}.units = ncchar('count');
nc{'burst'}.name = ncchar('burst');
nc{'burst'}.generic_name = ncchar('record');
nc{'burst'}.epic_code = nclong(1207);
nc{'burst'}.long_name = ncchar('Burst Number');
nc{'burst'}.DELTA_T = ncfloat(metadata.dt);

%4022:P   :INTERVAL PRESSURE        :pres:mbar:f10.3:interval pressure measurement, Seadata tripod
nc{'P_4022'} = ncfloat('time', 'sample'); ncobj = nc{'P_4022'};
ncobj.name = ncchar('P');
ncobj.long_name = ncchar('PRESSURE            ');
ncobj.generic_name = ncchar('pres');
ncobj.FORTRAN_format = ncchar('f10.3');
ncobj.epic_code = nclong(4022);
ncobj.units = ncchar('mbar');
ncobj.sensor_type = ncchar(metadata.press.sensor_type);
ncobj.serial_number = ncchar(metadata.press.serial);
ncobj.valid_range = ncfloat([0 1000]);

return

%% --------- Wave Burst Auto-spectrum Statistics ------------------------ %
function nc = definewas(ncFile, metadata)

nc = netcdf(ncFile, 'clobber');

nc.SeaGuageCfg.WaveIntegrationTime = ncfloat(metadata.dt);
nc.SeaGuageCfg.SamplesPerBurst = ncfloat(metadata.nsamples);
nc.SeasoftWavesCfg.MaxNFrequencies = ncfloat(metadata.maxfbands);
nc.SeasoftWavesCfg.FirstFrequency = ncfloat(metadata.f1);
nc.SeasoftWavesCfg.FrequencyInterval = ncfloat(metadata.df);

% now record the settings used in the wave burst processing data setup
nc.SeasoftWavesCfg.AverageWaterColumnTemperature = metadata.SeasoftWaves.temperature; % average water temperature above the pressure sensor, C
nc.SeasoftWavesCfg.AverageWaterColumnSalinity = metadata.SeasoftWaves.salinity; % average salinity above the pressure sensor, PSU
nc.SeasoftWavesCfg.NSpectralEstimates = metadata.SeasoftWaves.estimates; % number of spectral estimates
nc.SeasoftWavesCfg.MinAttenutation = metadata.SeasoftWaves.min_attenuation; % minimum attenuation
nc.SeasoftWavesCfg.MinAutospectrumPeriod = metadata.SeasoftWaves.min_period; % minimum period in autospectrum, s
nc.SeasoftWavesCfg.MaxAutospectrumPeriod = metadata.SeasoftWaves.max_period; % maximum period in autospectrum, s
nc.SeasoftWavesCfg.HanningCutoff = metadata.SeasoftWaves.Hanning_cutoff; % Hannin window cutoff

nc('time') = 0;
nc('depth') = 1;
nc('lat') = 1;
nc('lon') = 1;
nc('frequency') = metadata.maxfbands;

%coordinate variables are the same for both the raw and
%processed NetCDF files
nc{'time'} = nclong('time');
nc{'time'}.FORTRAN_format = ncchar('F10.2');
nc{'time'}.units = ncchar('True Julian Day');
nc{'time'}.type = ncchar('EVEN');
nc{'time'}.name = ncchar('time');
nc{'time'}.generic_name = ncchar('time');
nc{'time'}.epic_code = nclong(624);
nc{'time'}.comment = ncchar('UT Julian days that begin at midnight; 1968-05-23 = 2440000');
nc{'time'}.comment1 = ncchar('time, taken from burst data, fixed at center of burst');

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

nc{'depth'} = ncfloat('depth') ;
nc{'depth'}.FORTRAN_format = ncchar('F10.2');
nc{'depth'}.units = ncchar('m');
nc{'depth'}.type = ncchar('EVEN');
nc{'depth'}.epic_code = nclong(3);

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

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

nc{'burst'} = ncfloat('time');
nc{'burst'}.units = ncchar('count');
nc{'burst'}.name = ncchar('burst');
nc{'burst'}.generic_name = ncchar('record');
nc{'burst'}.epic_code = nclong(1207);
nc{'burst'}.long_name = ncchar('Burst Number');
nc{'burst'}.DELTA_T = ncfloat(metadata.dt);


nc{'wh_4061'} = ncdouble('time','lat','lon');  ncobj = nc{'wh_4061'};
ncobj.long_name=ncchar('Significant Wave Height (m)');
ncobj.name = ncchar('wh');
ncobj.generic_name=ncchar('wave_height');
ncobj.units=ncchar('m');
ncobj.epic_code=nclong(4061);
ncobj.FORTRAN_format=ncchar('F10.2');
ncobj.serial_number = ncchar(metadata.press.serial);
ncobj.valid_range = ncfloat([0 1000]);

%   3:D  :DEPTH (M)                :depth:m:f10.1:
nc{'D_3'} = ncfloat('time','lat','lon'); ncobj = nc{'D_3'};
ncobj.FORTRAN_format = ncchar('F10.1');
ncobj.units = ncchar('m');
ncobj.name = ncchar('D');
ncobj.type = ncchar('EVEN');
ncobj.generic_name = ncchar('depth');
ncobj.epic_code = nclong(3);
ncobj.long_name = ncchar('DEPTH (M)');
ncobj.sensor_type = ncchar(metadata.press.sensor_type);
ncobj.serial_number = ncchar(metadata.press.serial);
ncobj.valid_range = ncfloat([0 10000]);

%  18:hght:height (m)              :height:m:f10.2:
nc{'hght_18'} = ncdouble('time','lat','lon'); ncobj = nc{'hght_18'};
ncobj.long_name=ncchar('Height of the Sea Surface (m)');
ncobj.generic_name=ncchar('height');
ncobj.units=ncchar('m');
ncobj.name = ncchar('hght');
ncobj.epic_code=nclong(18);
ncobj.FORTRAN_format=ncchar('F10.2');
ncobj.NOTE=ncchar('height of sea surface relative to sensor');
ncobj.sensor_type = ncchar(metadata.press.sensor_type);
ncobj.serial_number = ncchar(metadata.press.serial);
ncobj.valid_range = ncfloat([0 10000]);

nc{'frequency'}=ncdouble('frequency','lat','lon'); ncobj = nc{'frequency'};
ncobj.name=ncchar('freq');
ncobj.long_name=ncchar('Frequency (Hz)');
ncobj.generic_name=ncchar('frequency');
ncobj.units=ncchar('Hz');
ncobj.FORTRAN_format=ncchar('F10.2');
ncobj.minimum = ncfloat(0);
ncobj.maximum = ncfloat(0);
nc{'nfreq'}=ncdouble('time','lat','lon'); ncobj = nc{'nfreq'};
ncobj.name=ncchar('nfreq');
ncobj.long_name=ncchar('Number of Frequency Bands Calculated');
ncobj.units=ncchar('count');
ncobj.FORTRAN_format=ncchar('F10.2');
ncobj.minimum = ncfloat(0);
ncobj.maximum = ncfloat(0);

% Seabird calls this the Signficant Period - the frequency band with the
% greatest variance
nc{'wp_sig'} = ncdouble('time','lat','lon');   ncobj = nc{'wp_sig'};
ncobj.long_name=ncchar('Significant Wave Period (s)');
ncobj.generic_name=ncchar('Tsig');
ncobj.units=ncchar('s');
ncobj.FORTRAN_format=ncchar('F10.2');
ncobj.sensor_type = ncchar(metadata.press.sensor_type);
ncobj.serial_number = ncchar(metadata.press.serial);
ncobj.valid_range = ncfloat([0 100000]);

nc{'pspec'}=ncshort('time','frequency','lat','lon'); ncobj = nc{'pspec'};
ncobj.name=ncchar('pspec');
ncobj.long_name=ncchar('Pressure-derived Non-directional Wave Height Spectrum (mm/sqrt(Hz))');
ncobj.generic_name=ncchar('pressure spectrum');
ncobj.units=ncchar('mm/sqrt(Hz)');
ncobj.FORTRAN_format=ncchar('F10.2');
ncobj.sensor_type = ncchar(metadata.press.sensor_type);
ncobj.serial_number = ncchar(metadata.press.serial);
ncobj.valid_range = ncfloat([0 100000]);

%84:RHO:DENSITY (KG/M**3)        :den:kg m-3:f10.2:
nc{'RHO_84'} = ncfloat('time', 'lat', 'lon');  ncobj = nc{'RHO_84'};
ncobj.epic_code = nclong(84);
ncobj.name = ncchar('RHO');
ncobj.long_name = ncchar('DENSITY (KG/M**3)');
ncobj.generic_name = ncchar('den');
ncobj.units = ncchar('kg m-3');
ncobj.FORTRAN_format = ncchar('f10.2');
ncobj.sensor_type = ncchar(metadata.press.sensor_type);
ncobj.serial_number = ncchar(metadata.press.serial);
ncobj.valid_range = ncfloat([0 100000]);

nc{'chi2'} = ncfloat('time', 'lat', 'lon');  ncobj = nc{'chi2'};
ncobj.long_name = ncchar('chi squared confidence interval (%)');
ncobj.units = ncchar('%');
ncobj.valid_range = ncfloat([0 100]);
nc{'chi2upper'} = ncfloat('time', 'lat', 'lon');  ncobj = nc{'chi2upper'};
ncobj.long_name = ncchar('multiplier for chi squared upper bound');
ncobj.units = ncchar('%');
ncobj.valid_range = ncfloat([0 100]);
nc{'chi2lower'} = ncfloat('time', 'lat', 'lon');  ncobj = nc{'chi2lower'};
ncobj.long_name = ncchar('multiplier for chi squared lower bound');
ncobj.units = ncchar('%');
ncobj.valid_range = ncfloat([0 100]);

nc{'tjoul'} = ncfloat('time', 'lat', 'lon');  ncobj = nc{'tjoul'};
ncobj.long_name = ncchar('total energy');
ncobj.units = ncchar('J/m^2');
ncobj.valid_range = ncfloat([0 10000]);
ncobj.serial_number = ncchar(metadata.press.serial);

nc{'tvar'} = ncfloat('time', 'lat', 'lon');  ncobj = nc{'tvar'};
ncobj.long_name = ncchar('total variance');
ncobj.units = ncchar('m^2');
ncobj.valid_range = ncfloat([0 10000]);
ncobj.serial_number = ncchar(metadata.press.serial);

return

%% --------- Surafce Wave Time Series Statistics ------------------------ %
function nc = definewts(ncFile, metadata)

nc = netcdf(ncFile, 'clobber');

nc.SeaGuageCfg.WaveIntegrationTime = ncfloat(metadata.dt);
nc.SeaGuageCfg.SamplesPerBurst = ncfloat(metadata.nsamples);

nc('time') = 0;
nc('depth') = 1;
nc('lat') = 1;
nc('lon') = 1;

%coordinate variables are the same for both the raw and
%processed NetCDF files
nc{'time'} = nclong('time');
nc{'time'}.FORTRAN_format = ncchar('F10.2');
nc{'time'}.units = ncchar('True Julian Day');
nc{'time'}.type = ncchar('EVEN');
nc{'time'}.name = ncchar('time');
nc{'time'}.generic_name = ncchar('time');
nc{'time'}.epic_code = nclong(624);
nc{'time'}.comment = ncchar('UT Julian days that begin at midnight; 1968-05-23 = 2440000');
nc{'time'}.comment1 = ncchar('time, taken from burst data, fixed at center of burst');

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

nc{'depth'} = ncfloat('depth') ;
nc{'depth'}.FORTRAN_format = ncchar('F10.2');
nc{'depth'}.units = ncchar('m');
nc{'depth'}.type = ncchar('EVEN');
nc{'depth'}.epic_code = nclong(3);

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

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

nc{'burst'} = ncfloat('time');
nc{'burst'}.units = ncchar('count');
nc{'burst'}.name = ncchar('burst');
nc{'burst'}.generic_name = ncchar('record');
nc{'burst'}.epic_code = nclong(1207);
nc{'burst'}.long_name = ncchar('Burst Number');
nc{'burst'}.DELTA_T = ncfloat(metadata.dt);

% 4061:wh  :significant wave height   :wave_height:m:f10.2:
nc{'wh_4061'} = ncdouble('time','lat','lon');  ncobj = nc{'wh_4061'};
ncobj.name = ncchar('wh');
ncobj.long_name=ncchar('Significant Wave Height (m)');
ncobj.generic_name=ncchar('Hs');
ncobj.units=ncchar('m');
ncobj.epic_code=nclong(4061);
ncobj.FORTRAN_format=ncchar('F10.2');
ncobj.sensor_type = ncchar(metadata.press.sensor_type);
ncobj.serial_number = ncchar(metadata.press.serial);

%  18:hght:height (m)              :height:m:f10.2:
nc{'hght_18'} = ncdouble('time','lat','lon'); ncobj = nc{'hght_18'};
ncobj.name=ncchar('hght');
ncobj.long_name=ncchar('Height of the Sea Surface (m)');
ncobj.generic_name=ncchar('height');
ncobj.units=ncchar('m');
ncobj.epic_code=nclong(18);
ncobj.FORTRAN_format=ncchar('F10.2');
ncobj.NOTE=ncchar('height of sea surface relative to sensor');
ncobj.sensor_type = ncchar(metadata.press.sensor_type);
ncobj.serial_number = ncchar(metadata.press.serial);

%   3:D  :DEPTH (M)                :depth:m:f10.1:
nc{'D_3'} = ncfloat('time','lat','lon'); ncobj = nc{'D_3'};
ncobj.FORTRAN_format = ncchar('F10.1');
ncobj.units = ncchar('m');
ncobj.name = ncchar('D');
ncobj.type = ncchar('EVEN');
ncobj.generic_name = ncchar('depth');
ncobj.epic_code = nclong(3);
ncobj.long_name = ncchar('DEPTH (M)');
ncobj.serial_number = ncchar(metadata.press.serial);

% 4064:mwh :Maximum Wave Height      :wave_height:m:f10.2:
nc{'mwh_4064'} = ncdouble('time','lat','lon'); ncobj = nc{'mwh_4064'};
ncobj.long_name=ncchar('Maximum Wave Height (m)');
ncobj.name=ncchar('mwh');
ncobj.generic_name=ncchar('wave_height');
ncobj.units=ncchar('m');
ncobj.epic_code=nclong(4064);
ncobj.FORTRAN_format=ncchar('F10.2');
ncobj.minimum = ncfloat(0);
ncobj.maximum = ncfloat(0);
ncobj.sensor_type = ncchar(metadata.press.sensor_type);
ncobj.serial_number = ncchar(metadata.press.serial);

% Seabird calls this the Signficant Period
% 4063:dwp :dominant wave period     :wave_period:s:f10.2:
nc{'dwp_4063'} = ncdouble('time','lat','lon');   ncobj = nc{'dwp_4063'};
ncobj.long_name=ncchar('dominant wave period (s)');
ncobj.name=ncchar('dwp');
ncobj.generic_name=ncchar('wave_period');
ncobj.units=ncchar('s');
ncobj.FORTRAN_format=ncchar('F10.2');
ncobj.minimum = ncfloat(0);
ncobj.maximum = ncfloat(0);
ncobj.sensor_type = ncchar(metadata.press.sensor_type);
ncobj.serial_number = ncchar(metadata.press.serial);

% 4060:wp  :average wave period       :wave_period:s:f10.2:
nc{'wp_4060'} = ncdouble('time','lat','lon');   ncobj = nc{'wp_4060'};
ncobj.name=ncchar('wp');
ncobj.long_name=ncchar('average wave period (s)');
ncobj.generic_name=ncchar('wave_period');
ncobj.units=ncchar('s');
ncobj.FORTRAN_format=ncchar('F10.2');
ncobj.minimum = ncfloat(0);
ncobj.maximum = ncfloat(0);
ncobj.sensor_type = ncchar(metadata.press.sensor_type);
ncobj.serial_number = ncchar(metadata.press.serial);

nc{'awh'} = ncdouble('time','lat','lon'); ncobj = nc{'awh'};
ncobj.long_name=ncchar('Average Wave Height (m)');
ncobj.name=ncchar('awh');
ncobj.generic_name=ncchar('wave_height');
ncobj.units=ncchar('m');
ncobj.epic_code=nclong(0);
ncobj.FORTRAN_format=ncchar('F10.2');
ncobj.sensor_type = ncchar(metadata.press.sensor_type);
ncobj.serial_number = ncchar(metadata.press.serial);

nc{'awh10'} = ncdouble('time','lat','lon'); ncobj = nc{'awh100'};
ncobj.long_name=ncchar('Average Height of largest 1/10 waves (m)');
ncobj.name=ncchar('awh10');
ncobj.generic_name=ncchar('wave_height');
ncobj.units=ncchar('m');
ncobj.epic_code=nclong(0);
ncobj.FORTRAN_format=ncchar('F10.2');
ncobj.sensor_type = ncchar(metadata.press.sensor_type);
ncobj.serial_number = ncchar(metadata.press.serial);

nc{'awh100'} = ncdouble('time','lat','lon'); ncobj = nc{'awh100'};
ncobj.long_name=ncchar('Average Wave Height (m)');
ncobj.name=ncchar('awh100');
ncobj.generic_name=ncchar('wave_height');
ncobj.units=ncchar('m');
ncobj.epic_code=nclong(0);
ncobj.FORTRAN_format=ncchar('F10.2');
ncobj.sensor_type = ncchar(metadata.press.sensor_type);
ncobj.serial_number = ncchar(metadata.press.serial);

nc{'nw'} = ncdouble('time','lat','lon'); ncobj = nc{'nw'};
ncobj.long_name=ncchar('Number of Waves');
ncobj.name=ncchar('nw');
ncobj.generic_name=ncchar('wave_count');
ncobj.units=ncchar('count');
ncobj.epic_code=nclong(0);
ncobj.FORTRAN_format=ncchar('F10.2');
ncobj.sensor_type = ncchar(metadata.press.sensor_type);
ncobj.serial_number = ncchar(metadata.press.serial);

return
%% How Seabird does it...
% Process Wave Burst Data Algorithm
%
% For each wave burst, Process Wave Burst Data performs an auto spectrum analysis:
%
% ·	Reads burst into an array.
%
% ·	Removes mean, saves mean value.
%
% ·	Uses mean value, and average water temperature and average salinity (input by user) to compute density.
%
% ·	Removes trend.
%
% ·	Makes array a power of two.
%
% ·	Applies Hanning window to suppress side-lobe leakage.
%
% ·	Adjusts scale factor to account for tapering by Hanning window.
%
% ·	Fast Fourier Transforms to create raw spectral estimates of subsurface pressure.
%
% ·	Computes maximum frequency to process.
%
% ·	Sets Fourier coefficients greater than maximum frequency or less than minimum frequency to 0; or
% (if Use filter cutoff is selected) applies a filter that ramps the Fourier coefficients down to 0 for frequencies greater than maximum frequency or less than minimum frequency.
%
% ·	Saves Fourier coefficients.
%
% ·	Band averages raw spectral estimates to create auto-spectrum.
%
% ·	Applies dispersion transfer function to band center frequencies.
%
% ·	Calculates wave statistics from auto-spectrum: variance, energy, significant wave height, and significant period.
%
% Using the saved (non-band-averaged) Fourier coefficients, Process Wave Burst Data performs a surface wave zero crossing analysis:
%
% ·	Applies dispersion transfer function to each frequency.
%
% ·	Inverse Fast Fourier Transforms to create surface wave time series.
%
% ·	Applies inverse Hanning window and adjusts scale factor.
%
% ·	Zeroes all elements where inverse Hanning factor is greater than 10.
%
% ·	Performs zero crossing analysis of surface wave time series to create an array of individual waves and their corresponding periods.
%
% ·	Sorts wave array in ascending order.
%
% ·	Calculates wave statistics from surface wave time series:
% average wave height, average period, maximum wave height, significant period, significant wave height H1/3 , H1/10 , H1/100.
%