0
votes

enter image description hereIt is well documented that to find the angle of intersection between two planes we use the dot product of the normal vectors perpendicular to each plane - which then gives the cosine of the angle. I tried to program this but realized that in some cases the calculation results in the supplement rather than the angle of intersection itself. Here are two example codes below that shows this:

x = linspace(-10,20, 12)
y1 = (0*x) + 9
y2 = -x + 19

figure
k = plotyy(x,y1, x,y2);
set(k(2),'YDir','reverse')

%vectors and normal vectors 
n1 = [0, 1];
v1 = [1, 0];

n2 = [1, 1];
v2 = [1, -1];

angle = (acos(dot(v1, v2) / (norm(v1) * norm(v2))) * 180/pi)

the second:

x = linspace(-10,10, 12);
y1 = -(0.5*x) + 1.333;

plot(x, y1); hold on

%2nd line
xd = 5;
plot(xd, x, 'o')

%vectors and normal vectors 
n1 = [0.5, 1]; v1 = [1, -0.5];
n2 = [-5, 0]; v2 = [0, 5];

angle = (acos(dot(v1, v2) / (norm(v1) * norm(v2))) * 180/pi)

Note that the 1st example calculates the correct angle (45 deg), but the 2nd example calculates the supplement (116.5651 deg). After several attempts to crack this, I realized that if one normal points in a positive orientation and the other points in a negative orientation (see fig. A). Then: angle = 180 -(acos(dot(n1, n2) / (norm(v1) * norm(v2)))) * 180/pi. However, if both n1 and n2 we're to point in the same direction (either positive or negative) as in fig. B, then: angle = acos(dot(n1, n2) / (norm(v1) * norm(v2)))) * 180/pi

I have tested this with several examples, and i am sure the convention will work in all cases. I am also very sure that this will be useful to some many folks out there. Nevertheless, the bugging issue for me is how to program this. Any advice/help/suggestions will greatly be appreciated. Thanks!

1
Can you give the output versus the expected output?Eli Sadoff

1 Answers

1
votes

Two planes form two pairs of angles (a, Pi-a, a, Pi-a). All of them are correct, of course. And arrcosine approach gives one from these angles in range 0..Pi.

If you have oriented planes, defined by normal directions, you can calculate angle in range -Pi..Pi between normal directions using function:

Angle = atan2(vectorproduct(normal1, normal2), dot(normal1, normal2))

Note that sign of this angle depends on order of normal vectors!

Delphi example for your 2D data gives

angle -45.00
angle 116.57

just as I expected from paper sketch - method gives angle needed to rotate v1 to make it collinear with v2

var
  an, v1x, v1y, v2x, v2y: Double;
begin
  v1x := 1;
  v1y := 0;
  v2x := 1;
  v2y := -1;
  an := RadToDeg(ArcTan2(v1x * v2y - v2x * v1y, v1x * v2x + v1y * v2y));
  Memo1.Lines.Add(Format('angle %5.2f', [an]));
  v1x := 1;
  v1y := -0.5;
  v2x := 0;
  v2y := 5;
  an := RadToDeg(ArcTan2(v1x * v2y - v2x * v1y, v1x * v2x + v1y * v2y));
  Memo1.Lines.Add(Format('angle %5.2f', [an]));