function varargout = desb_durosplus_results2netcdf(varargin)
%DESB_DUROSPLUS_RESULTS2NETCDF  One line description goes here.
%
%   More detailed description goes here.
%
%   Syntax:
%   varargout = desb_durosplus_results2netcdf(varargin)
%
%   Input:
%   varargin  =
%
%   Output:
%   varargout =
%
%   Example
%   desb_durosplus_results2netcdf
%
%   See also 

%% Copyright notice
%   --------------------------------------------------------------------
%   Copyright (C) 2011 Delft University of Technology
%       Kees den Heijer
%
%       C.denHeijer@TUDelft.nl	
%
%       Faculty of Civil Engineering and Geosciences
%       P.O. Box 5048
%       2600 GA Delft
%       The Netherlands
%
%   This library is free software: you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation, either version 3 of the License, or
%   (at your option) any later version.
%
%   This library is distributed 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.  See the
%   GNU General Public License for more details.
%
%   You should have received a copy of the GNU General Public License
%   along with this library.  If not, see <http://www.gnu.org/licenses/>.
%   --------------------------------------------------------------------

% This tool is part of <a href="http://www.OpenEarth.eu">OpenEarthTools</a>.
% OpenEarthTools is an online collaboration to share and manage data and 
% programming tools in an open source, version controlled environment.
% Sign up to recieve regular updates of this function, and to contribute 
% your own tools.

%% Version <http://svnbook.red-bean.com/en/1.5/svn.advanced.props.special.keywords.html>
% Created: 07 Sep 2011
% Created with Matlab version: 7.7.0.471 (R2008b)

% $Id$
% $Date$
% $Author$
% $Revision$
% $HeadURL$
% $Keywords: $

%%
OPT = struct(...
    'id', 7003925,...
    'irun', 0,...
    'maindir', fullfile(fileparts(which(mfilename)), '..', 'simulations'));

OPT = setproperty(OPT, varargin{:});

%%
for itransect = 1:length(OPT.id)
    basepath = fullfile(OPT.maindir, sprintf('%08i', OPT.id(itransect)), 'durosplus', sprintf('run%02i', OPT.irun));
    matfiles = dir(fullfile(basepath, '*.mat'));
    nT = length(matfiles);
    
    X = [];
    [T WL xmin xtransition xmax Vtot Vssl] = deal(NaN(nT,1));
    [x z z2] = deal({});
    
    for iT = 1:nT
        T(iT) = str2double(matfiles(iT).name(1:4));
        load(fullfile(basepath, matfiles(iT).name))
        WL(iT) = result(1).info.input.WL_t;
        xmin(iT) = min(result(1).xActive);
        cs = cumsum(result(1).Volumes.volumes);
        transitionid = cs == min(cs);
        xtransition(iT) = mean(result(1).xActive(transitionid(1:end) | transitionid([1 1:end-1])));
        xmax(iT) = max(result(1).xActive);
        Vtot(iT) = result(1).Volumes.Erosion;
        Vssl(iT) = result(2).Volumes.Erosion;
        x{iT} = [result(1).xLand; result(1).xActive; result(1).xSea];
        z{iT} = [result(1).zLand; result(1).zActive; result(1).zSea];
        z2{iT} = [result(1).zLand; result(1).z2Active; result(1).zSea];
        X = unique([X; x{iT}(:)]);
    end
    
    Z = NaN(nT,2,length(X));
    for iT = 1:nT
        Z(iT,1,:) = interp1(x{iT}, z{iT}, X);
        Z(iT,2,:) = interp1(x{iT}, z2{iT}, X);
    end
    
    results2netcdf(basepath, OPT.id(itransect), T, xmin, xtransition, xmax, Vtot, WL, Vssl, X, Z)
end

%%
function results2netcdf(basepath, id, T, xmin, xtransition, xmax, Vtot, ssl, Vssl, X, Z)
outputfile = fullfile(basepath, sprintf('%08i_pp.nc', id));

nc_create_empty(outputfile);

nc_add_dimension(outputfile, 'time', length(T))
nc_add_dimension(outputfile, 'state', size(Z,2))
nc_add_dimension(outputfile, 'cross_shore', size(Z,3))

% add variable and its attributes
s.Name      = 'time';
s.Nctype    = nc_double;
s.Dimension = {'time'};
s.Attribute = struct('Name', {'standard_name', 'units', 'axis'}, ...
    'Value', {'time', 'days since 1970-01-01 00:00 +1:00', 'T'});
nc_addvar(outputfile, s);

s.Name      = 'x_min';
s.Nctype    = nc_double;
s.Dimension = {'time'};
s.Attribute = struct('Name', {'long_name', 'units'}, ...
    'Value', {'minimum x-coordinate of erosion/sedimentation zone', 'm'});
nc_addvar(outputfile, s);

s.Name      = 'x_transition';
s.Nctype    = nc_double;
s.Dimension = {'time'};
s.Attribute = struct('Name', {'long_name', 'units'}, ...
    'Value', {'x-coordinate of transition between erosion/sedimentation zone', 'm'});
nc_addvar(outputfile, s);

s.Name      = 'x_max';
s.Nctype    = nc_double;
s.Dimension = {'time'};
s.Attribute = struct('Name', {'long_name', 'units'}, ...
    'Value', {'maximum x-coordinate of erosion/sedimentation zone', 'm'});
nc_addvar(outputfile, s);

s.Name      = 'erosion_volume';
s.Nctype    = nc_double;
s.Dimension = {'time'};
s.Attribute = struct('Name', {'long_name', 'units', 'comment'}, ...
    'Value', {'total erosion volume', 'm^3/m', 'total erosion volume per running meter alongshore'});
nc_addvar(outputfile, s);

s.Name      = 'ssl';
s.Nctype    = nc_double;
s.Dimension = {'time'};
s.Attribute = struct('Name', {'long_name', 'units', 'comment'}, ...
    'Value', {'maximum storm surge level', 'm', 'maximum storm surge level above geoid (NAP)'});
nc_addvar(outputfile, s);

s.Name      = 'erosion_volume_ssl';
s.Nctype    = nc_double;
s.Dimension = {'time'};
s.Attribute = struct('Name', {'long_name', 'units', 'comment'}, ...
    'Value', {'erosion volume above maximum storm surge level', 'm^3/m', 'erosion volume above maximum storm surge level per running meter alongshore'});
nc_addvar(outputfile, s);

s.Name      = 'cross_shore';
s.Nctype    = nc_double;
s.Dimension = {'cross_shore'};
s.Attribute = struct('Name', {'long_name', 'units', 'comment'}, ...
    'Value', {'cross-shore coordinate', 'm', 'cross-shore coordinate relative to the rsp (rijks strand paal)'});
nc_addvar(outputfile, s);

s.Name      = 'altitude';
s.Nctype    = nc_double;
s.Dimension = {'time' 'state' 'cross_shore'};
s.Attribute = struct('Name', {'long_name', 'units', 'comment'}, ...
    'Value', {'surface_altitude', 'm', 'altitude above geoid (NAP), before (pre-storm) and after (post-storm) DUROS+ calculation'});
nc_addvar(outputfile, s);

% put values into variable
if isnan(T)
    nc_varput(outputfile, 'time',  T)
else
    nc_varput(outputfile, 'time',  datenum([T repmat([ 1 1 0 0 0],length(T),1)])-datenum(1970,1,1))
end
nc_varput(outputfile, 'x_min', xmin)
nc_varput(outputfile, 'x_transition', xtransition)
nc_varput(outputfile, 'x_max', xmax)
nc_varput(outputfile, 'erosion_volume', Vtot)
nc_varput(outputfile, 'ssl', ssl)
nc_varput(outputfile, 'erosion_volume_ssl', Vssl)
nc_varput(outputfile, 'cross_shore', X)
nc_varput(outputfile, 'altitude', Z)

% add attributes
history = sprintf('created on %s', datestr(now));
nc_attput(outputfile, nc_global, 'history', history);

producer = sprintf('$HeadURL$ $Revision$');
nc_attput(outputfile, nc_global, 'producer', producer);

description = sprintf('Postprocessed results of DUROS+ simulations');
nc_attput(outputfile, nc_global, 'description', description);