function [X Y Z alfa propertyVar] = XBeach_GridOrientation(xw, yw, zw, varargin)
%XBEACH_GRIDORIENTATION  One line description goes here.
%
%   More detailed description goes here.
%
%   Syntax:
%   varargout = XBeach_GridOrientation(varargin)
%
%   Input:
%   varargin  =
%
%   Output:
%   varargout =
%
%   Example
%   XBeach_GridOrientation
%
%   See also

%   --------------------------------------------------------------------
%   Copyright (C) 2009 Deltares
%       Dano Roelvink / Ap van Dongeren / C.(Kees) den Heijer
%
%       Kees.denHeijer@Deltares.nl	
%
%       Deltares
%       P.O. Box 177
%       2600 MH Delft
%       The Netherlands
%
%   This library is free software; you can redistribute it and/or
%   modify it under the terms of the GNU Lesser General Public
%   License as published by the Free Software Foundation; either
%   version 2.1 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
%   Lesser General Public License for more details.
%
%   You should have received a copy of the GNU Lesser General Public
%   License along with this library; if not, write to the Free Software
%   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307
%   USA
%   or http://www.gnu.org/licenses/licenses.html, http://www.gnu.org/, http://www.fsf.org/
%   --------------------------------------------------------------------

% Created: 03 Feb 2009
% Created with Matlab version: 7.6.0.324 (R2008a)

% $Id$
% $Date$
% $Author$
% $Revision$
% $HeadURL$

%% default properties
OPT = struct(...
    'manual', false,...
    'dx', 2,...
    'dy', 2,...
    'xori', 0,...
    'yori', 0,...
    'xend_y0', [max(max(xw)) 0],...
    'x_yend', [0 max(max(yw))]);

% apply custom properties
OPT = setProperty(OPT, varargin{:});
OPT.grid = 1;
if size(xw,2)>1
    OPT.grid = 2;
end
    
%% show data, for selection
if OPT.manual
    if OPT.grid==1
        figure;
        scatter(xw, yw, 5, zw, 'filled');
        axis([min(xw)-.5*(max(xw)-min(xw)) ...
            max(xw)+.5*(max(xw)-min(xw)) ...
            min(yw)-.5*(max(yw)-min(yw)) ...
            max(yw)+.5*(max(yw)-min(yw)) ])
        axis equal;
        axis manual
        colorbar
        hold on

        % Loop, picking up the points.
        disp('Click grid corner x=0,y=0')
        disp('Then click point x=xn,y=0')
        disp('Finally click to select extent of y')
        [xi yi] = select_oblique_rectangle;
    else
        figure;
        ind = ceil(size(xw)./100);
        xwt = xw(1:ind(1):end,1:ind(2):end); ywt = yw(1:ind(1):end,1:ind(2):end); zwt = zw(1:ind(1):end,1:ind(2):end);
        % limit number of points plotted
        
        scatter(xwt(:), ywt(:), 5, zwt(:),'filled');
        axis([min(xwt(:))-.5*(max(xwt(:))-min(xwt(:))) ...
            max(xwt(:))+.5*(max(xwt(:))-min(xwt(:))) ...
            min(ywt(:))-.5*(max(ywt(:))-min(ywt(:))) ...
            max(ywt(:))+.5*(max(ywt(:))-min(ywt(:))) ])
        axis equal;
        axis manual
        colorbar
        hold on

        % Loop, picking up the points.
        disp('Click grid corner x=0,y=0')
        disp('Then click point x=xn,y=0')
        disp('Finally click to select extent of y')
        [xi yi] = select_oblique_rectangle;
    end
    OPT.xori = xi(1);
    OPT.yori = yi(1);
    OPT.xend_y0 = [xi(2) yi(2)];
    OPT.x_yend = [xi(3) yi(3)];
end

alfa = atan2(OPT.xend_y0(2)-OPT.yori, OPT.xend_y0(1)-OPT.xori);


%
A = [cos(alfa) sin(alfa); -sin(alfa) cos(alfa)];
iA = inv(A);

xn = A(1,1)*(OPT.xend_y0(1)-OPT.xori)+A(1,2)*(OPT.xend_y0(2)-OPT.yori);
yn = A(2,1)*(OPT.x_yend(1)-OPT.xori)+A(2,2)*(OPT.x_yend(2)-OPT.yori);
xx = (0:OPT.dx:xn)';
yy = 0:OPT.dy:yn;
X = repmat(xx, 1, length(yy));
Y = repmat(yy, length(xx), 1);

if OPT.grid==1
    % in case grid =1
    Xt = A(1,1)*(xw-OPT.xori)+A(1,2)*(yw-OPT.yori);
    Yt = A(2,1)*(xw-OPT.xori)+A(1,1)*(yw-OPT.yori);
else
    % in case grid = 2 transform X and Y to coordinate system bathymetry data
    % for interpolation
    Xt = iA(1,1)*X + iA(1,2)*Y + OPT.xori;
    Yt = iA(2,1)*X + iA(2,2)*Y + OPT.yori;
end


if OPT.grid==2
    Z = interp2(xw, yw, zw, Xt, Yt);
    Xbathy = A(1,1)*(xwt-OPT.xori)+A(1,2)*(ywt-OPT.yori); 
    Ybathy = A(2,1)*(xwt-OPT.xori)+A(1,1)*(ywt-OPT.yori);
    Zbathy = zwt;
else
    try
        Z = griddata(Xt, Yt, zw, X, Y);
    catch
        Err = lasterror;
        if strcmp(Err.identifier, 'MATLAB:qhullmx:UndefinedError')
            Z = griddata(Xt,Yt,zw,X,Y,'linear',{'QJ'});
        else
            rethrow(Err)
        end
    end
    Xbathy = Xt; 
    Ybathy = Yt;
    Zbathy = zwt;
end
   
    

propertyVar = {...
    'dx', OPT.dx,...
    'dy', OPT.dy,...
    'xori', OPT.xori,...
    'yori', OPT.yori,...
    'bathy', {Xbathy Ybathy Zbathy},...
    'grid',OPT.grid};
