function [x_rotated, y_rotated]=rotatevector(x,y,theta);
% define the x- and y-data for the original line we would like to rotate
% x = 0:1; % 1:10;
% y = [0 0]; 1:10;
% create a matrix of these points, which will be useful in future calculations
v = [x;y];
% choose a point which will be the center of rotation
x_center = 0; % x(3);
y_center = 0; % y(3);
% create a matrix which will be used later in calculations
center = repmat([x_center; y_center], 1, length(x));
% define a 60 degree counter-clockwise rotation matrix
% theta = pi/3;       % pi/3 radians = 60 degrees
R = [cos(theta) -sin(theta); sin(theta) cos(theta)];
% do the rotation...
s = v - center;     % shift points in the plane so that the center of rotation is at the origin
so = R*s;           % apply the rotation about the origin
vo = so + center;   % shift again so the origin goes back to the desired center of rotation
% this can be done in one line as:
% vo = R*(v - center) + center
% pick out the vectors of rotated x- and y-data
x_rotated = vo(1,:);
y_rotated = vo(2,:);

% % make a plot
% plot(x, y, 'k-', x_rotated, y_rotated, 'r-', x_center, y_center, 'bo');
% axis equal
%%                 
%                 % Rotation matrix:
%                 theta=phi(n);
%                 R = [cos(theta) -sin(theta); sin(theta) cos(theta)];
%                 % if V=[u; v] (is a column vector), then R*V gives the rotated vector

%
% % From Patagonia paper
% % Converting winds to rotated coordinate along a fixed angle
%                 D=90-ones(81,81)*phi(n); % rotation angle from the U axis, measured positive counter-clockwise.
%                 RperD=pi/180; % to convert radians to degres
%                 
%                 Ustream(:,:,n)=U(:,:,n).*cos(D.*RperD) + V(:,:,n).*sin(D.*RperD);
%                 Vstream(:,:,n)=-U(:,:,n).*sin(D.*RperD) + V(:,:,n).*cos(D.*RperD);
%                 % As expected, if D=0 then Ustream = U, and Vstream = V.
%                 % Rotate coordinate system, relative to mean wind?