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

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

%%
OPT = struct(...
    'maindir', fullfile(fileparts(which(mfilename)), '..', 'simulations'),...
    'regexp_select', '\d{8}_\d{8}.*?run02.*?xb(output)\.nc',...
    'regexp_id1', '(\d{8})(?=_)',...
    'regexp_id2', '(?<=_)\d{8}',...
    'save', false ...
    );

OPT = setproperty(OPT, varargin);

%%
[dirs dircont files] = dirlisting(OPT.maindir);

fid = ~cellfun(@isempty, (cellfun(@(x) regexp(x, OPT.regexp_select), files, 'uniformoutput', false)));

ncfiles = files(fid);

for incfile = 1:length(ncfiles)
    xb(incfile) = xb_read_output(ncfiles{incfile}, 'vars', {'globalx', 'zb'});
    x{incfile} = xs_get(xb(incfile), 'DIMS.x');
    zb{incfile} = squeeze(xs_get(xb(incfile), 'zb'));
%     figure;
%     plot(x{incfile}, zb{incfile}([1 end], :))
    [dummy result(incfile)] = jarkus_getVolume(-x{incfile}, zb{incfile}(1, :),...
        'x2', -x{incfile}, 'z2', zb{incfile}(end, :),...
        'LowerBoundary', 5);
    EV(incfile) = result(incfile).Volumes.Erosion;
    id1(incfile) = str2double(regexp(ncfiles{incfile}, OPT.regexp_id1, 'match'));
    id2(incfile) = str2double(regexp(ncfiles{incfile}, OPT.regexp_id2, 'match'));
end

refids = id1 == id2;

%% derive corresponding DUROS erosion volumes
for i = 1:length(id1)
    result = DUROS(-x{i}, zb{i}(1, :)', 225e-6, 5, 9, 12);
    EVd(i) = result(2).Volumes.Erosion;
end

%% derive relative erosion volumes (both of upper part and lower part influences)
% [sens1 sens2 sensd1 sensd2] = deal(ones(size(id1)));
% figure
% hold on
% for id = unique([id1(:); id2(:)])'
%     sens1(id1 == id) = EV(id1 == id) / EV(id1 == id & refids);
%     plot(sens1(id1 == id), 'b')
%     sens2(id2 == id) = EV(id2 == id) / EV(id1 == id & refids);
%     plot(sens2(id2 == id), 'r')
%     sensd1(id1 == id) = EVd(id1 == id) / EVd(id1 == id & refids);
%     plot(sensd1(id1 == id), ':b')
%     sensd2(id2 == id) = EVd(id2 == id) / EVd(id1 == id & refids);
%     plot(sensd2(id2 == id), ':r')
% end

array2Rstr = @(label,data,format) sprintf(['%s=c(' repmat([format ','],1,length(data)-1), format ')'], label, data);

sprintf('a <- data.frame(\n%s,\n%s,\n%s)',...
     array2Rstr('landward_part', id1, '''%i'''),...
     array2Rstr('seaward_part', id2, '''%i'''),...
     array2Rstr('erosion_volume', EV, '%.2f') ...
     )

%% plot series of profiles with identical upper profiles
c = 'rgb';
for id = unique(id1)
    figure;
    hold on
    
    for i = find(id1 == id)
        plot(x{i}, zb{i}(1, :), c(mod(i-1,length(c))+1),...
            'DisplayName', sprintf('%i', id2(i)))
        plot(x{i}, zb{i}(end, :), [c(mod(i-1,length(c))+1) ':'],...
            'DisplayName', sprintf('%i; EV = %.0f', id2(i), EV(i)))
    end
    title(sprintf('%i', id))
    legend show
    set(gca, 'xdir', 'reverse')
%     EV(id1 == id) / EV(id1 == id & id1 == id2)
end

%% plot matrix diagram with circles indicating the erosion volumes
id = unique([id1(:); id2(:)]);
[n1 n2] = deal(length(unique(id1)), length(unique(id2)));
n = n1*n2;
% create combinations
[id2id id1id] =  ind2sub([n2 n1], 1:n);
r = sqrt(EV/pi);
r = r/max(r)/2.2;

figure;
set(gca,...
    'xlim',  minmax(id1id) + [-.5 .5],...
    'ylim',  minmax(id2id) + [-.5 .5],...
    'xtick', unique(id1id),...
    'ytick', unique(id2id),...
    'xticklabel', cellfun(@num2str, num2cell(unique(id1)), 'UniformOutput', false),...
    'yticklabel', cellfun(@num2str, num2cell(unique(id2)), 'UniformOutput', false))
hold on
xlabel('Transect id landward part')
ylabel('Transect id seaward part')
for i = 1:n
    [x y] = cylinder(r(i),200);
    ph(i) = plot(x(1,:)+id1id(i), y(1,:)+id2id(i));
    th(i) = text(id1id(i), id2id(i), num2str(EV(i), '%.0f'),...
        'horizontalalignment', 'center');
end
axis square

figureStyle(@styleThesis,...
    {'Figure' 'Position' [1 1 15 8]});

%% plot relative erosion as function of the seaward part (id1), for distict landward parts
rgb = cellfun(@str2rgb, strread(c(1:3), '%1s'), 'UniformOutput', false);
refEV = EV(refids);
set(ph, {'color'}, rgb(id1id))
set(th, {'string'}, cellfun(@(x) sprintf('%.0f/%.0f', x), num2cell([EV; refEV(id1id)], 1), 'UniformOutput', false)')

figure;
set(gca,...
    'xlim',  minmax(id1id) + [-.5 .5],...
    ...'ylim',  minmax(id2id) + [-.5 .5],...
    'xtick', unique(id1id),...
    ...'ytick', unique(id2id),...
    'xticklabel', cellfun(@num2str, num2cell(id2(unique(id1id))), 'UniformOutput', false))%,...
    ...'yticklabel', cellfun(@num2str, num2cell(unique(id2)), 'UniformOutput', false))
xlabel('Transect id seaward part')
hold on
phs = plot(unique(id2id), reshape(EV./ refEV(id1id), 3,3), 'o');
set(phs, {'color'}, rgb)

legtxt = cellfun(@(x) sprintf('Landward part %i', x), num2cell(unique(id2)), 'UniformOutput', false);
legend(legtxt{:})

%% plot relative erosion as function of the landward part (id2), for distinct seaward parts
set(ph, {'color'}, rgb(id2id))
set(th, {'string'}, cellfun(@(x) sprintf('%.0f/%.0f', x), num2cell([EV; refEV(id2id)], 1), 'UniformOutput', false)')

figure
set(gca,...
    'xlim',  minmax(id2id) + [-.5 .5],...
    ...'ylim',  minmax(id2id) + [-.5 .5],...
    'xtick', unique(id2id),...
    ...'ytick', unique(id2id),...
    'xticklabel', cellfun(@num2str, num2cell(id2(unique(id2id))), 'UniformOutput', false))%,...
    ...'yticklabel', cellfun(@num2str, num2cell(unique(id2)), 'UniformOutput', false))
xlabel('Transect id landward part')
hold on
phs = plot(unique(id1id), reshape(EV./ refEV(id2id), 3,3)', 'o');
set(phs, {'color'}, rgb)

% phs = cellfun(@(x,y) plot(x,y,'o'), num2cell(id1id), num2cell(EV./ refEV(id2id)));
% set(phs, {'color'}, rgb(id2id))

legtxt = cellfun(@(x) sprintf('Seaward part %i', x), num2cell(unique(id1)), 'UniformOutput', false);
legend(legtxt{:})

%%
if OPT.save
    fpath = fileparts(which(mfilename));
    fname = mfilename;
    figurepath = fpath;
    set(gcf, 'tag', '$HeadURL$ $Revision$')
    print(gcf, '-depsc2', fullfile(figurepath, [fname '.eps']))
    if ispc
        saveas(gcf, fullfile(figurepath, [fname '.emf']))
    end
    saveas(gcf, fullfile(figurepath, [fname '.fig']))
end