function create_cbathy_ncfile(bathy,xs)

% Matlab script generated by netCDF kickstarter
% at 2014-02-10T11:08Z
% http://publicwiki.deltares.nl/display/OET/netCDF%20kickstarter

oetsettings;

%% GET PARAMS
station = regexp(bathy(1).sName,'(?<=\.[0-9]{4}\.)[A-Za-z0-9]*(?=\.cx\.)','match');
station = station{1};
time = char(bathy(:).sName);
time = str2num(time(:,1:10));

begstr = datestr(str2num(bathy(1).sName(1:10))/(24*3600)+datenum(1970,1,1),'yyyy_mm_dd');
endstr = datestr(str2num(bathy(end).sName(1:10))/(24*3600)+datenum(1970,1,1),'yyyy_mm_dd');

[xRD,yRD] = argus2rd(station,bathy(1).xm,bathy(1).ym);
[lon,lat] = convertCoordinates(xRD,yRD,'CS1.code',28992,'CS2.code',4326);

%% CREATE EMPTY NETCDF FILE
ncfile = ['cbathy_' station '_from_' begstr '_till_' endstr '.nc'];
nc_create_empty(ncfile);

%% ADD DIMENSIONS
nc_adddim(ncfile, 'x', length(bathy(1).xm));
nc_adddim(ncfile, 'y', length(bathy(1).ym));
nc_adddim(ncfile, 'time', length(bathy));
nc_adddim(ncfile, 'fB', bathy(1).params.nKeep);
nc_adddim(ncfile, 'nv', 2);
nc_adddim(ncfile, 'nv2', 4);

%% ADD GLOBAL ATTRIBUTES
% see http://www.unidata.ucar.edu/software/thredds/current/netcdf-java/formats/DataDiscoveryAttConvention.html
nc_attput(ncfile, nc_global, 'Conventions', 'CF-1.6');
nc_attput(ncfile, nc_global, 'Metadata_Conventions', 'Unidata Dataset Discovery v1.0');
nc_attput(ncfile, nc_global, 'featureType', 'grid');
nc_attput(ncfile, nc_global, 'cdm_data_type', 'grid');
nc_attput(ncfile, nc_global, 'standard_name_vocabulary', 'CF-1.6');
nc_attput(ncfile, nc_global, 'title', 'cBathy output');
nc_attput(ncfile, nc_global, 'summary', 'Output of the cBathy algorithm for Argus imagery');
nc_attput(ncfile, nc_global, 'source', ['Argus station ' station]);
nc_attput(ncfile, nc_global, 'id', ['cBathy_' station '_' ]);
nc_attput(ncfile, nc_global, 'naming_authority', 'tudelft.nl');
nc_attput(ncfile, nc_global, 'time_coverage_start', datestr(bathy(1).epoch/24/3600 + datenum(1970,1,1),'yyyy-mm-ddTHH:MM:SSZ'));
nc_attput(ncfile, nc_global, 'time_coverage_end', datestr(bathy(end).epoch/24/3600 + datenum(1970,1,1),'yyyy-mm-ddTHH:MM:SSZ'));
nc_attput(ncfile, nc_global, 'time_coverage_resolution', [int2str(round(min(diff([bathy(:).epoch]))/3600)) ' hours']);
nc_attput(ncfile, nc_global, 'geospatial_lat_min', 0);
nc_attput(ncfile, nc_global, 'geospatial_lat_max', 0);
nc_attput(ncfile, nc_global, 'geospatial_lat_units', 'degrees_north');
nc_attput(ncfile, nc_global, 'geospatial_lat_resolution', '');
nc_attput(ncfile, nc_global, 'geospatial_lon_min', 0);
nc_attput(ncfile, nc_global, 'geospatial_lon_max', 0);
nc_attput(ncfile, nc_global, 'geospatial_lon_units', 'degrees_east');
nc_attput(ncfile, nc_global, 'geospatial_lon_resolution', '');
nc_attput(ncfile, nc_global, 'geospatial_vertical_min', 0);
nc_attput(ncfile, nc_global, 'geospatial_vertical_max', 0);
nc_attput(ncfile, nc_global, 'geospatial_vertical_units', '');
nc_attput(ncfile, nc_global, 'geospatial_vertical_resolution', '');
nc_attput(ncfile, nc_global, 'geospatial_vertical_positive', '');
nc_attput(ncfile, nc_global, 'institution', 'Delft University of Technology');
nc_attput(ncfile, nc_global, 'creator_name', 'Max Radermacher');
nc_attput(ncfile, nc_global, 'creator_url', 'http://www.naturecoast.nl');
nc_attput(ncfile, nc_global, 'creator_email', 'm.radermacher@tudelft.nl');
nc_attput(ncfile, nc_global, 'project', 'STW NatureCoast');
nc_attput(ncfile, nc_global, 'processing_level', 'preliminary');
nc_attput(ncfile, nc_global, 'references', '');
nc_attput(ncfile, nc_global, 'keywords_vocabulary', 'NASA/GCMD Earth Science Keywords. Version 6.0');
nc_attput(ncfile, nc_global, 'keywords', 'Argus, Depth inversion, Coastal imagery, Bathymetry');

if strcmpi(station,'kijkduin')
    ack = 'Delft University of Technology and Coastal Imaging Lab at Oregon State University';
    lic = ['These data can be used freely for research purposes provided that the following source is acknowledged: ',...
        'TU Delft. disclaimer: This data is made available in the hope that it will be useful, but ',...
        'WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.'];
elseif strcmpi(station,'zandmotor')
    ack = 'Rijkswaterstaat, STW NatureCoast and Coastal Imaging Lab at Oregon State University';
    lic = ['These data can be used freely for research purposes as from January 1st 2015 and within the Zandmotor consortium before that date, ',...
        'provided that the following sources are acknowledged: RIJKSWATERSTAAT and STW NatureCoast. disclaimer: This data is made available ',...
        'in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.'];
end

nc_attput(ncfile, nc_global, 'acknowledgment', ack);
nc_attput(ncfile, nc_global, 'comment', '');
nc_attput(ncfile, nc_global, 'contributor_name', 'Meagan Wengrove');
nc_attput(ncfile, nc_global, 'contributor_role', 'Initiator of cBathy processing at Kijkduin and Zandmotor stations');
nc_attput(ncfile, nc_global, 'date_created', '2014-02-10T11:08Z');
nc_attput(ncfile, nc_global, 'date_modified', '2014-02-10T11:08Z');
nc_attput(ncfile, nc_global, 'date_issued', '2014-02-10T11:08Z');
nc_attput(ncfile, nc_global, 'publisher_name', 'Max Radermacher');
nc_attput(ncfile, nc_global, 'publisher_email', 'm.radermacher@tudelft.nl');
nc_attput(ncfile, nc_global, 'publisher_url', 'http://www.naturecoast.nl');
nc_attput(ncfile, nc_global, 'history', '');
nc_attput(ncfile, nc_global, 'license', lic);
nc_attput(ncfile, nc_global, 'metadata_link', '0');
nc_attput(ncfile, nc_global, 'options',var2evalstr(bathy(1).params, 'basevarname', 'params'));

%% ADD VARIABLES
var = struct(...
    'Name', 'x',...
    'Datatype', 'double',...
    'Dimension', {{'x'}}); % note the double brackets
nc_addvar(ncfile, var);
nc_attput(ncfile, 'x', 'long_name', 'x-coordinate');
nc_attput(ncfile, 'x', 'standard_name', 'projection_x_coordinate');
nc_attput(ncfile, 'x', 'units', 'm');
nc_attput(ncfile, 'x', 'axis', 'X');
nc_attput(ncfile, 'x', 'valid_min', 0);
nc_attput(ncfile, 'x', 'valid_max', 0);
nc_attput(ncfile, 'x', 'bounds', 'x_bounds');
% nc_attput(ncfile, 'x', 'grid_mapping', 'crs');
nc_attput(ncfile, 'x', 'comment', 'Local Argus coordinate system. Origin: Orientation:');
nc_varput(ncfile, 'x', bathy(1).xm);

var = struct(...
    'Name', 'y',...
    'Datatype', 'double',...
    'Dimension', {{'y'}}); % note the double brackets
nc_addvar(ncfile, var);
nc_attput(ncfile, 'y', 'long_name', 'y-coordinate');
nc_attput(ncfile, 'y', 'standard_name', 'projection_y_coordinate');
nc_attput(ncfile, 'y', 'units', 'm');
nc_attput(ncfile, 'y', 'axis', 'Y');
nc_attput(ncfile, 'y', 'valid_min', 0);
nc_attput(ncfile, 'y', 'valid_max', 0);
nc_attput(ncfile, 'y', 'bounds', 'y_bounds');
% nc_attput(ncfile, 'y', 'grid_mapping', 'crs');
nc_attput(ncfile, 'y', 'comment', 'Local Argus coordinate system. Origin: Orientation:');
nc_varput(ncfile, 'y', bathy(1).ym);

var = struct(...
    'Name', 'xRD',...
    'Datatype', 'double',...
    'Dimension', {{'y','x'}}); % note the double brackets
nc_addvar(ncfile, var);
nc_attput(ncfile, 'xRD', 'long_name', 'x-coordinate RD system');
nc_attput(ncfile, 'xRD', 'standard_name', 'projection_x_coordinate');
nc_attput(ncfile, 'xRD', 'units', 'm');
nc_attput(ncfile, 'xRD', 'valid_min', 0);
nc_attput(ncfile, 'xRD', 'valid_max', 0);
nc_attput(ncfile, 'xRD', 'grid_mapping', 'crs');
nc_attput(ncfile, 'xRD', 'comment', '');
nc_varput(ncfile, 'xRD', xRD);

var = struct(...
    'Name', 'yRD',...
    'Datatype', 'double',...
    'Dimension', {{'y','x'}}); % note the double brackets
nc_addvar(ncfile, var);
nc_attput(ncfile, 'yRD', 'long_name', 'y-coordinate RD system');
nc_attput(ncfile, 'yRD', 'standard_name', 'projection_y_coordinate');
nc_attput(ncfile, 'yRD', 'units', 'm');
nc_attput(ncfile, 'yRD', 'valid_min', 0);
nc_attput(ncfile, 'yRD', 'valid_max', 0);
nc_attput(ncfile, 'yRD', 'grid_mapping', 'crs');
nc_attput(ncfile, 'yRD', 'comment', '');
nc_varput(ncfile, 'yRD', yRD);

var = struct(...
    'Name', 'lat',...
    'Datatype', 'double',...
    'Dimension', {{'y','x'}}); % note the double brackets
nc_addvar(ncfile, var);
nc_attput(ncfile, 'lat', 'long_name', 'latitude');
nc_attput(ncfile, 'lat', 'standard_name', 'latitude');
nc_attput(ncfile, 'lat', 'units', 'degrees_north');
nc_attput(ncfile, 'lat', 'valid_min', 0);
nc_attput(ncfile, 'lat', 'valid_max', 0);
nc_attput(ncfile, 'lat', 'bounds', 'lat_bounds');
nc_attput(ncfile, 'lat', 'ancillary_variables', '');
nc_attput(ncfile, 'lat', 'comment', '');
nc_varput(ncfile, 'lat', lat);

var = struct(...
    'Name', 'lon',...
    'Datatype', 'double',...
    'Dimension', {{'y','x'}}); % note the double brackets
nc_addvar(ncfile, var);
nc_attput(ncfile, 'lon', 'long_name', 'longitude');
nc_attput(ncfile, 'lon', 'standard_name', 'longitude');
nc_attput(ncfile, 'lon', 'units', 'degrees_east');
nc_attput(ncfile, 'lon', 'valid_min', 0);
nc_attput(ncfile, 'lon', 'valid_max', 0);
nc_attput(ncfile, 'lon', 'bounds', 'lon_bounds');
nc_attput(ncfile, 'lon', 'ancillary_variables', '');
nc_attput(ncfile, 'lon', 'comment', '');
nc_varput(ncfile, 'lon', lon);

var = struct(...
    'Name', 'time',...
    'Datatype', 'double',...
    'Dimension', {{'time'}}); % note the double brackets
nc_addvar(ncfile, var);
nc_attput(ncfile, 'time', 'long_name', 'time');
nc_attput(ncfile, 'time', 'standard_name', 'time');
nc_attput(ncfile, 'time', 'units', 'seconds since 1970-01-01 00:00:00 0:00');
nc_attput(ncfile, 'time', 'calendar', 'julian');
nc_attput(ncfile, 'time', 'axis', 'T');
nc_attput(ncfile, 'time', 'bounds', 'time_bounds');
nc_attput(ncfile, 'time', 'ancillary_variables', '');
nc_attput(ncfile, 'time', 'comment', '');
nc_varput(ncfile, 'time', time);

var = struct(...
    'Name', 'x_bounds',...
    'Datatype', 'double',...
    'Dimension', {{'x','nv'}}); % note the double brackets
nc_addvar(ncfile, var);
nc_attput(ncfile, 'x_bounds', 'units', 'm');
nc_attput(ncfile, 'x_bounds', 'comment', 'x-coordinate values at the upper and lower bounds of each pixel.');

var = struct(...
    'Name', 'y_bounds',...
    'Datatype', 'double',...
    'Dimension', {{'y','nv'}}); % note the double brackets
nc_addvar(ncfile, var);
nc_attput(ncfile, 'y_bounds', 'units', 'm');
nc_attput(ncfile, 'y_bounds', 'comment', 'y-coordinate values at the left and right bounds of each pixel.');

var = struct(...
    'Name', 'lat_bounds',...
    'Datatype', 'double',...
    'Dimension', {{'y','x','nv2'}}); % note the double brackets
nc_addvar(ncfile, var);
nc_attput(ncfile, 'lat_bounds', 'units', 'degrees_north');
nc_attput(ncfile, 'lat_bounds', 'comment', 'latitude values at the north and south bounds of each pixel.');

var = struct(...
    'Name', 'lon_bounds',...
    'Datatype', 'double',...
    'Dimension', {{'y','x','nv2'}}); % note the double brackets
nc_addvar(ncfile, var);
nc_attput(ncfile, 'lon_bounds', 'units', 'degrees_east');
nc_attput(ncfile, 'lon_bounds', 'comment', 'longitude values at the west and east bounds of each pixel.');

var = struct(...
    'Name', 'time_bounds',...
    'Datatype', 'double',...
    'Dimension', {{'time','nv'}}); % note the double brackets
nc_addvar(ncfile, var);
nc_attput(ncfile, 'time_bounds', 'units', 'seconds since 1970-01-01 00:00:00 0:00');
nc_attput(ncfile, 'time_bounds', 'comment', 'time bounds for each time value');

var = struct(...
    'Name', 'x0',...
    'Datatype', 'double',...
    'Dimension', {{'y'}}); % note the double brackets
nc_addvar(ncfile, var);
nc_attput(ncfile, 'x0', 'long_name', 'x coordinate of breaker bar');
nc_attput(ncfile, 'x0', 'standard_name', '');
nc_attput(ncfile, 'x0', 'units', 'm');
nc_attput(ncfile, 'x0', 'scale_factor', 1.0);
nc_attput(ncfile, 'x0', 'add_offset', 0.0);
nc_attput(ncfile, 'x0', 'valid_min', 0);
nc_attput(ncfile, 'x0', 'valid_max', 0);
nc_attput(ncfile, 'x0', 'coordinates', 'y');
%nc_attput(ncfile, 'x0', 'grid_mapping', 'crs');
nc_attput(ncfile, 'x0', 'source', '');
nc_attput(ncfile, 'x0', 'references', '');
nc_attput(ncfile, 'x0', 'cell_methods', '');
nc_attput(ncfile, 'x0', 'ancillary_variables', '');
nc_attput(ncfile, 'x0', 'comment', '');
nc_varput(ncfile, 'x0', xs);

var = struct(...
    'Name', 'depth',...
    'Datatype', 'double',...
    'Dimension', {{'time','y','x'}}); % note the double brackets
nc_addvar(ncfile, var);
nc_attput(ncfile, 'depth', 'long_name', 'water depth below NAP');
nc_attput(ncfile, 'depth', 'standard_name', 'sea_floor_depth_below_sea_level');
nc_attput(ncfile, 'depth', 'units', 'm');
nc_attput(ncfile, 'depth', 'scale_factor', 1.0);
nc_attput(ncfile, 'depth', 'add_offset', 0.0);
nc_attput(ncfile, 'depth', 'valid_min', 0);
nc_attput(ncfile, 'depth', 'valid_max', 0);
nc_attput(ncfile, 'depth', 'coordinates', 'time y x');
%nc_attput(ncfile, 'depth', 'grid_mapping', 'crs');
nc_attput(ncfile, 'depth', 'source', '');
nc_attput(ncfile, 'depth', 'references', '');
nc_attput(ncfile, 'depth', 'cell_methods', '');
nc_attput(ncfile, 'depth', 'ancillary_variables', '');
nc_attput(ncfile, 'depth', 'comment', '');

h = nan(length(bathy),length(bathy(1).ym),length(bathy(1).xm));
herr = nan(length(bathy),length(bathy(1).ym),length(bathy(1).xm));
hkal = nan(length(bathy),length(bathy(1).ym),length(bathy(1).xm));
hkalerr = nan(length(bathy),length(bathy(1).ym),length(bathy(1).xm));
K = nan(length(bathy),length(bathy(1).ym),length(bathy(1).xm));
pr = cell(length(bathy),1);
zt = nan(length(bathy),1);
Hs = nan(length(bathy),1);
for i = 1:length(bathy)
    h(i,:,:) = bathy(i).fCombined.h;
    herr(i,:,:) = bathy(i).fCombined.hErr;
    hkal(i,:,:) = bathy(i).runningAverage.h;
    hkalerr(i,:,:) = bathy(i).runningAverage.hErr;
    K(i,:,:) = bathy(i).runningAverage.K;
    pr{i} = bathy(i).runningAverage.prior;
    zt(i) = bathy(i).tide.zt;
    Hs(i) = bathy(i).runningAverage.Hs;
end
nc_varput(ncfile, 'depth', h);

var = struct(...
    'Name', 'kalman_depth',...
    'Datatype', 'double',...
    'Dimension', {{'time','y','x'}}); % note the double brackets
nc_addvar(ncfile, var);
nc_attput(ncfile, 'kalman_depth', 'long_name', 'kalman filtered water depth below NAP');
nc_attput(ncfile, 'kalman_depth', 'standard_name', 'sea_floor_depth_below_sea_level');
nc_attput(ncfile, 'kalman_depth', 'units', 'm');
nc_attput(ncfile, 'kalman_depth', 'scale_factor', 1.0);
nc_attput(ncfile, 'kalman_depth', 'add_offset', 0.0);
nc_attput(ncfile, 'kalman_depth', 'valid_min', 0);
nc_attput(ncfile, 'kalman_depth', 'valid_max', 0);
nc_attput(ncfile, 'kalman_depth', 'coordinates', 'time y x');
%nc_attput(ncfile, 'kalman_depth', 'grid_mapping', 'crs');
nc_attput(ncfile, 'kalman_depth', 'source', '');
nc_attput(ncfile, 'kalman_depth', 'references', '');
nc_attput(ncfile, 'kalman_depth', 'cell_methods', '');
nc_attput(ncfile, 'kalman_depth', 'ancillary_variables', '');
nc_attput(ncfile, 'kalman_depth', 'comment', '');
nc_varput(ncfile, 'kalman_depth', hkal);

var = struct(...
    'Name', 'depth_error',...
    'Datatype', 'double',...
    'Dimension', {{'time','y','x'}}); % note the double brackets
nc_addvar(ncfile, var);
nc_attput(ncfile, 'depth_error', 'long_name', 'error estimate water depth');
nc_attput(ncfile, 'depth_error', 'standard_name', '');
nc_attput(ncfile, 'depth_error', 'units', 'm');
nc_attput(ncfile, 'depth_error', 'scale_factor', 1.0);
nc_attput(ncfile, 'depth_error', 'add_offset', 0.0);
nc_attput(ncfile, 'depth_error', 'valid_min', 0);
nc_attput(ncfile, 'depth_error', 'valid_max', 0);
nc_attput(ncfile, 'depth_error', 'coordinates', 'time y x');
%nc_attput(ncfile, 'depth_error', 'grid_mapping', 'crs');
nc_attput(ncfile, 'depth_error', 'source', '');
nc_attput(ncfile, 'depth_error', 'references', '');
nc_attput(ncfile, 'depth_error', 'cell_methods', '');
nc_attput(ncfile, 'depth_error', 'ancillary_variables', '');
nc_attput(ncfile, 'depth_error', 'comment', '');
nc_varput(ncfile, 'depth_error', herr);

var = struct(...
    'Name', 'kalman_error',...
    'Datatype', 'double',...
    'Dimension', {{'time','y','x'}}); % note the double brackets
nc_addvar(ncfile, var);
nc_attput(ncfile, 'kalman_error', 'long_name', 'error estimate kalman filtered water depth');
nc_attput(ncfile, 'kalman_error', 'standard_name', '');
nc_attput(ncfile, 'kalman_error', 'units', 'm');
nc_attput(ncfile, 'kalman_error', 'scale_factor', 1.0);
nc_attput(ncfile, 'kalman_error', 'add_offset', 0.0);
nc_attput(ncfile, 'kalman_error', 'valid_min', 0);
nc_attput(ncfile, 'kalman_error', 'valid_max', 0);
nc_attput(ncfile, 'kalman_error', 'coordinates', 'time y x');
%nc_attput(ncfile, 'kalman_error', 'grid_mapping', 'crs');
nc_attput(ncfile, 'kalman_error', 'source', '');
nc_attput(ncfile, 'kalman_error', 'references', '');
nc_attput(ncfile, 'kalman_error', 'cell_methods', '');
nc_attput(ncfile, 'kalman_error', 'ancillary_variables', '');
nc_attput(ncfile, 'kalman_error', 'comment', '');
nc_varput(ncfile, 'kalman_error', hkalerr);

var = struct(...
    'Name', 'tidal_elevation',...
    'Datatype', 'double',...
    'Dimension', {{'time'}}); % note the double brackets
nc_addvar(ncfile, var);
nc_attput(ncfile, 'tidal_elevation', 'long_name', 'tidal elevation linearly interpolated from Hoek van Holland and Scheveningen');
nc_attput(ncfile, 'tidal_elevation', 'standard_name', 'sea_surface_height_above_sea_level');
nc_attput(ncfile, 'tidal_elevation', 'units', 'm');
nc_attput(ncfile, 'tidal_elevation', 'scale_factor', 1.0);
nc_attput(ncfile, 'tidal_elevation', 'add_offset', 0.0);
nc_attput(ncfile, 'tidal_elevation', 'valid_min', 0);
nc_attput(ncfile, 'tidal_elevation', 'valid_max', 0);
nc_attput(ncfile, 'tidal_elevation', 'coordinates', 'time y x');
%nc_attput(ncfile, 'tidal_elevation', 'grid_mapping', 'crs');
nc_attput(ncfile, 'tidal_elevation', 'source', '');
nc_attput(ncfile, 'tidal_elevation', 'references', '');
nc_attput(ncfile, 'tidal_elevation', 'cell_methods', '');
nc_attput(ncfile, 'tidal_elevation', 'ancillary_variables', '');
nc_attput(ncfile, 'tidal_elevation', 'comment', '');
nc_varput(ncfile, 'tidal_elevation', zt);

var = struct(...
    'Name', 'fB',...
    'Datatype', 'double',...
    'Dimension', {{'time','y','x','fB'}}); % note the double brackets
nc_addvar(ncfile, var);
nc_attput(ncfile, 'fB', 'long_name', 'coherence ranked frequency bands');
nc_attput(ncfile, 'fB', 'standard_name', 'wave_frequency');
nc_attput(ncfile, 'fB', 'units', 's-1');
nc_attput(ncfile, 'fB', 'scale_factor', 1.0);
nc_attput(ncfile, 'fB', 'add_offset', 0.0);
nc_attput(ncfile, 'fB', 'valid_min', 0);
nc_attput(ncfile, 'fB', 'valid_max', 0);
nc_attput(ncfile, 'fB', 'coordinates', 'time y x');
%nc_attput(ncfile, 'fB', 'grid_mapping', 'crs');
nc_attput(ncfile, 'fB', 'source', '');
nc_attput(ncfile, 'fB', 'references', '');
nc_attput(ncfile, 'fB', 'cell_methods', '');
nc_attput(ncfile, 'fB', 'ancillary_variables', '');
nc_attput(ncfile, 'fB', 'comment', '');
fB = nan(length(bathy),length(bathy(1).ym),length(bathy(1).xm),bathy(1).params.nKeep);
k = nan(length(bathy),length(bathy(1).ym),length(bathy(1).xm),bathy(1).params.nKeep);
alpha = nan(length(bathy),length(bathy(1).ym),length(bathy(1).xm),bathy(1).params.nKeep);
for i = 1:length(bathy)
    fB(i,:,:,:) = bathy(i).fDependent.fB;
    k(i,:,:,:) = bathy(i).fDependent.k;
    alpha(i,:,:,:) = bathy(i).fDependent.a;
end
nc_varput(ncfile, 'fB', fB);

var = struct(...
    'Name', 'k',...
    'Datatype', 'double',...
    'Dimension', {{'time','y','x','fB'}}); % note the double brackets
nc_addvar(ncfile, var);
nc_attput(ncfile, 'k', 'long_name', 'wavenumbers found for coherence ranked frequency bands');
nc_attput(ncfile, 'k', 'standard_name', '');
nc_attput(ncfile, 'k', 'units', 'm-1');
nc_attput(ncfile, 'k', 'scale_factor', 1.0);
nc_attput(ncfile, 'k', 'add_offset', 0.0);
nc_attput(ncfile, 'k', 'valid_min', 0);
nc_attput(ncfile, 'k', 'valid_max', 0);
nc_attput(ncfile, 'k', 'coordinates', 'time y x');
%nc_attput(ncfile, 'k', 'grid_mapping', 'crs');
nc_attput(ncfile, 'k', 'source', '');
nc_attput(ncfile, 'k', 'references', '');
nc_attput(ncfile, 'k', 'cell_methods', '');
nc_attput(ncfile, 'k', 'ancillary_variables', '');
nc_attput(ncfile, 'k', 'comment', '');
nc_varput(ncfile, 'k', k);

var = struct(...
    'Name', 'alpha',...
    'Datatype', 'double',...
    'Dimension', {{'time','y','x','fB'}}); % note the double brackets
nc_addvar(ncfile, var);
nc_attput(ncfile, 'alpha', 'long_name', 'wave angle found for coherence ranked frequency bands with respect to local Argus x-coordinate axis');
nc_attput(ncfile, 'alpha', 'standard_name', '');
nc_attput(ncfile, 'alpha', 'units', 'radians');
nc_attput(ncfile, 'alpha', 'scale_factor', 1.0);
nc_attput(ncfile, 'alpha', 'add_offset', 0.0);
nc_attput(ncfile, 'alpha', 'valid_min', 0);
nc_attput(ncfile, 'alpha', 'valid_max', 0);
nc_attput(ncfile, 'alpha', 'coordinates', 'time y x');
%nc_attput(ncfile, 'alpha', 'grid_mapping', 'crs');
nc_attput(ncfile, 'alpha', 'source', '');
nc_attput(ncfile, 'alpha', 'references', '');
nc_attput(ncfile, 'alpha', 'cell_methods', '');
nc_attput(ncfile, 'alpha', 'ancillary_variables', '');
nc_attput(ncfile, 'alpha', 'comment', '');
nc_varput(ncfile, 'alpha', alpha);

var = struct(...
    'Name', 'kalman_gain',...
    'Datatype', 'double',...
    'Dimension', {{'time','y','x'}}); % note the double brackets
nc_addvar(ncfile, var);
nc_attput(ncfile, 'kalman_gain', 'long_name', 'kalman gain');
nc_attput(ncfile, 'kalman_gain', 'standard_name', '');
nc_attput(ncfile, 'kalman_gain', 'units', '1');
nc_attput(ncfile, 'kalman_gain', 'scale_factor', 1.0);
nc_attput(ncfile, 'kalman_gain', 'add_offset', 0.0);
nc_attput(ncfile, 'kalman_gain', 'valid_min', 0);
nc_attput(ncfile, 'kalman_gain', 'valid_max', 0);
nc_attput(ncfile, 'kalman_gain', 'coordinates', 'time y x');
%nc_attput(ncfile, 'kalman_gain', 'grid_mapping', 'crs');
nc_attput(ncfile, 'kalman_gain', 'source', '');
nc_attput(ncfile, 'kalman_gain', 'references', '');
nc_attput(ncfile, 'kalman_gain', 'cell_methods', '');
nc_attput(ncfile, 'kalman_gain', 'ancillary_variables', '');
nc_attput(ncfile, 'kalman_gain', 'comment', '');
nc_varput(ncfile, 'kalman_gain', K);

var = struct(...
    'Name', 'kalman_prior_stack',...
    'Datatype', 'char',...
    'Dimension', {{'time'}}); % note the double brackets
nc_addvar(ncfile, var);
nc_attput(ncfile, 'kalman_prior_stack', 'long_name', 'file name of pixel stack prior to present stack in kalman filtering procedure');
nc_attput(ncfile, 'kalman_prior_stack', 'standard_name', '');
nc_attput(ncfile, 'kalman_prior_stack', 'units', '');
nc_attput(ncfile, 'kalman_prior_stack', 'scale_factor', 1.0);
nc_attput(ncfile, 'kalman_prior_stack', 'add_offset', 0.0);
nc_attput(ncfile, 'kalman_prior_stack', 'valid_min', 0);
nc_attput(ncfile, 'kalman_prior_stack', 'valid_max', 0);
nc_attput(ncfile, 'kalman_prior_stack', 'coordinates', 'time y x');
%nc_attput(ncfile, 'kalman_prior_stack', 'grid_mapping', 'crs');
nc_attput(ncfile, 'kalman_prior_stack', 'source', '');
nc_attput(ncfile, 'kalman_prior_stack', 'references', '');
nc_attput(ncfile, 'kalman_prior_stack', 'cell_methods', '');
nc_attput(ncfile, 'kalman_prior_stack', 'ancillary_variables', '');
nc_attput(ncfile, 'kalman_prior_stack', 'comment', '');
nc_varput(ncfile, 'kalman_prior_stack', pr);

var = struct(...
    'Name', 'wave_height',...
    'Datatype', 'double',...
    'Dimension', {{'time'}}); % note the double brackets
nc_addvar(ncfile, var);
nc_attput(ncfile, 'wave_height', 'long_name', 'wave height interpolated from WaveWatch3 for process error computation');
nc_attput(ncfile, 'wave_height', 'standard_name', 'sea_surface_wave_significant_height');
nc_attput(ncfile, 'wave_height', 'units', 'm');
nc_attput(ncfile, 'wave_height', 'scale_factor', 1.0);
nc_attput(ncfile, 'wave_height', 'add_offset', 0.0);
nc_attput(ncfile, 'wave_height', 'valid_min', 0);
nc_attput(ncfile, 'wave_height', 'valid_max', 0);
nc_attput(ncfile, 'wave_height', 'coordinates', 'time y x');
%nc_attput(ncfile, 'wave_height', 'grid_mapping', 'crs');
nc_attput(ncfile, 'wave_height', 'source', '');
nc_attput(ncfile, 'wave_height', 'references', '');
nc_attput(ncfile, 'wave_height', 'cell_methods', '');
nc_attput(ncfile, 'wave_height', 'ancillary_variables', '');
nc_attput(ncfile, 'wave_height', 'comment', '');
nc_varput(ncfile, 'wave_height', Hs);

var = struct(...
    'Name', 'crs',...
    'Datatype', 'int32',...
    'Dimension', {{}}); % note the double brackets
nc_addvar(ncfile, var);
nc_attput(ncfile, 'crs', 'grid_mapping_name', 'oblique_stereographic');
nc_attput(ncfile, 'crs', 'epsg_code', 'EPSG:28992');
nc_attput(ncfile, 'crs', 'semi_major_axis', 6377397.155);
nc_attput(ncfile, 'crs', 'semi_minor_axis', 6356078.96282);
nc_attput(ncfile, 'crs', 'inverse_flattening', 299.1528128);
nc_attput(ncfile, 'crs', 'latitude_of_projection_origin', 52.0922178);
nc_attput(ncfile, 'crs', 'longitude_of_projection_origin', 5.23155);
nc_attput(ncfile, 'crs', 'scale_factor_at_projection_origin', 0.9999079);
nc_attput(ncfile, 'crs', 'false_easting', 155000.0);
nc_attput(ncfile, 'crs', 'false_northing', 463000.0);
nc_attput(ncfile, 'crs', 'proj4_params', '+proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.999908 +x_0=155000 +y_0=463000 +ellps=bessel +units=m +towgs84=565.4174,50.3319,465.5542,-0.398957388243134,0.343987817378283,-1.87740163998045,4.0725 +no_defs');
end

function [xrd,yrd] = argus2rd(station,xa,ya)
switch lower(station)
    case 'kijkduin'
        orig = [75130.3,454181];
        rot = 138;
    case 'zandmotor'
        orig = [72502,452074];
        rot = 140;
    otherwise
        error('Station not supported');
end

if isvector(xa)
    [xa,ya] = meshgrid(xa,ya);
end

alpha = rot*pi/180;

xrd = xa.*cos(alpha) - ya.*sin(alpha);
yrd = xa.*sin(alpha) + ya.*cos(alpha);

xrd = xrd - orig(1);
yrd = yrd - orig(2);
end