function [x y z phi] = decc_curved_grid(varargin)
%DECC_MODEL_SETUP  XBeach model setup for alongshore uniform curved coastline model
%
%   More detailed description goes here.
%
%   Syntax:
%   varargout = decc_model_setup(varargin)
%
%   Input:
%   varargin  =
%
%   Output:
%   varargout =
%
%   Example
%   decc_model_setup
%
%   See also

%% Copyright notice
%   --------------------------------------------------------------------
%   Copyright (C) 2011 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: 10 Nov 2011
% Created with Matlab version: 7.7.0.471 (R2008b)

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

%%
OPT = struct(...
    'x', [0 2549.97 2724.88 2774.97 2804.97 3030.6]',...
    'y', [0 400*pi],...
    'z', [-20 -3 0 3 15 15]',...
    'contour', 5,...
    'radius', 500,...
    'dymin', 20,...
    'vardy', false,...
    'phi_unit', 'rad',...
    'phi_mean', pi,...
    'curvi', true ...
    );

OPT = setproperty(OPT, varargin{:});

%%
x_contour = findCrossings(OPT.x, OPT.z, minmax(OPT.x')', ones(2,1)*OPT.contour);

% distance to center of circle (amplitude)
a = OPT.radius + x_contour - OPT.x(:);

if OPT.curvi
    
    % express phi in degrees
    if strcmpi(OPT.phi_unit, 'deg')
        phi_mean = deg2rad(OPT.phi_mean);
    else
        phi_mean = OPT.phi_mean;
    end
    
    Y = 2*pi*OPT.radius;
    phi_range = diff(OPT.y)/Y * 2 * pi * [-.5 .5] + phi_mean;
    
    % derive the phi-grid
    if OPT.vardy
        y_length = diff(phi_range) * OPT.radius;
        ygr = xb_grid_ygrid([0 y_length],...
            'dymin', OPT.dymin);
        phi = ygr/max(ygr) * diff(phi_range) + phi_range(1);
    else
        d_phi = atan(OPT.dymin/OPT.radius);
        ny = max( [ceil(abs(diff(phi_range)/d_phi)) 3] );
        if mod(ny,2) == 0
            % make sure ny is odd, in order to truely have the middle transect
            % available
            ny = ny + 1;
        end
        phi = linspace(phi_range(1), phi_range(2), ny);
    end
    
    phi = sort(phi, 'descend');
    
    % calculate x, y and z
    x = a*cos(phi);
    y = sign(OPT.radius)*a*sin(phi);
    
    % lateral extend
    % nextra = 5;
    % dy1 = repmat(linspace(-OPT.dymin*nextra, -OPT.dymin, nextra)*cos(phi(1)), size(x,1), 1);
    % dx1 = repmat(linspace(-OPT.dymin*nextra, -OPT.dymin, nextra)*sin(phi(1)), size(x,1), 1);
    % dy2 = repmat(linspace(-OPT.dymin*nextra, -OPT.dymin, nextra)*cos(phi(end)), size(x,1), 1);
    % dx2 = repmat(linspace(-OPT.dymin*nextra, -OPT.dymin, nextra)*sin(phi(end)), size(x,1), 1);
    % x = [repmat(x(:,1), 1, nextra) + dx1 x repmat(x(:,end), 1, nextra) - dx2];
    % y = [repmat(y(:,1), 1, nextra) - dy1 y repmat(y(:,end), 1, nextra) + dy2];
    
    z = repmat(OPT.z(:), 1, size(x,2));
    
    [x y z] = xb_grid_finalise(x,y,z,'actions',{'lateral_extend_curvi'});
    
else
    % rectalinear grid
    dymin = min(diff(OPT.x));
    radius = OPT.radius +  dymin - mod((max(a) - max(OPT.x)), dymin);
    a = radius + x_contour - OPT.x(:);
    x  = [OPT.x(:); (max(OPT.x) + dymin:dymin:max(a))'];
    y = unique([x; sort(-x) + 2*max(x)]);
    [x y] = meshgrid(x, y);
    x0 = x(mean(1:size(x,1)), end);
    y0 = y(mean(1:size(y,1)), end);
    dx = x - x0;
    dy = y - y0;
    r = sqrt(dy.^2 + dx.^2);
    warning('rectalinear grid creation not finalised')
end


%%
figure;
subplot(211)
surf(x,y,z)
shading interp
subplot(212)
plot(x,y, '-o')
axis square