function varargout = desb_xbeach_pp_results(varargin)
%DESB_XBEACH_PP_RESULTS  One line description goes here.
%
%   More detailed description goes here.
%
%   Syntax:
%   varargout = desb_xbeach_pp_results(varargin)
%
%   Input:
%   varargin  =
%
%   Output:
%   varargout =
%
%   Example
%   desb_xbeach_pp_results
%
%   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: 09 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'),...
    'threshold', .95);

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

%%
T = NaN;
if OPT.id ~= 0 | OPT.id ~= 0 & OPT.irun ~= 1
    T = jarkus_transects('url', fullfile(fileparts(which(mfilename)), '..', 'data', 'jarkus.nc'),...
        'output', {'time'});
    T = T.time;
end

%%
for itransect = 1:length(OPT.id)
    basepath = fullfile(OPT.maindir, sprintf('%08i', OPT.id(itransect)), 'xbeach', sprintf('run%02i', OPT.irun));
    ncfile = fullfile(basepath, sprintf('%08i.nc', OPT.id(itransect)));
    
    [X ix] = sort(-nc_varget(ncfile, 'globalx'));
    Z = nc_varget(ncfile, 'zb');
    if ndims(Z) == 2
        Z = reshape(Z, 1, size(Z,1), size(Z,2));
    end
    Z = Z(:,:,ix);
    nT = size(Z,1);
    
    zs0 = nc_varget(ncfile, 'zs0', [0 0 0 0], [-1 1 1 1]);
    
    if size(Z,1) == 1
        cumvol = cumtrapz(X, abs(squeeze(diff(Z,1,2))))';
    else
        cumvol = cumtrapz(X, abs(squeeze(diff(Z,1,2))), 2);
    end
    
    activeid = cumvol > (1-OPT.threshold)/2 * repmat(cumvol(:,end),1,size(cumvol,2)) &...
        cumvol < mean([1 OPT.threshold]) * repmat(cumvol(:,end),1,size(cumvol,2));
    
    [Vtot Vssl xmin xtransition xmax] = deal(NaN(nT,1));
    
    for iT = 1:nT
        [V Vres] = jarkus_getVolume(X, squeeze(Z(iT,1,:)), 'x2', X, 'z2', squeeze(Z(iT,end,:)));
        Vtot(iT) = Vres.Volumes.Erosion;
        cs = cumsum(Vres.Volumes.volumes);
        transitionid = cs == min(cs);
        xtransition(iT) = mean(Vres.xActive(transitionid(1:end) | transitionid([1 1:end-1])));
        [V Vres] = jarkus_getVolume(X, squeeze(Z(iT,1,:)), 'x2', X, 'z2', squeeze(Z(iT,end,:)),...
            'LowerBoundary', zs0(iT));
        Vssl(iT) = Vres.Volumes.Erosion;
        xmin(iT) = X(find(activeid(iT,:), 1, 'first'));
        xmax(iT) = X(find(activeid(iT,:), 1, 'last'));
    end
    
   
    results2netcdf(basepath, OPT.id(itransect), T, xmin, xtransition, xmax, Vtot, zs0, 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', size(Z,1))
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
nc_varput(outputfile, 'time',  T)
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 XBeach simulations');
nc_attput(outputfile, nc_global, 'description', description);