function [V, z1] = jarkus_verticalVolume(x, z, x1, x2, zref)
% jarkus_relativeVolume
%
%   Input:
%     x    = vector with specified distances
%     z    = vector with depths at locations x
%     x1   = landward boundary
%     x2   = seaward boundary
%     zref = lower boundary
%
%   Output:
%     V = volume in m2 relative to zref level (volume below this level is
%     treated as 'negative' volume)
%
% See also: JARKUS, snctools, jarkus_depthatX

%% Copyright notice
%   --------------------------------------------------------------------
%   Copyright (C) 2010 Alterra
%       Tommer Vermaas
%       tommer.vermaas@wur.nl | 0317-481654
%
%       Postbus 47
%       6700 AA Wageningen
%
%   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/>.
%   --------------------------------------------------------------------

%% Get position of upper and lower boundaries
z1 = jarkus_depthatX(x1,z,x);
z2 = jarkus_depthatX(x2,z,x);

if isnan(z1) || isnan(z2) %check if profile covers asked depth zone
    V = NaN;
else
    b  = find(x>=x1,1,'first');
    e  = find(x<=x2,1,'last');
    zm = min([z(b:e)' zref]);
    X  = [x1 x1 x(b:e)' x2 x2]';
    Z  = [zm z1 z(b:e)' z2 zm]';
    
    V=polyarea(X,Z)-((zref-zm)*(x2-x1));
end