Application Center - Maplesoft
  Submit your work
Application Center Applications The Intersection of a Line and Cone Preview

App Preview:

The Intersection of a Line and Cone

You can switch back to the summary page by clicking here.

Learn about Maple
Download Application

 

Ray and Object Intersections: Truncated Cone


by Otto Wilke
otto_wilke@hotmail.com

I am vaguely aware that graphics is normally done with vector operations, generic

solids positioned at the origin, and transformation matrices to move rays to and fro.

I thought it would be interesting to use rectangular coordinates and objects located

anywhere in space and oriented in any direction.

This is one of four files covering the plane, the sphere, the cylinder, and the cone.

>

Let C be a truncated right circular cone.  Let the larger circular end be centered at

P1(x1,y1,z1) and have radius R.  Let the smaller circular end be centered at

P2(x2,y2,z2) and have radius r.  Let P(x,y,z) be some point on the surface of the cone.

Let P3(x3,y3,z3) be a point on line segment P1 P2 such that line segment P P3 is

perpendicular to P1 P2.  Let d1 be distance P1 P2.  Let d2 be distance P P1.

Let d3 be distance P P2.  Let d4 be distance P P3.  Let d5 be distance P1 P3.

Let d6 be distance P2 P3.  By the distance formula.

> restart;with(plots):
draw1:=plot([[0,0],[1,4],[3,4],[4,0],[0,0],[2,0],[2,4],[3.5,2],
[2,0],[2,2],[3.5,2]],axes=none):
text1:=textplot([[2,-.3,"P1"],[2,4.3,"P2"],[3.6,2,"P"],
[2.6,3.8,"r"],[3,.2,"R"],[1.8,2,"P3"]]):
display({draw1,text1});
d1=sqrt((x2-x1)^2+(y2-y1)^2+(z2-z1)^2);
d2=sqrt((x-x1)^2+(y-y1)^2+(z-z1)^2);
d3=sqrt((x-x2)^2+(y-y2)^2+(z-z2)^2);

Warning, the name changecoords has been redefined

[Plot]

d1 = (x2^2-2*x2*x1+x1^2+y2^2-2*y2*y1+y1^2+z2^2-2*z2*z1+z1^2)^(1/2)

d2 = (x^2-2*x*x1+x1^2+y^2-2*y*y1+y1^2+z^2-2*z*z1+z1^2)^(1/2)

d3 = (x^2-2*x*x2+x2^2+y^2-2*y*y2+y2^2+z^2-2*z*z2+z2^2)^(1/2)

I think I remember that

> d4=r+(d6/d1)*(R-r);solve(d4=r+(d6/d1)*(R-r),d6);%^2;

d4 = r+d6*(R-r)/d1

-(-d4+r)*d1/(R-r)

(-d4+r)^2*d1^2/(R-r)^2

By the Pythagorean Theorem

> d4=sqrt(d3^2-d6^2);d4=sqrt(d3^2-((-d4+r)^2*d1^2/(-R+r)^2));

d4 = (d3^2-d6^2)^(1/2)

d4 = (d3^2-(-d4+r)^2*d1^2/(-R+r)^2)^(1/2)

Solving for d4

> {solve(d4 = (d3^2-(-d4+r)^2*d1^2/(-R+r)^2)^(1/2),d4)};

{(d1^2*r+(d1^2*d3^2*R^2-2*d1^2*d3^2*R*r+d1^2*d3^2*r^2+R^4*d3^2-4*R^3*d3^2*r+6*R^2*d3^2*r^2-R^2*d1^2*r^2-4*R*r^3*d3^2+2*R*r^3*d1^2+r^4*d3^2-r^4*d1^2)^(1/2))/(d1^2+R^2-2*R*r+r^2), (d1^2*r-(d1^2*d3^2*R^2...{(d1^2*r+(d1^2*d3^2*R^2-2*d1^2*d3^2*R*r+d1^2*d3^2*r^2+R^4*d3^2-4*R^3*d3^2*r+6*R^2*d3^2*r^2-R^2*d1^2*r^2-4*R*r^3*d3^2+2*R*r^3*d1^2+r^4*d3^2-r^4*d1^2)^(1/2))/(d1^2+R^2-2*R*r+r^2), (d1^2*r-(d1^2*d3^2*R^2...

Knowing d4 must be positive, d4 is

> op(1,%);

(d1^2*r+(d1^2*d3^2*R^2-2*d1^2*d3^2*R*r+d1^2*d3^2*r^2+R^4*d3^2-4*R^3*d3^2*r+6*R^2*d3^2*r^2-R^2*d1^2*r^2-4*R*r^3*d3^2+2*R*r^3*d1^2+r^4*d3^2-r^4*d1^2)^(1/2))/(d1^2+R^2-2*R*r+r^2)

The area of triangle P P1 P2 is half of the product of base and height

> aArea=1/2*d1*d4;areaSquared=(rhs(%))^2;

aArea = 1/2*d1*d4

areaSquared = 1/4*d1^2*d4^2

By Heron's formula, the area squared is also

> areaSquared=(d1+d2+d3)/2*((d1+d2+d3)/2-d1)*((d1+d2+d3)/2-d2)
*((d1+d2+d3)/2-d3);

areaSquared = 1/2*(d1+d2+d3)*(-1/2*d1+1/2*d2+1/2*d3)*(1/2*d1-1/2*d2+1/2*d3)*(1/2*d1+1/2*d2-1/2*d3)

Setting the squares of the area equal to each other

> (1/2*d1*d4)^2=(d1+d2+d3)/2*((d1+d2+d3)/2-d1)*((d1+d2+d3)/2-d2)
*((d1+d2+d3)/2-d3);

1/4*d1^2*d4^2 = 1/2*(d1+d2+d3)*(-1/2*d1+1/2*d2+1/2*d3)*(1/2*d1-1/2*d2+1/2*d3)*(1/2*d1+1/2*d2-1/2*d3)

Inserting the desired value of d4

> (1/2*d1*(1/2*(2*r*d1^2-2*(d1^2*d3^2*r^2+d1^2*d3^2*R^2-2*d1^2*
d3^2*R*r+6*R^2*d3^2*r^2-R^2*d1^2*r^2+R^4*d3^2-4*R^3*d3^2*r
-4*R*r^3*d3^2+2*R*r^3*d1^2+r^4*d3^2-r^4*d1^2)^(1/2))/
(d1^2+R^2-2*R*r+r^2)))^2=(d1+d2+d3)/2*((d1+d2+d3)/2-d1)*
((d1+d2+d3)/2-d2)*((d1+d2+d3)/2-d3);

1/16*d1^2*(2*d1^2*r-2*(d1^2*d3^2*R^2-2*d1^2*d3^2*R*r+d1^2*d3^2*r^2+R^4*d3^2-4*R^3*d3^2*r+6*R^2*d3^2*r^2-R^2*d1^2*r^2-4*R*r^3*d3^2+2*R*r^3*d1^2+r^4*d3^2-r^4*d1^2)^(1/2))^2/(d1^2+R^2-2*R*r+r^2)^2 = 1/2*...1/16*d1^2*(2*d1^2*r-2*(d1^2*d3^2*R^2-2*d1^2*d3^2*R*r+d1^2*d3^2*r^2+R^4*d3^2-4*R^3*d3^2*r+6*R^2*d3^2*r^2-R^2*d1^2*r^2-4*R*r^3*d3^2+2*R*r^3*d1^2+r^4*d3^2-r^4*d1^2)^(1/2))^2/(d1^2+R^2-2*R*r+r^2)^2 = 1/2*...

> 4*%;

1/4*d1^2*(2*d1^2*r-2*(d1^2*d3^2*R^2-2*d1^2*d3^2*R*r+d1^2*d3^2*r^2+R^4*d3^2-4*R^3*d3^2*r+6*R^2*d3^2*r^2-R^2*d1^2*r^2-4*R*r^3*d3^2+2*R*r^3*d1^2+r^4*d3^2-r^4*d1^2)^(1/2))^2/(d1^2+R^2-2*R*r+r^2)^2 = 2*(d1...1/4*d1^2*(2*d1^2*r-2*(d1^2*d3^2*R^2-2*d1^2*d3^2*R*r+d1^2*d3^2*r^2+R^4*d3^2-4*R^3*d3^2*r+6*R^2*d3^2*r^2-R^2*d1^2*r^2-4*R*r^3*d3^2+2*R*r^3*d1^2+r^4*d3^2-r^4*d1^2)^(1/2))^2/(d1^2+R^2-2*R*r+r^2)^2 = 2*(d1...

Substituting for d1, d2, and d3

> subs(d1=sqrt((x2-x1)^2+(y2-y1)^2+(z2-z1)^2),
d2=sqrt((x-x1)^2+(y-y1)^2+(z-z1)^2),
d3=sqrt((x-x2)^2+(y-y2)^2+(z-z2)^2),%);

1/4*(x2^2-2*x2*x1+x1^2+y2^2-2*y2*y1+y1^2+z2^2-2*z2*z1+z1^2)*(2*(x2^2-2*x2*x1+x1^2+y2^2-2*y2*y1+y1^2+z2^2-2*z2*z1+z1^2)*r-2*((x2^2-2*x2*x1+x1^2+y2^2-2*y2*y1+y1^2+z2^2-2*z2*z1+z1^2)*(x^2-2*x*x2+x2^2+y^2...1/4*(x2^2-2*x2*x1+x1^2+y2^2-2*y2*y1+y1^2+z2^2-2*z2*z1+z1^2)*(2*(x2^2-2*x2*x1+x1^2+y2^2-2*y2*y1+y1^2+z2^2-2*z2*z1+z1^2)*r-2*((x2^2-2*x2*x1+x1^2+y2^2-2*y2*y1+y1^2+z2^2-2*z2*z1+z1^2)*(x^2-2*x*x2+x2^2+y^2...1/4*(x2^2-2*x2*x1+x1^2+y2^2-2*y2*y1+y1^2+z2^2-2*z2*z1+z1^2)*(2*(x2^2-2*x2*x1+x1^2+y2^2-2*y2*y1+y1^2+z2^2-2*z2*z1+z1^2)*r-2*((x2^2-2*x2*x1+x1^2+y2^2-2*y2*y1+y1^2+z2^2-2*z2*z1+z1^2)*(x^2-2*x*x2+x2^2+y^2...1/4*(x2^2-2*x2*x1+x1^2+y2^2-2*y2*y1+y1^2+z2^2-2*z2*z1+z1^2)*(2*(x2^2-2*x2*x1+x1^2+y2^2-2*y2*y1+y1^2+z2^2-2*z2*z1+z1^2)*r-2*((x2^2-2*x2*x1+x1^2+y2^2-2*y2*y1+y1^2+z2^2-2*z2*z1+z1^2)*(x^2-2*x*x2+x2^2+y^2...1/4*(x2^2-2*x2*x1+x1^2+y2^2-2*y2*y1+y1^2+z2^2-2*z2*z1+z1^2)*(2*(x2^2-2*x2*x1+x1^2+y2^2-2*y2*y1+y1^2+z2^2-2*z2*z1+z1^2)*r-2*((x2^2-2*x2*x1+x1^2+y2^2-2*y2*y1+y1^2+z2^2-2*z2*z1+z1^2)*(x^2-2*x*x2+x2^2+y^2...1/4*(x2^2-2*x2*x1+x1^2+y2^2-2*y2*y1+y1^2+z2^2-2*z2*z1+z1^2)*(2*(x2^2-2*x2*x1+x1^2+y2^2-2*y2*y1+y1^2+z2^2-2*z2*z1+z1^2)*r-2*((x2^2-2*x2*x1+x1^2+y2^2-2*y2*y1+y1^2+z2^2-2*z2*z1+z1^2)*(x^2-2*x*x2+x2^2+y^2...1/4*(x2^2-2*x2*x1+x1^2+y2^2-2*y2*y1+y1^2+z2^2-2*z2*z1+z1^2)*(2*(x2^2-2*x2*x1+x1^2+y2^2-2*y2*y1+y1^2+z2^2-2*z2*z1+z1^2)*r-2*((x2^2-2*x2*x1+x1^2+y2^2-2*y2*y1+y1^2+z2^2-2*z2*z1+z1^2)*(x^2-2*x*x2+x2^2+y^2...1/4*(x2^2-2*x2*x1+x1^2+y2^2-2*y2*y1+y1^2+z2^2-2*z2*z1+z1^2)*(2*(x2^2-2*x2*x1+x1^2+y2^2-2*y2*y1+y1^2+z2^2-2*z2*z1+z1^2)*r-2*((x2^2-2*x2*x1+x1^2+y2^2-2*y2*y1+y1^2+z2^2-2*z2*z1+z1^2)*(x^2-2*x*x2+x2^2+y^2...1/4*(x2^2-2*x2*x1+x1^2+y2^2-2*y2*y1+y1^2+z2^2-2*z2*z1+z1^2)*(2*(x2^2-2*x2*x1+x1^2+y2^2-2*y2*y1+y1^2+z2^2-2*z2*z1+z1^2)*r-2*((x2^2-2*x2*x1+x1^2+y2^2-2*y2*y1+y1^2+z2^2-2*z2*z1+z1^2)*(x^2-2*x*x2+x2^2+y^2...1/4*(x2^2-2*x2*x1+x1^2+y2^2-2*y2*y1+y1^2+z2^2-2*z2*z1+z1^2)*(2*(x2^2-2*x2*x1+x1^2+y2^2-2*y2*y1+y1^2+z2^2-2*z2*z1+z1^2)*r-2*((x2^2-2*x2*x1+x1^2+y2^2-2*y2*y1+y1^2+z2^2-2*z2*z1+z1^2)*(x^2-2*x*x2+x2^2+y^2...1/4*(x2^2-2*x2*x1+x1^2+y2^2-2*y2*y1+y1^2+z2^2-2*z2*z1+z1^2)*(2*(x2^2-2*x2*x1+x1^2+y2^2-2*y2*y1+y1^2+z2^2-2*z2*z1+z1^2)*r-2*((x2^2-2*x2*x1+x1^2+y2^2-2*y2*y1+y1^2+z2^2-2*z2*z1+z1^2)*(x^2-2*x*x2+x2^2+y^2...1/4*(x2^2-2*x2*x1+x1^2+y2^2-2*y2*y1+y1^2+z2^2-2*z2*z1+z1^2)*(2*(x2^2-2*x2*x1+x1^2+y2^2-2*y2*y1+y1^2+z2^2-2*z2*z1+z1^2)*r-2*((x2^2-2*x2*x1+x1^2+y2^2-2*y2*y1+y1^2+z2^2-2*z2*z1+z1^2)*(x^2-2*x*x2+x2^2+y^2...

The equation immediately above give the coordinates of points on the truncated cone.

Let's graph a specific cone.

> x1:=4;y1:=2;z1:=-2;R:=2;x2:=0;y2:=6;z2:=2;r:=1;

x1 := 4

y1 := 2

z1 := -2

R := 2

x2 := 0

y2 := 6

z2 := 2

r := 1

> d1:=sqrt((x2-x1)^2+(y2-y1)^2+(z2-z1)^2);
d2:=sqrt((x-x1)^2+(y-y1)^2+(z-z1)^2);
d3:=sqrt((x-x2)^2+(y-y2)^2+(z-z2)^2);

d1 := 4*3^(1/2)

d2 := (x^2-8*x+24+y^2-4*y+z^2+4*z)^(1/2)

d3 := (x^2+y^2-12*y+40+z^2-4*z)^(1/2)

> with(plots):

> plot1:=implicitplot3d(1/4*d1^2*(2*r*d1^2-2*(d1^2*d3^2*r^2+d1^2
*d3^2*R^2-2*d1^2*d3^2*R*r+6*R^2*d3^2*r^2-R^2*d1^2*r^2+R^4*
d3^2-4*R^3*d3^2*r-4*R*r^3*d3^2+2*R*r^3*d1^2+r^4*d3^2-r^4*d1^2)
^(1/2))^2/(d1^2+R^2-2*R*r+r^2)^2 = 2*(d1+d2+d3)*(-1/2*d1+1/2*d2+1/2*d3)*(1/2*d1-1/2*d2+1/2*d3)
*(1/2*d1+1/2*d2-1/2*d3), x=-5..10, y=-5..10, z=-5..10, axes=normal,grid=[40,40,40]):
plot2:=pointplot3d({[x1,y1,z1],[x2,y2,z2] },color=red,
symbol=box,symbolsize=20):
display({plot1,plot2});

[Plot]

What we see is a series of cones on the same axis, with some plotted points beyond

the ends of the original truncated cone.  To find the intersection of a ray and the truncated

cone, we need to find only points on the truncated cone.  This is ensured by the following

inequalities.

> d2<=sqrt(r^2+d1^2);d3<=sqrt(R^2+d1^2);

(x^2-8*x+24+y^2-4*y+z^2+4*z)^(1/2) <= 7

(x^2+y^2-12*y+40+z^2-4*z)^(1/2) <= 2*13^(1/2)

Let L be a line on point P4(x4,y4,z4) and let L have direction numbers a, b, and c.  The parametric equations for L are

> x=x4+a*t;y=y4+b*t;z=z4+c*t;

x = x4+a*t

y = y4+b*t

z = z4+c*t

Let's pick P4 as (0,0,0) and a=.5, b=2, c=.5.

> x4:=0;y4:=0;z4:=0;a:=.5;b:=2;c:=.5;

x4 := 0

y4 := 0

z4 := 0

a := .5

b := 2

c := .5

Now that we have equations for points on a line and a cone and two inequalities, we'll ask Maple to solve the simlutaneous equations and inequalities.  

> solve({x=x4+a*t,y=y4+b*t,z=z4+c*t,

> d2<=sqrt(r^2+d1^2),d3<=sqrt(R^2+d1^2),
1/4*d1^2*(2*r*d1^2-2*(d1^2*d3^2*r^2+d1^2*d3^2*R^2-2*d1^2*
d3^2*R*r+6*R^2*d3^2*r^2-R^2*d1^2*r^2+R^4*d3^2-4*R^3*d3^2
*r-4*R*r^3*d3^2+2*R*r^3*d1^2+r^4*d3^2-r^4*d1^2)^(1/2))^2/
(d1^2+R^2-2*R*r+r^2)^2 = 2*(d1+d2+d3)*(-1/2*d1+1/2*d2+1/2*d3)*(1/2*d1-1/2*d2+1/2*d3)
*(1/2*d1+1/2*d2-1/2*d3)},{x,y,z,t});

{y = 5.671263838, x = 1.417815960, t = 2.835631919, z = 1.417815960}, {t = 2.084722063, y = 4.169444126, x = 1.042361032, z = 1.042361032}

> plot3:=arrow(<0,0,0>, <2.5,10,2.5>,
width=[0.02, relative], head_length=[0.2, relative], color=red):
plot4:=pointplot3d({[1.417815960,5.671263838,1.417815960],
[1.042361032,4.169444126,1.042361032] },
color=blue,symbol=box,symbolsize=20):
display({plot1,plot2,plot3,plot4});

[Plot]

The ray pierces the cone, entering and exiting.  To determine when a ray hits one or both

circular ends, use the same procedure as was used for the cylinder.

>