Angle between two vectors

Question

So, I'm trying to get the angle between two TPoints in Delphi, and it turns out to be harder then what I expected. The result I'm getting I can't explain (seems to be some problem with "to degrees"-part, or ArcTan2 does not return a sum in the form I expected. - Delpi-v7:

function Modulo(x,y:Extended): Extended;
var d: Extended;
begin
  d := x / y;
  Result := (d - floor(d)) * y;
end;

function Degrees(Rads: Extended): Extended;
begin
  Result := Rads*(180/Pi);
end;

function GetPointAngle(P1, P2: TPoint): Extended;
begin
  Result := Modulo(Degrees(ArcTan2(-(P1.Y - P2.Y), P1.X - P2.X)) - 90, 360);
end;

Yet, when I port the code to Python, or test it in another Pascal-variant, the above works. But now, it seems to return a sum that's static (not changing if I "move" the second TPoint).

In case your wondering; I created "modulo"-function simply because the divide-operator used in the "mod"-operator rounds to 0, and not down (so negative numbers don't work).

Edit: I noted that the value (angle) returned from GetPointAngle() increases when p gets further away from the other point c (and vice versa), even tho the TPoint (p) is dragged along the X-axis of the second TPoint (c).

EDIT:

You guys have outdone your self, I've looked over most of the answers, and it seems to be hard to choose best answer! And since you guys wrote everything with such detail, I will go trough everything with the same detail :-)

Also: what I did not share in my initial post, is that my function is being exported as a DLL to be reached from another pascal-interpretor (which is delphi-compatible).

Solution at last (changed):

GetPointAngle(P1, P2: TPoint) To: GetPointAngle(const P1, P2: TPoint)

^ I don't understand the need of declaring constants...

Answer 1

I assume you want to calculate the angle relative to the X-axis of the line which is formed between those two points.

For this situation, the following formula applies:

Tan(a) = (P2.Y - P1.Y) / (P2.X - P1.X)

Which translates to:

a = ArcTan((P2.Y - P1.Y) / (P2.X - P1.X))

When the two points have the same X coordinate, this will obviously result in a EDivByZero exception, so you have to take care of that yourself. Furthermore, ArcTan results in an angle within the range 0°..90° (i.e. 0..π/2) and thus disregards the correct quadrant, while ArcTan2 results in an angle within -180°..180°. Add 360° to the result to convert a negative angle to positive:

function AngleOfLine(const P1, P2: TPoint): Double;
begin
  if P2.X = P1.X then
    if P2.Y > P1.Y then
      Result := 90
    else
      Result := 270
  else
    Result := RadToDeg(ArcTan2(P2.Y - P1.Y, P2.X - P1.X));
  if Result < 0 then
    Result := Result + 360;
end;

Which results in:

  A := AngleOfLine(Point(10, 10), Point(20, 10)); // 0
  A := AngleOfLine(Point(10, 10), Point(20, 20)); // 45
  A := AngleOfLine(Point(10, 10), Point(10, 20)); // 90
  A := AngleOfLine(Point(10, 10), Point(0, 20));  // 135
  A := AngleOfLine(Point(10, 10), Point(0, 10));  // 180
  A := AngleOfLine(Point(10, 10), Point(0, 0));   // 225
  A := AngleOfLine(Point(10, 10), Point(10, 0));  // 270
  A := AngleOfLine(Point(10, 10), Point(20, 0));  // 315

Now, this is relative to the world coordinate system which has its positive Y-axis pointed upwards by default. If you want to convert the result to the device coordinate system wherein the positive Y-axis points downwards, then subtract the result from 360°:

Result := 360 - Result;

Update:

It seems ArcTan2 dóes take care of division by zero, (even in D7 inspite of the documentation) so the routine becomes much simpler:

function AngleOfLine(const P1, P2: TPoint): Double;
begin
  Result := RadToDeg(ArcTan2((P2.Y - P1.Y),(P2.X - P1.X)));
  if Result < 0 then
    Result := Result + 360;
end;

Edit:

I noted that the value returned from GetPointAngle() increases when p gets furter away from the other point c (and vice versa).

That depends. Looking at the diagram above, if the second point moves further along the x-axis, the angle decreases. If the second point moves further along the y-axis, the angle increases. Of course, this depends on which quadrant both points are in.

Furthermore, your code negates the first parameter of ArcTan2 and subtracts another 90° from the result. I do not know what you mean by that and whether it is intentional, but it could be the source of unexpected results.