function theta = findAngle2d(x1, y1, x2, y2)

% findangle(x1, y1, x2, y2) -- Find the angle between two points, in degrees.
%
% Given two points, with coordinates (x1, y1) and (x2, y2), the 
% findangle() function will compute the angle in degrees.
% This function uses the trigonometric identity: 
% COS(THETA) = A / C
% where A and B are the lengths of the short sides of the right
% triangle for which the line [(x1, y1), (x2, y2)] is the hypotenuse.
% Then we use
% ARCCOS(COS(THETA)) = THETA
% to find the angle, theta. Note that theta is reported as degrees
% counterclockwise, with due east = 0 (thus, North is 90, West
% is 180, and south is 270).

  a = x2 - x1;             % length of side a (the x-axis side)
  b = y2 - y1;             % length of side b (the y-axis side)
  c = sqrt(a^2 + b^2);     % length of the hypotenuse


  CosTheta = a / c;             % Now we know the tangent of theta

  theta = acos(CosTheta);       % Take the inverse of cos(theta) 

% Figure out which quadrant point 2 is in, relative to point 1

  if ((x2 > x1) & (y2 > y1))
    % point 2 is in Quadrant I (0-90 deg. range)
    theta = rad2deg(theta);          

  elseif ((x2 > x1) & (y2 < y1))
    % point 2 is in Quadrant II (270-360 range)
    theta = 360 - rad2deg(theta);

  elseif ((x2 < x1) & (y2 < y1))
    % point 2 is in Quadrant III (270-360 range)
    theta = 360 - rad2deg(theta);

  else
    % point 2 is in Quadrant IV (90-180 range)
    theta = rad2deg(theta);
  end