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

%% Copyright notice
%   --------------------------------------------------------------------
%   Copyright (C) 2014 Deltares
%       Kees den Heijer
%
%       kees.denheijer@deltares.nl
%
%       P.O. Box 177
%       2600 MH 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: 26 Sep 2014
% Created with Matlab version: 8.2.0.701 (R2013b)

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

%%
OPT = struct(...
    );
% return defaults (aka introspection)
if nargin==0;
    varargout = {OPT};
    %    return
end
% overwrite defaults with user arguments
OPT = setproperty(OPT, varargin);

%% code
resultsdir = fullfile(fileparts(mfilename('fullpath')), '..', 'results');
figuredir = fullfile(fileparts(mfilename('fullpath')), '..', 'figures');

D = dir2(resultsdir,...
    'file_incl', '\.nc$',...
    'no_dirs', true,...
    'dir_incl', '\d{7,8}_\d+$');
fullfilenames = strcat({D.pathname}, {D.name})';
% fullfilenames = [fullfilenames fullfilenames];
tokennames = regexp(fullfilenames, '(?<id>\d{7,8})_(?<year>\d+)', 'names');
metainfo = cell2mat(tokennames);
ids = cellfun(@str2double, {metainfo.id});
years = cellfun(@str2double, {metainfo.year});
uids = unique(ids','rows','stable');
for i = 1:length(uids)
    fnames = fullfilenames(ids==uids(i));
    hFig = figure(i);
    n = length(fnames);
    hAx1 = subplot(3,1,1);
    %     hAx2 = subplot(3,1,2);
    %     hAx3 = subplot(3,1,3);
    hAxD = subplot(3,1,2:3);
    hold(hAx1,'on');
    hold(hAxD,'on');
    %     hold(hAx2,'on');
    %     hold(hAx3,'on');
    czb0 = ['k--';'r--';'c--'];
    czb = ['k-';'r-';'c-'];
    ch  = ['r:';'r-';'k-'];
    czs = ['b:';'b-';'k-'];
    for j = 1:n
        xb = xb_read_output(fnames{j}, 'vars', {'H', 'zb', 'zs', 'globalx'});
        x = nc_varget(fnames{j},'globalx');%xs_get(xb, 'DIMS:globalx');
        H = xs_get(xb, 'H');
        zb = xs_get(xb, 'zb');
        zs = xs_get(xb, 'zs');
        tide = load([D(n*(i-1)+j).pathname 'tide.txt']);
        zb0 = squeeze(zb(1,1,:));
        zb1 = squeeze(zb(end,1,:));
        plot(hAx1,x,zb0,czb0(j,:),'DisplayName',['Initial bed ' num2str(metainfo(n*(i-1)+j).year)])
        %         plot(hAx1,x,zb1,czb(j,:),'DisplayName',['Final bed ' num2str(metainfo(2*(i-1)+j).year)])
        
        % Determine retreat point in the same waty Duros+ additional
        % erosion is determined
        plot(hAxD,x,zb0,czb0(j,:),'DisplayName',['Initial bed ' num2str(metainfo(n*(i-1)+j).year)])
        plot(hAxD,x,zb1,czb(j,:),'DisplayName',['Final bed ' num2str(metainfo(n*(i-1)+j).year)])
        hline(tide(1,2),'b:')
        [~, volres, ~] = getVolume('x',x,'z',zb0,'x2',x,'z2',zb1,'LowerBoundary',tide(1,2));
        TargetVolume = -volres.Volumes.Erosion;
        if ~isempty(zb1(zb1>tide(1,2)))
            ErosionResult = getAdditionalErosion(x,zb0, ...
                'TargetVolume', TargetVolume, ...
                'poslndwrd',true, ...
                'x0min', min(x), ...
                'x0max', max(findCrossings(x,zb1,[min(x),max(x)],ones(1,2)*tide(1,2))), ...
                'zmin',tide(1,2));
            
            plot(hAxD,ErosionResult.VTVinfo.Xr,ErosionResult.VTVinfo.Zr,[czb(j,1) 'o'],'MarkerFaceColor',czb(j,1),'DisplayName',['R point ' num2str(metainfo(n*(i-1)+j).year) ', V_{e} = ' num2str(roundoff(-TargetVolume,-1)) 'm^{3}/m'])
            plot(hAxD,ErosionResult.xActive,ErosionResult.z2Active,[czb(j,1) ':'],'DisplayName','Proxy erosion volume Duros+')
            hline(tide(1,2),'b:')
            % get retreat distance
            RD(j) = diff(minmax(ErosionResult.xActive'));
        end
        %         plot(hAx2,x,squeeze(H(end,1,:)),ch(j,:),'DisplayName',['Wave height ' num2str(metainfo(2*(i-1)+j).year)])
        %
        %         plot(hAx3,x,squeeze(zs(end,1,:)),czs(j,:),'DisplayName',['Water level ' num2str(metainfo(2*(i-1)+j).year)])
    end
    % print to screen
    fprintf('transect %8s\n', metainfo(n*(i-1)+j).id)
    fprintf('\tyear %s: retreat distance = %.0f\n', metainfo(n*(i-1)+1).year, RD(1))
    fprintf('\tyear %s: retreat distance = %.0f\n', metainfo(n*(i-1)+2).year, RD(2))
    fprintf('\tyear %s: retreat distance = %.0f\n', metainfo(n*(i-1)+3).year, RD(3))
    fprintf('\t\tverschil: %.0f\n', diff(RD))
    if isempty(ErosionResult)
        xlim(hAxD,[max(x)-800 max(x)])
    else
        xlim(hAxD,[ErosionResult.VTVinfo.Xr-300 ErosionResult.VTVinfo.Xr+500])
    end
    ylim(hAxD,[0 max(zb0)+2])
    title(hAx1, ['Transect ' num2str(uids(i))])
    xlabel(hAx1,'Cross-shore distance [m]')
    xlabel(hAxD,'Cross-shore distance [m]')
    ylabel(hAx1,'Elevation [m + NAP]', 'rotation', 270, 'VerticalAlignment', 'top')
    ylabel(hAxD,'Elevation [m + NAP]', 'rotation', 270, 'VerticalAlignment', 'top')
    hL1 = legend(hAx1,'show','location','NorthWest');
    hL2 = legend(hAxD,'show','location','NorthEast');
    set(hL2,'Position',[0.80 0.40 0.1 0.1])
    print(fullfile(figuredir,num2str(uids(i))),'-dpng','-r600')
    close(hFig)
    ErosionResult = [];
    %     legend(hAx2,'show','location','SouthEast')
    %     legend(hAx3,'show','location','SouthEast')
end