function [X Y Z alfa propertyVar] = XBeach_GridOrientation(xyz,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(...
    'dx', 200,...
    'dy', 200);
    
% apply custom properties
OPT = setProperty(OPT, varargin{:});

%% setup regular grid
xmin = min(min(xyz{1}(:,1)),min(xyz{2}(:,1)));
xmax = max(max(xyz{1}(:,1)),max(xyz{2}(:,1)));
ymin = min(min(xyz{1}(:,2)),min(xyz{2}(:,2)));
ymax = max(max(xyz{1}(:,2)),max(xyz{2}(:,2)));

%% show data, for selection
m = ['sv+o.><xd'];
figure;
for i = 1:length(xyz)
    ind = ceil(length(xyz{i})/10000);
    scatter(xyz{i}(1:ind:end,1), xyz{i}(1:ind:end,2), 5, xyz{i}(1:ind:end,3),m(i)); hold on;
end
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;

OPT.xori = xi(1);
OPT.yori = yi(1);
OPT.xend_yori = [xi(2) yi(2)];
OPT.xend_yend = [xi(3) yi(3)];

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

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

xn = cos(alfa)*(OPT.xend_yori(1)-OPT.xori)+sin(alfa)*(OPT.xend_yori(2)-OPT.yori);
yn = -sin(alfa)*(OPT.xend_yend(1)-OPT.xori)+cos(alfa)*(OPT.xend_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);

for i = 1:length(xyz);
    xyzxb{i}(:,1) = A(1,1)*(xyz{i}(:,1)-OPT.xori)+A(1,2)*(xyz{i}(:,2)-OPT.yori);
    xyzxb{i}(:,2) = A(2,1)*(xyz{i}(:,1)-OPT.xori)+A(2,2)*(xyz{i}(:,2)-OPT.yori);
    xyzxb{i}(:,3) = xyz{i}(:,3);
    filename = ['xyz',num2str(i),'.xyz'];
    temp0 = xyzxb{i};
    save(filename,'temp0','-ascii');
end

xyzmerge = NaN*ones(size(X));
ind = 100;
for i = 1:length(xyzxb)
    temp{i} = griddata(xyzxb{i}(1:ind:end,1),xyzxb{i}(1:ind:end,2),xyzxb{i}(1:ind:end,3),X,Y);
    xyzmerge(~isnan(temp{i})) = temp{i}(~isnan(temp{i}));   
end

Z = xyzmerge;
Xbathy = A(1,1)*(X-OPT.xori)+A(1,2)*(Y-OPT.yori);
Ybathy = A(2,1)*(X-OPT.xori)+A(2,2)*(Y-OPT.yori);

% figure; 
% subplot(131); pcolor(X,Y,temp{1}); shading flat; 
% subplot(132); pcolor(X,Y,temp{2}); shading flat;
% subplot(133); pcolor(X,Y,xyzmerge); shading flat;

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