function varargout = debn_generate_decc_dewa_data(varargin)
%DEBN_GENERATE_DECC_DEWA_DATA  One line description goes here.
%
%   More detailed description goes here.
%
%   Syntax:
%   varargout = debn_generate_decc_dewa_data(varargin)
%
%   Input: For <keyword,value> pairs call debn_generate_decc_dewa_data() without arguments.
%   varargin  =
%
%   Output:
%   varargout =
%
%   Example
%   debn_generate_decc_dewa_data
%
%   See also

%% Copyright notice
%   --------------------------------------------------------------------
%   Copyright (C) 2012 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: 19 Sep 2012
% Created with Matlab version: 7.12.0.635 (R2011a)

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

%%
OPT = struct(...
    ...'maindir', 'F:\Sources\OpenEarthModels\tudelft\dune_erosion_coastline_curvature\simulations\',...
    'maindir_decc', 'd:\heijer\sources\dune_erosion_coastline_curvature\simulations\00000000\xbeach\',...
    'maindir_dewa', 'd:\heijer\sources\dune_erosion_wave_angle\simulations\00000000\xbeach\',...
    'destdir', fullfile(fileparts(which(mfilename)), '..', '..', 'dune_erosion_bayesian_network', 'netica'),...
    'regexp_select_decc', 'run(40).*?0\\xb_sum\.nc',...
    'regexp_select_dewa', 'run(09).*?0\\xb_sum\.nc',...
    'save', false ...
    );

%% wave angle
[~, ~, files] = dirlisting(OPT.maindir_dewa);

ncfiles = files(~cellfun(@isempty, regexp(files, OPT.regexp_select_dewa)));

zb1 = cellfun(@(ncfile) nc_varget(ncfile, 'zb', [0 0 0], [2 -1 -1], [90 1 1]), ncfiles, 'uniformoutput', false);
[tdim ydim xdim] = cellfun(@(z) size(z), zb1);
zb0 = cellfun(@(z) repmat(z(1,:,:), [size(z,1) 1 1]), zb1, 'uniformoutput', false);
[x y] = cellfun(@(ncfile) deal(nc_varget(ncfile, 'globalx'), nc_varget(ncfile, 'globaly')), ncfiles, 'uniformoutput', false); 
g = cellfun(@xb_stagger, x, y);
A = cellfun(@(x,y,tdim) repmat(shiftdim(x'./y', -1), [tdim 1 1]), {g.dsdnz}, {g.dnz}, num2cell(tdim), 'uniformoutput', false);
Ero = cellfun(@(zb0,zb1,A) sum((max(zb1, 5) - max(zb0, 5)) .* A,3), zb0, zb1, A, 'uniformoutput', false);

idx_zonebounds = cellfun(@(y) interp1(y(:,1)/y(end,1), 1:size(y,1), [.3 .7], 'nearest'), y, 'uniformoutput', false);

Ero = cellfun(@(Ero,idx) -Ero(end,idx(1):idx(2)), Ero, idx_zonebounds, 'uniformoutput', false);

m = regexp(ncfiles, '(?<=angle)\d+', 'match');
angles = cellfun(@(x,y) repmat(str2double(x), 1, length(y)), [m{:}], Ero, 'uniformoutput', false);
phi_dewa = 270-[angles{:}];
ero_dewa = [Ero{:}];
radii_dewa = ones(size(ero_dewa)) * 2e4;

reference_erosion = mean(ero_dewa(phi_dewa==0));

%% curvature
[~, ~, files] = dirlisting(OPT.maindir_decc);

ncfiles = files(~cellfun(@isempty, regexp(files, OPT.regexp_select_decc)));

for iradii = 1:length(ncfiles)
    ncfile = ncfiles{iradii};
    
    result(iradii) = decc_analyse(ncfile,...
        'phi_resolution', nan,...
        's', false,...
        'zb', true,...
...        'Bal', false,...
...        'Fy_mean', false,...
...        'v_mean', false,...
...        'u_mean', false,...
...        'Svsg_cs', false,...
...        'ee_mean', true,...
        'contour', -2.2);
end

m = regexp(ncfiles, '(?<=radius)\d+', 'match');
radii_decc = cellfun(@(x,y) repmat(str2double(x), 1, y), [m{:}], num2cell(cellfun(@length, {result.Ero})), 'UniformOutput', false);
radii_decc = [radii_decc{:}];
phi_decc = abs([result.phi]);
ero_decc = [result.Ero];

%%
netica_write_case_file(fullfile(OPT.destdir, 'curvature_wave_angle.cas'),...
    'phi', [phi_dewa phi_decc],...
    'erosionfactor', [ero_dewa ero_decc]/reference_erosion,...
    'radius', [radii_dewa radii_decc]);
