function result = decc_plot_alongshore(varargin)
%DECC_PLOT_ALONGSHORE  One line description goes here.
%
%   More detailed description goes here.
%
%   Syntax:
%   varargout = decc_plot_alongshore(varargin)
%
%   Input:
%   varargin  =
%
%   Output:
%   varargout =
%
%   Example
%   decc_plot_alongshore
%
%   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: 20 Jan 2012
% Created with Matlab version: 7.7.0.471 (R2008b)

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

%%
OPT = struct(...
    'durosta', fullfile(fileparts(which(mfilename)), '..', '..', 'heijer', 'matlab', 'durosta', 'durosta_coastal_curvature.mat'),...
    'xbeach_dewa', fullfile(fileparts(which(mfilename)), '..', '..', 'dune_erosion_wave_angle', 'matlab', 'wave_angle_erosion_xbeach.mat'),...
    'maindir', fullfile(fileparts(which(mfilename)), '..', 'simulations'),...
    'regexp_select', 'run(40).*xb_sum.nc',...
    'regexp_legend', 'radius\d{4,5}(.*?rev\d{4}\D?)?',...
    'strrep_legend', {{{'\' '; '} {'s0' 's '} {'00' '00 m'}}},...
    'save', true ...
    );

OPT = setproperty(OPT, varargin);

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

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

%%
hsv1 = rgb2hsv([0 0 1]);
hsv2 = rgb2hsv([1 0 1]);
h = interp1(0:1, [hsv1(1) hsv2(1)], linspace(0,1,length(ncfiles)));
colors = cellfun(@(h) hsv2rgb([h .8 hsv1(3)]), num2cell(h),...
    'uniformoutput', false);

%%
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

%%

for iradii = 1:length(ncfiles)
    
    ncfile = ncfiles{iradii};
    
    result(iradii) = decc_analyse(ncfile,...
        's', true,...
        'zb', true,...
        'Bal', true,...
        'Fy_mean', false,...
        'v_mean', true,...
        'u_mean', true,...
        'ee_mean', false,...
        'contour', -2.2);
    % normalize alongshore transport gradient
%     result(iradii).dSvsg_dn = result(iradii).dSvsg_dn / max(result(iradii).dSvsg_dn);
end

if OPT.save
    % save result structure to mat file
    tag = '$HeadURL$ $Revision$';
    maindir =  abspath(OPT.maindir);
    regexp_select = OPT.regexp_select;
    save(fullfile(fileparts(which(mfilename)), 'decc_result.mat'), 'result', 'tag', 'maindir', 'regexp_select')
end

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

%%
% array2Rstr = @(label,data,format) sprintf(['%s=c(' repmat([format ','],1,length(data)-1), format ')'], label, data);
% phi = [result(1).phi];
% radius = [repmat(1900, 1, length(result(1).phi))] repmat(2400, 1, length(result(2).phi))];
% EV = [result(1).Ero];
% id = phi > 0;
% sprintf('a <- data.frame(\n%s,\n%s)',...
% ...     array2Rstr('radius', radius(id), '''%i'''),...
%      array2Rstr('angle', phi, '%.2f'),...
%      array2Rstr('erosion_volume', EV, '%.2f') ...
%      )

fh = [];
%%
erolimits = [min(cellfun(@min, {result.Ero})) max(cellfun(@max, {result.Ero}))];
sedlimits = [min(cellfun(@min, {result.Sed})) max(cellfun(@max, {result.Sed}))];
ballimits = [min(cellfun(@min, {result.Bal})) max(cellfun(@max, {result.Bal}))];
ydiff = max(cellfun(@diff, {erolimits sedlimits ballimits}));
erolimits = mean(erolimits) + ydiff/2 * [-1 1];
sedlimits = mean(sedlimits) + ydiff/2 * [-1 1];
ballimits = mean(ballimits) + ydiff/2 * [-1 1];
fh(end+1) = plot_engine({result.phi}, {result.Ero}, 'Erosion', 'm^3/m',...
    'figure', figure,...
    'axis', [-90 90 erolimits],...
    'colors', colors,...
    'ylabelFcn', @(varname,units) sprintf('%s volume [%s] above SSL', varname, units),...
    'xlabel', 'angle between shore normal and incident wave direction [degrees]',...
    'displayname', legtxt,...
    'title', false);
if ~isempty(OPT.durosta)
    load(OPT.durosta)
    curvature = arrayfun(@(x) x.settings.Params.RK, D);
    id = { mean(curvature,2)~=0,':'};
    wave_angle = mean(arrayfun(@(x) x.Input.WaveAngle, D(id{:})));
    erosion = arrayfun(@(x) x.Output.Volumes.Data(2), D(id{:}));
    erosion = [fliplr(erosion(:,2:end)) erosion];
    erosion = reshape(minmax(erosion'), 1, []);
    wave_angle = [fliplr(-wave_angle(:,2:end)) wave_angle];
    ax1 = findobj(fh(end), 'type', 'axes');
    set(ax1, 'nextplot', 'add');
    pah = patch([wave_angle fliplr(wave_angle)], erosion, ones(1,length(erosion)),...
        'parent', ax1,...
        'displayname', 'DurosTA',...
        'EdgeColor', [.6 .6 .6],...
        'FaceColor', [.6 .6 .6] ...
        );
    uistack(pah, 'bottom')
end
if ~isempty(OPT.xbeach_dewa)
    load(OPT.xbeach_dewa)
    hold on
%     % mirror wave directions
%     angles_mean_mir = unique([-angles_mean angles_mean]);
%     Ero_mean_mir = interp1(angles_mean, Ero_mean, abs(angles_mean_mir));
%     Ero_std_mir = interp1(angles_mean, Ero_std, abs(angles_mean_mir));
    errorbar(angles_mean, Ero_mean, Ero_std*2,...
        'color', 'r',...
        'marker', 'x',...
        'color', zeros(1,3),...
        'linestyle', 'none',...
        'displayname', 'no curvature');
    refid = angles_mean==0;
end
legend(ax1,...
    'location', 'southwest')
set(ax1,...
    'ylim', [100 260])

% ax2 = axes('Position',get(ax1,'Position'),...
%     'XAxisLocation','top',...
%     'YAxisLocation','right',...
%     'Color','none',...
%     'XColor','k','YColor','k');
% set(ax2, 'YLim', get(ax1, 'YLim')/Ero_mean(refid)*100,...
%     'XTick', [])
% ylabel(sprintf('relative erosion volume\nw.r.t. reference [%%]'))
% fix primary axis limits
% set([ax1 ax2], 'ylimmode', 'manual', 'xlimmode', 'manual')

if OPT.save
    figureStyle(@styleThesis, ...
        {'figure', 'Position', [1 1 15 9]});
    tightenfigure
    tightenfigure
%     set(fh(end),...
%         'Units', 'centimeters',...
%         'Position', [1 1 15 9],...
%         'PaperPositionMode', 'auto');
%     set(findobj(ax1, 'type', 'line'), 'linewidth', 2)
%     set(ax1,...
%         'Position', get(ax1, 'Position') .*[1 1 .98 1],...
%         'FontSize', 12)
%     set(ax2,...
%         'Position', get(ax1, 'Position'),...
%         'FontSize', 12)
%     set(findall(fh(end), 'type', 'text'),...
%         'FontSize', 12)
    [fpath fname] = fileparts(which(mfilename));
    fname = sprintf('%s_%s', fname, 'erosion');
    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

%%
fh(end+1) = plot_engine({result.phi}, {result.Ero}, 'Erosion', 'm^3/m',...
    'figure', figure,...
    'colors', colors,...
    'axes', subplot(3,1,1),...
    'axis', [-90 90 erolimits],...
    'xticklabel', false,...
    'displayname', legtxt,...
    'title', 'Erosion');
fh(end+1) = plot_engine({result.phi}, {result.Sed}, 'Accretion', 'm^3/m',...
    'figure', gcf,...
    'colors', colors,...
    'title', 'Accretion',...
    'axis', [-90 90 sedlimits],...
    'axes', subplot(3,1,2),...
    'xticklabel', false,...
    'displayname', legtxt,...
    'legend', true);
fh(end+1) = plot_engine({result.phi}, {result.Bal}, 'Balance', 'm^3/m',...
    'figure', gcf,...
    'colors', colors,...
    'axes', subplot(3,1,3),...
    'axis', [-90 90 ballimits],...
    'title', 'Accretion - Erosion',...
    'xlabel', 'angle between shore normal and incident wave direction [degrees]',...
    'save', OPT.save,...
    'suffix', 'Erosion Accretion');

%%
fh(end+1) = plot_engine({result.phi}, {result.Svsg_cs}, 'S', 'm^3/s',...
    'figure', figure,...
    'colors', colors,...
    'axes', subplot(2,1,1),...
    'axis', 'tight',...
    'xticklabel', false,...
    'title', sprintf('Integrated alongshore transport; above NAP%+.1f m', result(1).contour),...
    'legend', legtxt,...
    'legendlocation', 'SouthEast');
fh(end+1) = plot_engine({result.phi}, {result.dSvsg_dn}, 'S gradient', 'm^3/s/m',...
    'figure', gcf,...
    'colors', colors,...
    'title', 'Alongshore transport gradient',...
    'axis', 'tight',...
    'axes', subplot(2,1,2),...
    ...'xticklabel', false,...
    'xlabel', 'angle between shore normal and incident wave direction [degrees]',...
    'save', OPT.save,...
    'suffix', 'Alongshore Transport gradient');
% fh(end+1) = plot_engine({result.phi}, {result.ccg_cs}, 'Conc', 'kg/m^3',...
%     'figure', gcf,...
%     'axes', subplot(3,1,3),...
%     'axis', 'tight',...
%     'title', 'Maximum of mean sediment concentration',...
%     'xlabel', 'angle between shore normal and incident wave direction [degrees]',...
%     'save', OPT.save,...
%     'suffix', 'Transport Concentration');

%%
% fh(end+1) = plot_engine({result.phi}, {result.zs_mean_offshore}, 'water level', 'm w.r.t. NAP',...
%     'figure', figure,...
%     'axes', subplot(2,1,1),...
%     'axis', 'tight',...
%     'xticklabel', false,...
%     'title', 'Mean water level @ offshore boundary',...
%     'legend', legtxt);
fh(end+1) = plot_engine({result.phi}, {result.zs_mean_nearshore},  'water level', 'm w.r.t. NAP',...
    'figure', figure,...
    'colors', colors,...
    'axes', subplot(2,1,1),...
    'axis', 'tight',...
    'title', sprintf('Mean water level nearshore (@ NAP%+.1f m)', result(1).zs_mean_contour_nearshore),...
    'xlabel', false,...
    'xticklabel', false,...
    'legend', legtxt,...
    'save', false,...
    'suffix', 'Water levels');

fh(end+1) = plot_engine({result.phi}, {result.dzs_mean_nearshore_dnz},  'water level', 'm w.r.t. NAP',...
    'figure', gcf,...
    'colors', colors,...
    'axes', subplot(2,1,2),...
    'axis', 'tight',...
    'title', sprintf('Mean water level nearshore gradient (@ NAP%+.1f m)', result(1).zs_mean_contour_nearshore),...
    'xlabel', 'angle between shore normal and incident wave direction [degrees]',...
    'save', OPT.save,...
    'suffix', 'Water level gradient');

fh(end+1) = plot_engine({result.phi}, {result.Susg_cs},  'max cross-shore transport', 'm^3/s',...
    'figure', figure,...
    'colors', colors,...
    'axes', subplot(2,1,1),...
    'axis', 'tight',...
    'title', 'Maximum cross-shore transport',...
    'xlabel', false,...
    'xticklabel', false,...
    'legend', legtxt,...
    'save', false,...
    'suffix', 'max crossshore transport');
fh(end+1) = plot_engine({result.phi}, {result.ccg_cs}, 'Concentration', 'kg/m^3',...
    'figure', gcf,...
    'colors', colors,...
    'axes', subplot(2,1,2),...
    'axis', 'tight',...
    'title', 'Maximum of mean suspended sediment concentration',...
    'xlabel', 'angle between shore normal and incident wave direction [degrees]',...
    'save', OPT.save,...
    'suffix', 'Crossshore Transport Concentration');

%%
fh(end+1) = plot_engine({result.phi}, {result.v_mean}, 'mean alongshore velocity', 'm/s',...
    'figure', figure,...
    'colors', colors,...
    'axes', subplot(2,1,1),...
    'axis', 'tight',...
    'xticklabel', false,...
    'legend', legtxt,...
    'legendlocation', 'NorthWest',...
    'title', 'Mean alongshore velocity @ offshore boundary');
fh(end+1) = plot_engine({result.phi}, {result.u_mean}, 'mean cross-shore velocity', 'm/s',...
    'figure', gcf,...
    'colors', colors,...
    'axes', subplot(2,1,2),...
    'axis', 'tight',...
    ...'xticklabel', false,...
    'save', OPT.save,...
    'suffix', 'Velocities',...
    'title', 'Mean cross-shore velocity @ offshore boundary');


%%
% fh(end+1) = plot_engine(mat2cell(result.theta, 17, ones(size(result.theta,2),1)), mat2cell(result.ee_mean, 17, ones(size(result.theta,2),1)), 'ee_mean', 'J/m2/rad',...
%     'figure', figure,...
%     'axes', subplot(2,1,1),...
%     'axis', 'tight',...
%     'xticklabel', false,...
%     'title', 'sum ee mean for all transects',...
%     'colors', repmat('rgbkmc', 1, ceil(size(result.ee_mean,2) / 6)),...
%     'legend', false);
% fh(end+1) = plot_engine({result.theta}, {result.ee_mean}, 'ee mean', 'J/m2/rad',...
%     'figure', figure,...
%     'axes', subplot(1,1,1),...
%     'axis', 'tight',...
%     'title', 'mean sum ee mean for all transects',...
%     'xlabel', 'wave angle [degrees]',...
%     'save', OPT.save,...
%     'legend', legtxt);

%%
linkaxes(findobj(fh, 'type', 'axes', '-and', '-not', 'tag', 'legend'), 'x')

%%
function fh = plot_engine(x, y, varname, varunit, varargin)

OPT = struct(...
    'legend', false,...
    'displayname', {repmat({''}, size(x))},...
    'legendlocation', '',...
    'figure', true,...
    'axes', true,...
    'axis', false,...
    'plotfunction', @plot,...
    'title', true,...
    'save', false,...
    'style', @stylePresentation,...
    'xlim', true,...
    'xtick', true,...
    'xticklabel', true,...
    'xlabel', false,...
    'ylim', true,...
    'ytick', true,...
    'ylabelFcn', @(varname, varunit) sprintf('%s [%s]', varname, varunit),...
    'suffix', strrep(varname, ' ', '_'),...
    'colors', 'rgbkmc');

OPT = setproperty(OPT, varargin);

if islogical(OPT.figure) && OPT.figure
    fh = figure;
elseif ishandle(OPT.figure)
    fh = OPT.figure;
else
    fh = gcf;
end

if islogical(OPT.axes) && OPT.axes
    ah = axes('parent', fh);
elseif ishandle(OPT.axes)
    ah = OPT.axes;
else
    ah = gca;
end

if ischar(OPT.colors)
    OPT.colors = regexp('rgbkmc', '\D', 'match');
end

for i = 1:length(x)
    ph(i) = feval(OPT.plotfunction, x{i}, y{i},...
        'displayname', OPT.displayname{i},...
        'color', OPT.colors{i},...
        'parent', ah);
    
    if i == 1 && length(x) > 1
        hold on
    end
end

if ~isempty(OPT.ylabelFcn)
    ylabel(OPT.ylabelFcn(varname, varunit))
end

if ischar(OPT.xlabel)
    xlabel(OPT.xlabel)
end

if isnumeric(OPT.xlim)
    set(ah, 'xlim', OPT.xlim)
elseif islogical(OPT.xlim) && ~OPT.xlim
    set(ah, 'xlim', [])
end

if isnumeric(OPT.xtick)
    set(ah, 'xtick', OPT.xtick)
elseif islogical(OPT.xtick) && ~OPT.xtick
    set(ah, 'xtick', [])
end

if isnumeric(OPT.xticklabel)
    set(ah, 'xticklabel', OPT.xticklabel)
elseif islogical(OPT.xticklabel) && ~OPT.xticklabel
    set(ah, 'xticklabel', [])
end

if isnumeric(OPT.ytick)
    set(ah, 'ytick', OPT.ytick)
elseif islogical(OPT.ytick) && ~OPT.ytick
    set(ah, 'ytick', [])
end

if islogical(OPT.legend) && OPT.legend
    lh = legend('show');
elseif ~islogical(OPT.legend)
    lh = legend(OPT.legend);
end

if exist('lh', 'var') && ~isempty(OPT.legendlocation)
    set(lh,...
        'location', OPT.legendlocation)
end    

if ischar(OPT.title)
    title(OPT.title)
end

if OPT.axis
    axis(OPT.axis);
end

if OPT.save
    figureStyle(OPT.style);
    [fpath fname] = fileparts(which(mfilename));
    fname = sprintf('%s_%s', fname, strrep(OPT.suffix, ' ', '_'));
    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

