function varargout = decc_plot_angle_time_erosion(varargin)
%DECC_PLOT_ANGLE_TIME_EROSION  One line description goes here.
%
%   More detailed description goes here.
%
%   Syntax:
%   varargout = decc_plot_angle_time_erosion(varargin)
%
%   Input:
%   varargin  =
%
%   Output:
%   varargout =
%
%   Example
%   decc_plot_angle_time_erosion
%
%   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: 29 Aug 2012
% Created with Matlab version: 7.13.0.564 (R2011b)

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

%%
OPT = struct(...
    'maindir', fullfile(fileparts(which(mfilename)), '..', 'simulations'),...
    'regexp_select', 'run(40|41).*1900.*xb_sum.nc',...
    'regexp_legend', 'run\d{2}.*?radius\d{4,5}(.*?rev\d{4}\D?)?',...
    'strrep_legend', {{{'\' '; '} {'s0' 's '} {'00' '00 m'}}},...
    'offset', 2,...
    'n_contours', 5,...
    'save', false,...
    'style', @styleThesis ...
    );

OPT = setproperty(OPT, varargin);

%%
[~, ~, files] = dirlisting(OPT.maindir);

ncfiles = files(~cellfun(@isempty, regexp(files, OPT.regexp_select)));

legtxt = cellfun(@(x) regexp(x, OPT.regexp_legend, 'match'), ncfiles);
for ilegrep = 1:length(OPT.strrep_legend)
    legtxt = cellfun(@(x) strrep(x, OPT.strrep_legend{ilegrep}{:}), legtxt,...
        'UniformOutput', false);
end
%%
if isscalar(ncfiles)
    fprintf('%s\n', ncfiles{:});
else
    fparts = cellfun(@(x) textscan(x, '%s', 'delimiter', repmat(filesep, 1, ispc+1)), ncfiles, 'UniformOutput', false);
    fparts = [fparts{:}];
    for i = 1:min(cellfun(@length, fparts))-1
        if ~isscalar(unique(cellfun(@(x) x{i}, fparts, 'UniformOutput', false)))
            fprintf('%s\n', fullfile(fparts{1}{1:i-1}));
            for j = 1:length(ncfiles)
                fprintf('\t%s\n', fullfile(' ', fparts{j}{i:end}));
            end
            break
        end
    end
end

%%
zb1 = cellfun(@(ncfile) nc_varget(ncfile, 'zb'), ncfiles, 'uniformoutput', false);
[tdim ydim xdim] = cellfun(@(z) size(z), zb1);
zb0 = cellfun(@(z) repmat(z(1,:,:), [size(z,1) 1 1]), zb1, 'uniformoutput', false);
[x y] = cellfun(@(ncfile) deal(nc_varget(ncfile, 'globalx'), nc_varget(ncfile, 'globaly')), ncfiles, 'uniformoutput', false); 
g = cellfun(@xb_stagger, x, y);
A = cellfun(@(x,y,tdim) repmat(shiftdim(x'./y', -1), [tdim 1 1]), {g.dsdnz}, {g.dnz}, num2cell(tdim), 'uniformoutput', false);
Ero = cellfun(@(zb0,zb1,A) sum(((max(zb1, 5) - max(zb0, 5))) .* A,3), zb0, zb1, A, 'uniformoutput', false);
% optional calculate the difference in erosion between 2 simulations
% Ero = cellfun(@(E0,E1) E0 - E1, Ero, fliplr(Ero), 'uniformoutput', false);
maxEro = cellfun(@(Ero) max(-Ero(:)), Ero);
phirad = cellfun(@(alfaz,tdim) repmat(pi-mean(alfaz,1), [tdim 1]), {g.alfaz}, num2cell(tdim), 'uniformoutput', false);
id = cellfun(@(phirad) phirad(1,:) > pi/2 & phirad(1,:) < 1.5*pi, phirad, 'uniformoutput', false);
for iradii = 1:length(ncfiles)
    ncfile = ncfiles{iradii};
    try
        t = nc_varget(ncfile, 'globaltime');
        t = t - min(t(:));
    catch
        t = linspace(0, 18000, tdim(iradii));
    end
    
    t = repmat(t(:), [1 ydim(iradii)])/3600+OPT.offset;
    [X,Y] = pol2cart(phirad{iradii},t);
    
    figure
    pcolor(X(:,id{iradii}), Y(:,id{iradii}), flipud(-Ero{iradii}(:,id{iradii})))
    title(legtxt{iradii})
    hold on
    contour(X(:,id{iradii}), Y(:,id{iradii}), flipud(-Ero{iradii}(:,id{iradii})), OPT.n_contours, 'w')
    shading interp
    c = colorbar;
    caxis([0 max(maxEro)])
    colormap(flipud(colormap('pink')))
    axis image
    yticklabel = (min(t(:)):max(t(:))) - OPT.offset;
    yticklabel = [yticklabel fliplr(yticklabel)];
    ytick = OPT.offset:max(t(:));
    ytick = sort([-ytick ytick]);
    set(gca,...
        'YAxisLocation', 'right',...
        'ytick', ytick,...
        'yticklabel', cellfun(@num2str, num2cell(yticklabel), 'uniformoutput', false),...
        'xtick', [])
    x = num2cell(-(max(t(:)):-1:min(t(:))));
    y = mean(minmax(Y(:)'));
    T = num2cell((min(t(:)):max(t(:)))-OPT.offset);
    cellfun(@(x,T) text(x,y, num2str(T),...
        'horizontalalignment', 'center',...
        'verticalalignment', 'bottom'), x, T)
    text(mean(-t(:)), y, 'time [h] \rightarrow',...
        'horizontalalignment', 'center',...
        'verticalalignment', 'top')
    set(get(c,'title'),'string','Erosion [m^3/m]');
    a = -90:45:90;
    r = minmax(t(:)')';
    x = -cosd(repmat(a,size(r))) .* repmat(r, size(a));
    y = sind(repmat(a,size(r))) .* repmat(r, size(a));
    plot(x,y,':k')
    cellfun(@(x,y,a) text(x,y,num2str(a),...
        'horizontalalignment', 'center',...
        'verticalalignment', 'middle'), num2cell(x(2,:)*1.07), num2cell(y(2,:)*1.07), num2cell(a))
    box off
    set(get(c, 'title'),'Position',[7.5 120 1],'rotation', 90)
    
    if OPT.save
    figureStyle(@styleThesis,...
        {'figure', 'Units', 'Centimeters', 'Position', [1 1 12 12]});
        [fpath fname] = fileparts(which(mfilename));
        fname = sprintf('%s_%s', fname, strrep(strrep(legtxt{iradii}, ' ', '_'), ';', ''));
        figurepath = fullfile(fpath, '..', 'figures');
        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
end