function z = debn_x2zfunction(varargin)
%DEBN_X2ZFUNCTION  One line description goes here.
%
%   More detailed description goes here.
%
%   Syntax:
%   varargout = debn_x2zfunction(varargin)
%
%   Input:
%   varargin  =
%
%   Output:
%   varargout =
%
%   Example
%   debn_x2zfunction
%
%   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: 14 Mar 2012
% Created with Matlab version: 7.7.0.471 (R2008b)

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

%%
OPT = struct(...
    'xInitial', [],...
    'zInitial', [],...
    'id', [],...
    'year', [],...
    'profileix', [],...
    'D50', [],...
    'WL_t', [],...
    'Hsig_t', [],...
    'Tp_t', [],...
    'filename', 'testcase.cas');

OPT = setproperty(OPT, varargin);

%%
[x z id yr] = getprofiles(OPT.id, OPT.year);
pfilter = cellfun(@isempty, x);
[x z id yr] = deal(x(~pfilter), z(~pfilter), id(~pfilter), yr(~pfilter));
idu = id(ceil(OPT.profileix*length(id)));
yru = yr(ceil(OPT.profileix*length(id)));

[Islope Ivolume duneHeight beachwidth EV EVsvp rWL] = deal(NaN(size(OPT.D50)));
[xInitial zInitial] = deal(cell(size(OPT.D50)));
for i = 1:length(x)
    [Islopetmp Ivolumetmp duneHeighttmp beachwidthtmp] = debn_profileindicators(x{i}, z{i});
    if ~isnan(yr(i))
        ix = ismember(idu, id(i)) & ismember(yru, yr(i));
    else
        ix = ismember(idu, id(i)) & isnan(yru);
    end
    xInitial(ix) = deal(x(i));
    zInitial(ix) = deal(z(i));
    [Islope(ix) Ivolume(ix) duneHeight(ix) beachwidth(ix)] = deal(Islopetmp, Ivolumetmp, duneHeighttmp, beachwidthtmp);
end
ws = getFallVelocity('D50', OPT.D50);
HpTw = OPT.Hsig_t ./ OPT.Tp_t ./ ws;

DuneErosionSettings('set',...
    'AdditionalErosion', false,...
    'BoundaryProfile', false,...
    'DuneBreachCalculation', false);

for isample = 1:length(OPT.D50)
    result = DUROS(xInitial{isample}, zInitial{isample}, OPT.D50(isample), OPT.WL_t(isample), OPT.Hsig_t(isample), OPT.Tp_t(isample));
    if ~isempty(result(1).Volumes.Erosion)
        EV(isample) = result(1).Volumes.Erosion;
    end
    if ~isscalar(result)
        EVsvp(isample) = result(2).Volumes.Erosion;
    else
        EVsvp(isample) = NaN;
    end
    rWL(isample) = interp1([3 nanmax(zInitial{isample})], 0:1, OPT.WL_t(isample), 'linear', 'extrap');
end

rEVsvp = EVsvp ./ (duneHeight - OPT.WL_t).^2;

netica_write_case_file(OPT.filename,...
...    'id', idu,...
...    'year', yru,...
    'slope', Islope,...
    'volume', Ivolume,...
    'duneheight', duneHeight,...
    'beachwidth', beachwidth,...
    'HpTw', HpTw,...
    'rWL', rWL,...
    'EVsvp', EVsvp,...
    'rEVsvp', rEVsvp);

z = NaN(size(OPT.D50));

%%
function [x z ido yro] = getprofiles(id, yr)

refprof = ismember(id, 0);
[x z] = deal({});
[yro ido] = deal([]);
if any(~refprof)
    tr = jarkus_transects(...
        'id', id(~refprof),...
        'year', yr);
    % interpolate nans in cross-shore direction (default)
    tr = jarkus_interpolatenans(tr);
    [yrm idm] = meshgrid(tr.time, tr.id);
    for ix = 1:numel(yrm)
        [i j] = ind2sub(size(yrm), ix);
        ztmp = squeeze(tr.altitude(j, i, :));
        nnid = ~isnan(ztmp);
        x{ix} = tr.cross_shore(nnid)';
        z{ix} = ztmp(nnid);
        ido(ix) = idm(ix);
        yro(ix) = year(yrm(ix) + datenum(1970,1,1));
    end
end

if any(refprof)
    [xtmp ztmp] = SimpleProfile;
    x = [{xtmp} x];
    z = [{ztmp} z];
    ido = [0 ido];
    yro = [NaN yro];
end