function varargout = decc_plot(varargin)
%DECC_PLOT  One line description goes here.
%
%   More detailed description goes here.
%
%   Syntax:
%   varargout = decc_plot(varargin)
%
%   Input:
%   varargin  =
%
%   Output:
%   varargout =
%
%   Example
%   decc_plot
%
%   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: 18 Jan 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', 'run(40).*1900.*xb_sum.nc',...
    'radiusfun', @(ncfile) str2double(regexp(ncfile, '(?<=radius)\d{5}',  'match')),...
    'phi_limits', [-90 90],...
    'save', true);

OPT = setproperty(OPT, varargin);

%%

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

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

ncfile = files{fid};

%% load data
x = nc_varget(ncfile, 'globalx');
y = nc_varget(ncfile, 'globaly');
g = xb_stagger(x,y);

phi = -rad2deg(mean(g.alfaz,1));
yid = phi > min(OPT.phi_limits) & phi < max(OPT.phi_limits);
x = x(yid,:);
y = y(yid,:);
g = structfun(@(x) x(:,yid), g, 'uniformoutput', false);

xbo.ue_mean = nc_varget(ncfile, 'ue_mean');
xbo.ue_mean = xbo.ue_mean(yid,:);
xbo.ve_mean = nc_varget(ncfile, 've_mean');
xbo.ve_mean = xbo.ve_mean(yid,:);
% Svsg = nc_varget(ncfile, 'Svsg_mean');
% dSvsg = diff(Svsg, 1, 2);
% dSvsg(:,end+1,:) = dSvsg(:,end,:);
% zs = nc_varget(ncfile, 'zs_mean');
xbo.zb = nc_varget(ncfile, 'zb');
idx = repmat({':'}, 1, ndims(xbo.zb));
idx{size(xbo.zb) == length(yid)} = yid;
xbo.zb = xbo.zb(idx{:});

% xb = xb_read_output(ncfile, 'vars', {'DIMS' 'zb' 'zs_mean' 'ue_mean' 've_mean'});
% xbo = xs_peel(xb);

%% x and y velocities
decc_plot_engine(g.xu, g.yu, xbo.ue_mean',...
    'Cross-shore velocity', 'm/s', OPT.radiusfun(ncfile),...
    'figure', figure,...
    'title', 'Cross-shore velocity',...
    'xtick', false,...
    'ytick', false,...
    'axes', subplot(1,2,1))
decc_plot_engine(g.xv, g.yv, xbo.ve_mean',...
    'Alongshore velocity', 'm/s', OPT.radiusfun(ncfile),...
    'figure', gcf,...
    'axes', subplot(1,2,2),...
    'title', 'Alongshore velocity',...
    'save', OPT.save,...
    'xtick', false,...
    'ytick', false,...
    'suffix', 'velocities')

decc_plot_engine(g.xz, g.yz, squeeze(xbo.zb(1,:,:))',...
    'Bed level', 'm', OPT.radiusfun(ncfile),...
    'plotfunction', @surf,...
    'suffix', 'Grid and Bathymetry',...
    'colorbar', false,...
    'title', 'Bathymetry',...
    'figure', figure,...
    'axes', subplot(4,2,[3 5]),...
    'axis', 'tight',...
    'xtick', false,...
    'ytick', false,...
    'ztick', false,...
    'shading', 'interp',...
    'save', false)
view(-36.5, 72)

xthinning = 5;
ythinning = 5;
% make sure first and last index are included and distribute remaining
% points "as evenly as possible"
[ny nx] = size(x);
xtid = round(linspace(1, nx, nx/xthinning));
ytid = round(linspace(1, ny, ny/ythinning)); 

decc_plot_engine(g.xz(xtid,ytid), g.yz(xtid,ytid), squeeze(xbo.zb(1,ytid,xtid))',...
    'Bed level', 'm', OPT.radiusfun(ncfile),...
    'suffix', 'Grid and Bathymetry',...
    'figure', gcf,...
    'axes', subplot(4,2,2:2:8),...
    'title', 'Curvilinear grid',...
    'caxis', false,...
    'xtick', false,...
    'ytick', false,...
    'shading', false,...
    'save', OPT.save)
% decc_plot_engine(g.xv, g.yv, squeeze(Svsg(end,:,:))',...
%     'Alongshore transport', 'm^3/s', OPT.radius)
% decc_plot_engine(g.xv, g.yv, squeeze(dSvsg(end,:,:))'./g.dnz,...
%     'Alongshore transport gradient', 'm^3/s/m', OPT.radius)
% decc_plot_engine(g.xv, g.yv, squeeze(mean(zs,1))',...
%     'Water level', 'm', OPT.radius, 'caxis', [4.5 5.5])

function decc_plot_engine(x, y, c, varname, varunit, radius, varargin)

OPT = struct(...
    'caxis', true,...
    'figure', true,...
    'axes', true,...
    'axis', 'image',...
    'plotfunction', @pcolor,...
    'colorbar', true,...
    'title', true,...
    'save', false,...
    'xtick', true,...
    'ytick', true,...
    'ztick', true,...
    'suffix', strrep(varname, ' ', '_'),...
    'shading', 'flat');
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
    
ph = feval(OPT.plotfunction, x, y, c,...
    'parent', ah);

if OPT.shading
    shading(OPT.shading);
end

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

if OPT.colorbar
    cb = colorbar;
    if islogical(OPT.caxis) && OPT.caxis
        climits = caxis;
        climits = [-1 1] * max(abs(climits));
        caxis(climits)
    elseif isnumeric(OPT.caxis)
        caxis(OPT.caxis)
    end
    set(get(cb, 'YLabel'), 'string', sprintf('%s [%s]', varname, varunit))
end

if islogical(OPT.title) && OPT.title
    title(sprintf('%s; radius = %i m', varname, radius))
elseif ischar(OPT.title)
    title(OPT.title)
end

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

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

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

if OPT.save
    figureStyle(@stylePresentation,...
        {'figure', 'Units', 'Points'});%, 'Position', [200 200 350 250]});
    [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