finding angle between two points

Join Bytes! wrote:[color=blue]
> Hi gang,
> I have a grid full of particles in my program, and I want to find an
> angle between two particles. I'm having trouble because it seems like
> the answer depends on whether or not the target particle is above or
> below, in front or behind the refernce particle.
>
> If I have a reference particle at (10,10), and another particle at
> (20,20), i'm currently finding the angle by:
>
> angle = atan((10-20)/(10-20)) = 45 degrees. When I draw this out, I can
> see that the (20,20) object is 45 degrees below the x-axis from the
> (10,10) object. Fine.
>
> But if I look at another particle and compare it's angle to (10,10),
> for example a particle located at (2,20), I get something kind of
> different:
>
> atan((10-20)/(10-2)) = -51.34 degrees
>
> When I take x and y components of these angles, things don't add up.
> Somtimes it gives me a negative value of x when it should be positive.
>
> I don't think I've been very clear here, but I'm hoping this is a
> common enough mistake that someone knows what I'm trying to say. I'm
> wondering also if there is a better way to find a angle between two
> points that might be relative to a certain axis, because it seems like
> my angle is relative to different axes at different times. Thanks guys[/color]

Sign errors of this nature are a common problem when first using any
trigonometry for computation. It's not C++-specific, but here's a few
pointers which, if you follow them, should lead you right:

1. Think carefully about quadrants and how they relate to the
trigonometric function you're using, and the signs involved.
2. Remember that an infinite number of angles have the same
trigonometric measures.
3. Familiarize yourself with the behavior of the trig function (atan,
in this case) you're using. Read the documentation. You should know a
priori what the sign will be given a particular input.

You've got to know your trig, and you've got to read the docs. The
rest is just fiddling bits.

Luke

* Join Bytes!:[color=blue]
>
> I have a grid full of particles in my program, and I want to find an
> angle between two particles. I'm having trouble because it seems like
> the answer depends on whether or not the target particle is above or
> below, in front or behind the refernce particle.
>
> If I have a reference particle at (10,10), and another particle at
> (20,20), i'm currently finding the angle by:[/color]

There's no such thing as angle between two points; if you mean the angle
between the vectors, then that's the exactly 0 degrees in this case.

At the mathematics side you can use the dot product, see e.g.
<url: http://mathworld.wolfram.com/DotProduct.html>.

C++ doesn't much support trigonometry, so you'll probably have to roll
your own arc cos function, based on what's there.

--
A: Because it messes up the order in which people normally read text.
Q: Why is it such a bad thing?
A: Top-posting.
Q: What is the most annoying thing on usenet and in e-mail?

On 1 Feb 2006 19:30:32 -0800, Join Bytes! wrote:
[color=blue]
>Hi gang,
>I have a grid full of particles in my program, and I want to find an
>angle between two particles. I'm having trouble because it seems like
>the answer depends on whether or not the target particle is above or
>below, in front or behind the refernce particle.
>
>If I have a reference particle at (10,10), and another particle at
>(20,20), i'm currently finding the angle by:
>
>angle = atan((10-20)/(10-20)) = 45 degrees. When I draw this out, I can
>see that the (20,20) object is 45 degrees below the x-axis from the
>(10,10) object. Fine.
>
>But if I look at another particle and compare it's angle to (10,10),
>for example a particle located at (2,20), I get something kind of
>different:
>
>atan((10-20)/(10-2)) = -51.34 degrees
>
>When I take x and y components of these angles, things don't add up.
>Somtimes it gives me a negative value of x when it should be positive.
>
>I don't think I've been very clear here, but I'm hoping this is a
>common enough mistake that someone knows what I'm trying to say. I'm
>wondering also if there is a better way to find a angle between two
>points that might be relative to a certain axis, because it seems like
>my angle is relative to different axes at different times. Thanks guys[/color]

Trig was a few years back.

Need to measure angles between lines, not points.
So must have a common reference.
The origin is typical: (0,0)
if in same quadrant >> Use: abs(arctan(angle #1) - arctan(angle #2))

arctan(20/2) - arctan(10/10)

remember >> rise over run

1.471128 - 0.785398 = 0.68573 radians

0.68573 radians x 180 degree / PI radians = 39.28941 degrees

10,10 is obviously 45 degrees
2,20 is close to the y-axis
answer slightly less than 45 degrees checks out.

if in quadrants #1 and #2 ... same as above
if in quadrants #1 and #4 ... draw the picture and add them correctly
your reference is the origin and the x-axis
if in quadrants #1 and #3 or #2 and #3 ... again .. draw the picture
look at each angle with the x-axis independently.
to add them or subtract them depends on the quadrants.
sometimes you have to subtract the calculated angle from 180 degrees
to get the working angle with the axis.
answer is typically 180 degrees or less .. unless direction to
traverse is specified. For example ... (10,1) & (10,0) could wrap
around and be near 360 degrees ... or could be very small ... depends
on problem specification.
good luck.

-- Tom

<cwsullivan@ucdavis.edu> wrote in message
news:1138851032.398293.4650@g47g2000cwa.googlegrou ps.com...[color=blue]
> Hi gang,
> I have a grid full of particles in my program, and I want to find an
> angle between two particles. I'm having trouble because it seems like
> the answer depends on whether or not the target particle is above or
> below, in front or behind the refernce particle.
>
> If I have a reference particle at (10,10), and another particle at
> (20,20), i'm currently finding the angle by:
>
> angle = atan((10-20)/(10-20)) = 45 degrees. When I draw this out, I can
> see that the (20,20) object is 45 degrees below the x-axis from the
> (10,10) object. Fine.
>
> But if I look at another particle and compare it's angle to (10,10),
> for example a particle located at (2,20), I get something kind of
> different:
>
> atan((10-20)/(10-2)) = -51.34 degrees
>
> When I take x and y components of these angles, things don't add up.
> Somtimes it gives me a negative value of x when it should be positive.
>
> I don't think I've been very clear here, but I'm hoping this is a
> common enough mistake that someone knows what I'm trying to say. I'm
> wondering also if there is a better way to find a angle between two
> points that might be relative to a certain axis, because it seems like
> my angle is relative to different axes at different times. Thanks guys[/color]

Yes, it is a bit of a hassle. Luckily I've needed this myself so provide
the code for you. The value is returned in radians instead of degrees
though, but you can modify it to return degrees if you want.

float CalcTheta( const JVEC2 Point1, const JVEC2 Point2 )
{
float Theta;
if ( Point2.x - Point1.x == 0 )
if ( Point2.y > Point1.y )
Theta = 0;
else
Theta = static_cast<float>( PI );
else
{
Theta = std::atan( (Point2.y - Point1.y) / (Point2.x - Point1.x) );
if ( Point2.x > Point1.x )
Theta = static_cast<float>( PI ) / 2.0f - Theta;
else
Theta = static_cast<float>( PI ) * 1.5f - Theta;
};
return Theta;
}

JVEC2 is simply a structure with a float x and a float y. You can make it
POD or as complex as you want.

Join Bytes! wrote:

[ ... ]
[color=blue]
> But if I look at another particle and compare it's angle to (10,10),
> for example a particle located at (2,20), I get something kind of
> different:
>
> atan((10-20)/(10-2)) = -51.34 degrees[/color]

Depending on your viewpoint, this angle is either -51.34, or 128.66
degrees. Unfortunately, when you divide the two numbers, their signs
(can) cancel out, so atan can only give answers in the first or fourth
quadrant (i.e. -90 to +90 degrees). As far as it cares, all positive
slopes are in the first quadrant, and all negative slopes are in the
fourth quadrant.

In most cases, atan2 works quite a bit better. You pass it the rise and
run separately, and it can take the signs of both into account to
determine the quadrant for the result. Keep in mind, however, that to
make real use of this, you need to pass them correctly -- for the
points above (for example):

atan2(10-20, 10-2) = -51.34 degrees
atan2(20-10, 2-10) = 128.66 degrees

Think of the first point you supply as a pivot, and the second as a
pointer. It's going to tell you the angle from the pivot to the
pointer, not vice versa.

--
Later,
Jerry.

Join Bytes! wrote:[color=blue]
> Hi gang,
> I have a grid full of particles in my program, and I want to find an
> angle between two particles. I'm having trouble because it seems like
> the answer depends on whether or not the target particle is above or
> below, in front or behind the refernce particle.
>
> If I have a reference particle at (10,10), and another particle at
> (20,20), i'm currently finding the angle by:
>
> angle = atan((10-20)/(10-20)) = 45 degrees. When I draw this out, I can
> see that the (20,20) object is 45 degrees below the x-axis from the
> (10,10) object. Fine.
>
> But if I look at another particle and compare it's angle to (10,10),
> for example a particle located at (2,20), I get something kind of
> different:
>
> atan((10-20)/(10-2)) = -51.34 degrees
>
> When I take x and y components of these angles, things don't add up.
> Somtimes it gives me a negative value of x when it should be positive.
>
> I don't think I've been very clear here, but I'm hoping this is a
> common enough mistake that someone knows what I'm trying to say. I'm
> wondering also if there is a better way to find a angle between two
> points that might be relative to a certain axis, because it seems like
> my angle is relative to different axes at different times. Thanks guys
>[/color]

Look up atan2

Alf P. Steinbach wrote:
[color=blue]
>
> C++ doesn't much support trigonometry, so you'll probably have to roll
> your own arc cos function, based on what's there.
>[/color]

Or just use acos. <g>

--

Pete Becker
Dinkumware, Ltd. (http://www.dinkumware.com)

Join Bytes! wrote:[color=blue]
> Hi gang,
> I have a grid full of particles in my program, and I want to find an
> angle between two particles.[/color]

What you proceed to describe is not "finding the angle between two
particles" but finding the angle between the vectors a and b:

a = p1 - p2
b = (1, 0)

where p1 and p2 are the particles.

You can calculate this angle simply using the dot product equation.

Recall:

a.b = |a||b|cos(theta)

where theta is the angle between the two vectors a and b. Rearrange this
equation to get:

theta = arccos((a.b)/(|a||b|))

As a c function:

double calcAngle(double p1_x, double p1_y, double p2_x, double p2_y)
{
double a_x = p2_x - p1_x;
double a_y = p2_y - p1_y;

double b_x = 1.0;
double b_y = 0.0;

return acos((a_x*b_x+a_y*b_y)/sqrt(a_x*a_x+a_y*a_y));
}

The only special case you have to worry about is when p1 == p2. The
angle is clearly undefined in this case and this will manifest itself in
a division by zero error when you try to divide by the magnitude. Apart
from that this function should work for all cases with no messy sign
checking.

--
Ben Radford
"Why is it drug addicts and computer aficionados are both called users?"

Ben Radford wrote:[color=blue]
> a = p1 - p2[/color]

Or 'a = p2 - p1', which is what the program calculates. It doesn't
really matter which way wround you do it, it will just switch your view
point in the same way Jerry said.

--
Ben Radford
"Why is it drug addicts and computer aficionados are both called users?"

In message <drt27k$56c$1@news.ox.ac.uk>, Ben Radford
<benjamin.radford@new.ox.ac.uk> writes[color=blue]
>cwsullivan@ucdavis.edu wrote:[color=green]
>> Hi gang,
>> I have a grid full of particles in my program, and I want to find an
>> angle between two particles.[/color]
>
>What you proceed to describe is not "finding the angle between two
>particles" but finding the angle between the vectors a and b:
>
>a = p1 - p2
>b = (1, 0)
>
>where p1 and p2 are the particles.
>
>You can calculate this angle simply using the dot product equation.
>
>Recall:
>
>a.b = |a||b|cos(theta)
>
>where theta is the angle between the two vectors a and b. Rearrange
>this equation to get:
>
>theta = arccos((a.b)/(|a||b|))
>
>As a c function:
>
>double calcAngle(double p1_x, double p1_y, double p2_x, double p2_y)
>{
> double a_x = p2_x - p1_x;
> double a_y = p2_y - p1_y;
>
> double b_x = 1.0;
> double b_y = 0.0;
>
> return acos((a_x*b_x+a_y*b_y)/sqrt(a_x*a_x+a_y*a_y));
>}
>
>The only special case you have to worry about is when p1 == p2. The
>angle is clearly undefined in this case and this will manifest itself
>in a division by zero error when you try to divide by the magnitude.
>Apart from that this function should work for all cases with no messy
>sign checking.[/color]

It won't work if the square of a_x or a_y can't be represented as a
double. There's a simple trick for evaluating the hypotenuse without
unnecessary overflow, but there's no call for such complications here.

What is wrong with using the standard function std::atan2(), as Jerry
Coffin suggested? In a single function call it takes care of all four
quadrants with no sign checking. And it really will work on all
combinations of inputs except the indeterminate case (0, 0).

--
Richard Herring

Richard Herring wrote:[color=blue]
> In message <drt27k$56c$1@news.ox.ac.uk>, Ben Radford
> <benjamin.radford@new.ox.ac.uk> writes
>[color=green]
>> Join Bytes! wrote:
>>[color=darkred]
>>> Hi gang,
>>> I have a grid full of particles in my program, and I want to find an
>>> angle between two particles.[/color]
>>
>>
>> What you proceed to describe is not "finding the angle between two
>> particles" but finding the angle between the vectors a and b:
>>
>> a = p1 - p2
>> b = (1, 0)
>>
>> where p1 and p2 are the particles.
>>
>> You can calculate this angle simply using the dot product equation.
>>
>> Recall:
>>
>> a.b = |a||b|cos(theta)
>>
>> where theta is the angle between the two vectors a and b. Rearrange
>> this equation to get:
>>
>> theta = arccos((a.b)/(|a||b|))
>>
>> As a c function:
>>
>> double calcAngle(double p1_x, double p1_y, double p2_x, double p2_y)
>> {
>> double a_x = p2_x - p1_x;
>> double a_y = p2_y - p1_y;
>>
>> double b_x = 1.0;
>> double b_y = 0.0;
>>
>> return acos((a_x*b_x+a_y*b_y)/sqrt(a_x*a_x+a_y*a_y));
>> }
>>
>> The only special case you have to worry about is when p1 == p2. The
>> angle is clearly undefined in this case and this will manifest itself
>> in a division by zero error when you try to divide by the magnitude.
>> Apart from that this function should work for all cases with no messy
>> sign checking.[/color]
>
>
> It won't work if the square of a_x or a_y can't be represented as a
> double. There's a simple trick for evaluating the hypotenuse without
> unnecessary overflow, but there's no call for such complications here.[/color]

Exactly. There are some changes that can be made to improve this
function but it was my intention to give something that matched the
original equation as closely as possible.
[color=blue]
> What is wrong with using the standard function std::atan2(), as Jerry
> Coffin suggested? In a single function call it takes care of all four
> quadrants with no sign checking. And it really will work on all
> combinations of inputs except the indeterminate case (0, 0).[/color]

The only case my function won't work on is when the square causes an
overflow as you pointed out. This is a limitation of the architecture
that can be worked around not a problem with the algorithm.

However there is nothing wrong with what Jerry suggested, I just didn't
examine his post clearly enough to realise exactly what he was saying.
The only difference is that atan2 gives an angle in the range -PI/2 <
theta < PI/2 whereas I believe the dot product method always gives a
positive angle. Got to run now or I'm going to be late for a lecture.

--
Ben Radford
"Why is it drug addicts and computer aficionados are both called users?"

In message <drt6m7$6pm$3@news.ox.ac.uk>, Ben Radford
<benjamin.radford@new.ox.ac.uk> writes[color=blue]
>Richard Herring wrote:[color=green]
>> In message <drt27k$56c$1@news.ox.ac.uk>, Ben Radford
>><benjamin.radford@new.ox.ac.uk> writes
>>[color=darkred]
>>> Join Bytes! wrote:
>>>
>>>> Hi gang,
>>>> I have a grid full of particles in my program, and I want to find an
>>>> angle between two particles.
>>>
>>>
>>> What you proceed to describe is not "finding the angle between two
>>>particles" but finding the angle between the vectors a and b:
>>>
>>> a = p1 - p2
>>> b = (1, 0)
>>>
>>> where p1 and p2 are the particles.
>>>
>>> You can calculate this angle simply using the dot product equation.
>>>
>>> Recall:
>>>
>>> a.b = |a||b|cos(theta)
>>>
>>> where theta is the angle between the two vectors a and b. Rearrange
>>>this equation to get:
>>>
>>> theta = arccos((a.b)/(|a||b|))
>>>
>>> As a c function:
>>>
>>> double calcAngle(double p1_x, double p1_y, double p2_x, double p2_y)
>>> {
>>> double a_x = p2_x - p1_x;
>>> double a_y = p2_y - p1_y;
>>>
>>> double b_x = 1.0;
>>> double b_y = 0.0;
>>>
>>> return acos((a_x*b_x+a_y*b_y)/sqrt(a_x*a_x+a_y*a_y));
>>> }
>>>
>>> The only special case you have to worry about is when p1 == p2. The
>>>angle is clearly undefined in this case and this will manifest itself
>>>division by zero error when you try to divide by the magnitude.
>>>Apart from that this function should work for all cases with no messy
>>>checking.[/color]
>> It won't work if the square of a_x or a_y can't be represented as a
>>double. There's a simple trick for evaluating the hypotenuse without
>>unnecessary overflow, but there's no call for such complications here.[/color]
>
>Exactly. There are some changes that can be made to improve this
>function but it was my intention to give something that matched the
>original equation as closely as possible.[/color]

That's _your_ equation. The OP just wanted the angle between two line
segments.[color=blue]
>[color=green]
>> What is wrong with using the standard function std::atan2(), as Jerry
>>Coffin suggested? In a single function call it takes care of all four
>>quadrants with no sign checking. And it really will work on all
>>combinations of inputs except the indeterminate case (0, 0).[/color]
>
>The only case my function won't work on is when the square causes an
>overflow as you pointed out.[/color]

That's a whole set of cases more than is necessary.
[color=blue]
>This is a limitation of the architecture that can be worked around not
>a problem with the algorithm.[/color]

If the algorithm design doesn't take proper account of domains and
bounds then it _is_ a problem with the algorithm.[color=blue]
>
>However there is nothing wrong with what Jerry suggested, I just didn't
>examine his post clearly enough to realise exactly what he was saying.
>The only difference is that atan2 gives an angle in the range -PI/2 <
>theta < PI/2 whereas I believe the dot product method always gives a
>positive angle.[/color]

I don't have to believe it, I know it. If a positive angle is what you
need, std::abs() is your friend.
[color=blue]
>Got to run now or I'm going to be late for a lecture.
>[/color]
Numerical Analysis 101?

--
Richard Herring

Richard Herring wrote:
<snip>
[color=blue][color=green][color=darkred]
>>> What is wrong with using the standard function std::atan2(), as Jerry
>>> Coffin suggested? In a single function call it takes care of all four
>>> quadrants with no sign checking. And it really will work on all
>>> combinations of inputs except the indeterminate case (0, 0).[/color]
>>
>>
>> The only case my function won't work on is when the square causes an
>> overflow as you pointed out.[/color]
>
>
> That's a whole set of cases more than is necessary.
>[color=green]
>> This is a limitation of the architecture that can be worked around not
>> a problem with the algorithm.[/color]
>
>
> If the algorithm design doesn't take proper account of domains and
> bounds then it _is_ a problem with the algorithm.[/color]

Indeed, but only when those bounds apply in the general case and are not
simply a limitation of the architecture. The bounds vary depending on
the platform and are therefore an implementation and not algorithmic
detail. By your argument an algorithm that is implemented differently on
two different platforms would actually be a different algorithm in each
case because it would have to take proper account of the domains and
bounds which apply for that case.
[color=blue][color=green]
>> However there is nothing wrong with what Jerry suggested, I just
>> didn't examine his post clearly enough to realise exactly what he was
>> saying. The only difference is that atan2 gives an angle in the range
>> -PI/2 < theta < PI/2 whereas I believe the dot product method always
>> gives a positive angle.[/color]
>
>
> I don't have to believe it, I know it.
> If a positive angle is what you
> need, std::abs() is your friend.[/color]

I don't see where you are going with this. I acknowledged that Jerry's
method was a valid way of calculating the answer the OP wanted. The
algorithm I gave was just an alternative; albeit with an overflow
problem in the hastily written implementation code - which you pointed out.
[color=blue][color=green]
>> Got to run now or I'm going to be late for a lecture.
>>[/color]
> Numerical Analysis 101?[/color]

Nice insult, but totally opportunistic and unnecessary in conveying your
point.

--
Ben Radford
"Why is it drug addicts and computer aficionados are both called users?"

In message <drtgfv$a8r$1@news.ox.ac.uk>, Ben Radford
<benjamin.radford@new.ox.ac.uk> writes[color=blue]
>Richard Herring wrote:
><snip>
>[color=green][color=darkred]
>>>> What is wrong with using the standard function std::atan2(), as
>>>>Jerry Coffin suggested? In a single function call it takes care of
>>>>all four quadrants with no sign checking. And it really will work on
>>>>all combinations of inputs except the indeterminate case (0, 0).
>>>
>>>
>>> The only case my function won't work on is when the square causes an
>>>overflow as you pointed out.[/color]
>> That's a whole set of cases more than is necessary.
>>[color=darkred]
>>> This is a limitation of the architecture that can be worked around
>>>not a problem with the algorithm.[/color]
>> If the algorithm design doesn't take proper account of domains and
>>bounds then it _is_ a problem with the algorithm.[/color]
>
>Indeed, but only when those bounds apply in the general case and are
>not simply a limitation of the architecture. The bounds vary depending
>on the platform and are therefore an implementation and not algorithmic
>detail.[/color]

There's always a limit to the magnitude of floating-point numbers. Its
value is an implementation-specific detail, which you can find out using
tools like numeric_limits. Its existence is certainly not.
[color=blue]
>By your argument an algorithm that is implemented differently on two
>different platforms would actually be a different algorithm in each
>case because it would have to take proper account of the domains and
>bounds which apply for that case.[/color]

No. If you design the algorithm with the existence of bounds in mind,
you can produce a solution that runs without change on all platforms,
for all input values up to the fundamental limits imposed by the
architecture, regardless of what the actual bounds are.[color=blue]
>[color=green][color=darkred]
>>> However there is nothing wrong with what Jerry suggested, I just
>>>didn't examine his post clearly enough to realise exactly what he was
>>>saying. The only difference is that atan2 gives an angle in the range
>>>-PI/2 < theta < PI/2 whereas I believe the dot product method always
>>>gives a positive angle.[/color]
>> I don't have to believe it, I know it. If a positive angle is what
>>you need, std::abs() is your friend.[/color]
>
>I don't see where you are going with this. I acknowledged that Jerry's
>method was a valid way of calculating the answer the OP wanted. The
>algorithm I gave was just an alternative; albeit with an overflow
>problem in the hastily written implementation code - which you pointed
>out.[/color]

I'm simply trying to make the point that there's no point in proposing
and elaborating a flawed solution when the standard library provides a
simpler, error-free way to achieve the same result.
[color=blue]
>[color=green][color=darkred]
>>> Got to run now or I'm going to be late for a lecture.
>>>[/color]
>> Numerical Analysis 101?[/color]
>
>Nice insult, but totally opportunistic and unnecessary in conveying
>your point.
>[/color]
--
Richard Herring

On Thu, 02 Feb 2006 13:42:36 +0000, Ben Radford
<benjamin.radford@new.ox.ac.uk> wrote:
[color=blue]
>cwsullivan@ucdavis.edu wrote:[color=green]
>> Hi gang,
>> I have a grid full of particles in my program, and I want to find an
>> angle between two particles.[/color]
>
>What you proceed to describe is not "finding the angle between two
>particles" but finding the angle between the vectors a and b:
>
>a = p1 - p2
>b = (1, 0)
>
>where p1 and p2 are the particles.
>
>You can calculate this angle simply using the dot product equation.
>
>Recall:
>
>a.b = |a||b|cos(theta)
>
>where theta is the angle between the two vectors a and b. Rearrange this
>equation to get:
>
>theta = arccos((a.b)/(|a||b|))
>
>As a c function:
>
>double calcAngle(double p1_x, double p1_y, double p2_x, double p2_y)
>{
> double a_x = p2_x - p1_x;
> double a_y = p2_y - p1_y;
>
> double b_x = 1.0;
> double b_y = 0.0;
>
> return acos((a_x*b_x+a_y*b_y)/sqrt(a_x*a_x+a_y*a_y));
>}
>
>The only special case you have to worry about is when p1 == p2. The
>angle is clearly undefined in this case and this will manifest itself in
>a division by zero error when you try to divide by the magnitude. Apart
>from that this function should work for all cases with no messy sign
>checking.[/color]

Ben -- I like your solution the best. I blew some dust off a math
reference and found similar solution in three dimensions. Its more
general in that you can always specify the z-component as zero and
collapse the problem back down to two dimensions.

double temp1 = p1_x * p2_x + p1_y * p2_y + p1_z * p2_z;
double temp2 = p1_x * p1_x + p1_y * p1_y + p1_z * p1_z;
double temp3 = p2_x * p2_x + p2_y * p2_y + p2_z * p2_z;
return acos( temp1 / sqrt(temp2 * temp3) );

-- Tom