function [z result] = x2zBreach(varargin)
%X2ZBREACH  limite state function for breaching of dune
%
%   The limit state is here defined as less than the boundary profile
%   volume is left above the storm surge level.
%
%   Syntax:
%   [z result] = x2zBreach(varargin)
%
%   Input:
%   varargin   =
%
%   Output:
%   z          = z-values based on the samples
%   result     = DUROS result
%
%   Example
%   x2zBreach
%
%   See also 

%% Copyright notice
%   --------------------------------------------------------------------
%   Copyright (C) 2010 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 Lesser General Public
%   License as published by the Free Software Foundation, either
%   version 2.1 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
%   Lesser General Public License for more details.
%
%   You should have received a copy of the GNU Lesser General Public
%   License along with this library. If not, see <http://www.gnu.org/licenses/>.
%   --------------------------------------------------------------------

% This tool is part of <a href="http://OpenEarth.nl">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: 26 Oct 2010
% Created with Matlab version: 7.10.0.499 (R2010a)

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

%%
OPT = struct(...
    'WL_t', [],...
    'Hsig_t', [],...
    'Tp_t', [],...
    'D50', [],...
    'Accuracy', [],...
    'Duration', [],...
    'xInitial', (-210:5:760)',...
    'zInitial', [3.66 3.21 4.42 6.86 8.73 11.34 11.33 11.86 13.65 15.41 17.34 19.11 18.18 16.69 15.51 14.12 13.36 13.51 13.62 13.29 12.82 12.03 11.05 10.14 8.71 6.89 3.63 2.86 2.33 1.96 1.62 1.38 1.16 0.99 0.8 0.62 0.58 0.45 0.29 0.12 0.06 -0.01 -0.11 -0.22 -0.38 -0.54 -0.57 -0.57 -0.61 -0.69 -0.76 -0.87 -0.97 -1.12 -1.27 -1.495 -1.72 -1.875 -2.03 -2.13 -2.23 -2.365 -2.5 -2.59 -2.68 -2.65 -2.62 -2.51 -2.4 -2.255 -2.11 -2.095 -2.08 -2.245 -2.41 -2.63 -2.85 -3.09 -3.33 -3.545 -3.76 -4.125 -4.49 -4.685 -4.88 -4.865 -4.85 -4.775 -4.7 -4.62 -4.54 -4.455 -4.37 -4.38 -4.39 -4.46 -4.53 -4.56 -4.59 -4.605 -4.62 -4.64 -4.66 -4.7 -4.74 -4.795 -4.85 -4.865 -4.88 -4.9 -4.92 -4.935 -4.95 -4.945 -4.94 -4.955 -4.97 -4.98 -4.99 -4.98 -4.97 -4.98 -4.99 -4.985 -4.98 -4.995 -5.01 -4.995 -4.98 -4.96 -4.94 -4.985 -5.03 -5.085 -5.14 -5.23 -5.32 -5.39 -5.46 -5.585 -5.71 -5.835 -5.96 -6.065 -6.17 -6.25 -6.33 -6.305 -6.28 -6.175 -6.07 -5.88 -5.69 -5.42 -5.15 -4.81 -4.47 -4.05 -3.63 -3.325 -3.02 -3.04 -3.06 -3.165 -3.27 -3.4 -3.53 -3.64 -3.75 -3.85 -3.95 -4.065 -4.18 -4.325 -4.47 -4.6 -4.73 -4.88 -5.03 -5.185 -5.34 -5.525 -5.71 -5.855 -6 -6.05 -6.1 -6.2 -6.3 -6.365 -6.43 -6.42 -6.41 -6.49 -6.57]',...
    'database', '',...
    'transectid', [],...
    'table', '');

OPT = setproperty(OPT, varargin{:});
    
%% cross-shore profile
xInitial = OPT.xInitial;
zInitial = OPT.zInitial;

%% samples
samples.WL_t = OPT.WL_t;
samples.Hsig_t = OPT.Hsig_t;
samples.Tp_t = OPT.Tp_t;
samples.D50 = OPT.D50;
samples.Accuracy = OPT.Accuracy;
samples.Duration = OPT.Duration;

writeerrors = ~isempty(OPT.database);

%% retrieve calculation values
z = ones(size(samples.D50,1),1); % pre-allocate z
[VolumeAboveSSL ErosionAboveSSL] = deal(zeros(size(z)));
restVolume = NaN(size(z));

for i = 1:size(samples.D50,1)
    try
        if samples.WL_t(i) < max(zInitial)
            %% set calculation values for additional volume
            DuneErosionSettings('set',...
                'AdditionalVolume', [num2str(samples.Duration(i)) '*Volume + ' num2str(samples.Accuracy(i)) '*Volume'],... string voor het bepalen van het toeslagvolume gedurende de berekening (afslagvolume is negatief)
                'BoundaryProfile', false,...       % Grensprofiel berekenen is niet nodig, gebruiken we niet
                'DuneBreachCalculation', false);
            
            %% calculate dune volume above current storm surge level
            VolumeAboveSSL(i) = getVolume(xInitial, zInitial, 'LowerBoundary', samples.WL_t(i));
            
            %% carry out DUROS+ computation
            [result messages] = DUROS(xInitial, zInitial, samples.D50(i), samples.WL_t(i), samples.Hsig_t(i), samples.Tp_t(i));
            if ~isscalar(result)
                samples.Tp_t(i) = result(1).info.input.Tp_t;
            end
            
            if writeerrors
                mksqlite('open', OPT.database);
                mksqlite('begin');
                nmessage = size(messages,1);
                for imessage = 1:nmessage
                    mksqlite(['insert into ' OPT.table ' (profile, sample, errorcode, message) values (' num2str(OPT.transectid) ', ' num2str(i) ', ' num2str(messages{imessage,1}) ', ''' messages{imessage,2} ''')']);
                end
                mksqlite('commit');
                mksqlite('close', OPT.database);
            end
            
            %% calculate erosion volume above storm surge level
            if ~isscalar(result)
                [Erosion Accretion] = deal(0);
                for iresult = 2:length(result)
                    Erosion = Erosion + result(iresult).Volumes.Erosion;
                    Accretion = Accretion + result(iresult).Volumes.Accretion;
                end
                ErosionAboveSSL(i) = Erosion - Accretion;
                restVolume(i) = VolumeAboveSSL(i) - ErosionAboveSSL(i);
            end
        else
            restVolume(i) = -Inf;
        end
        
        % calculate boundary profile volume
        ProfileHeight = max([2.5 0.12 * samples.Tp_t(i) * sqrt(samples.Hsig_t(i))]);
        boundaryprofilevolume = 3/2 * ProfileHeight^2 + 3*ProfileHeight;
        
%         figure;
%         plotDuneErosion(result)

        %% Derive whether dune is breached
        % breaching is here defined as "less than the boundaryprofile
        % volume is left above SSL"
        Breach(i) = restVolume(i) <= boundaryprofilevolume;
        if Breach(i)
            z(i,:) = -Inf;
        else
            z(i,:) = restVolume(i);
        end

    catch
        ME = lasterror;
        if strcmp(ME.identifier, 'DUROSCHECKCONDITIONS:lowprofile')
            Breach(i) = true;
        else
            Breach(i) = false;
            if writeerrors
                mksqlite('open', OPT.database);
                mksqlite('begin');
                mksqlite(['insert into ' OPT.table ' (profile, sample, errorcode, message) values (' num2str(OPT.transectid) ', ' num2str(i) ', -999, ''Fatal error'')']);
                mksqlite('commit');
                mksqlite('close', OPT.database);
            end
        end
    end

end