Intersecting cone surfaces

Question

Given two cone surfaces with common vertex O in 3D:

Given are: Coordinates of point O Vectors OP, OL (OL is height for the green cone), OR (height of the blue cone), OC

I would like to know the angles P1RC and CRP2

My current approach is the following: Normalize vector OL (length = 1) Build a plane through L perpendicular to OL. That plane intersects the blue cone which results in creating an ellipse with center R and intersects the green cone in a circle with radius LP and center L. Now I have to find the intersecting points P1 and P2 in the built plane and knowing them I can find angles P1RC and CRP2. How to find the intersecting points of the circle and the ellipse P1 and P2 using Mathematica's functionalities?

N.B: The two cone planes will always be intersecting this way.

Here is Mathematica code that does the set-up in the picture above:

\[CapitalOmega] = {0, 0, 2};
rC1 = 1; rC2 = 2;
C1 = Cone[{{0, 0, -1}, \[CapitalOmega]}, rC1];
C2 = Cone[{{1, 1, -1}, \[CapitalOmega]}, rC2];
S = InfinitePlane[{0, 0, 0}, {{1, 0, 0}, {0, 1, 0}}];
L = RegionIntersection[InfiniteLine[{\[CapitalOmega], {0, 0, -1}}], S];
R = RegionIntersection[InfiniteLine[{\[CapitalOmega], {1, 1, -1}}], S];
iCP = RegionIntersection[C1, C2, S];
iCPD = iCP // DiscretizeRegion;
P1 = {-(Sqrt[2]/3), Sqrt[2]/3, 0}; P2 = {Sqrt[2]/3, -(Sqrt[2]/3), 0};

Graphics3D[{
  {Opacity[0.3], Darker@Green, C1},
  {Opacity[0.3], Darker@Blue, C2},
  {Opacity[0.5], LightRed, S},
  PointSize[0.01],
  Blue, L, Black, Text["L", L[[1]], {1, 1}],
  Blue, R, Black, Text["R", R[[1]], {1, 1}],
  Blue, Point[\[CapitalOmega]], Black, 
  Text["O", \[CapitalOmega], -{1, 1}],
  Red, Point[P1], Black, Text["P1", P1, {1, 1}],
  Red, Point[P2], Black, Text["P2", P2, {1, 1}],
  White, Line[{\[CapitalOmega], R[[1]]}], 
  Line[{\[CapitalOmega], L[[1]]}], Line[{\[CapitalOmega], P1}], 
  Line[{\[CapitalOmega], P2}], Line[{R[[1]], P1}], Line[{R[[1]], P2}]
  }, Lighting -> "Neutral", Boxed -> False]

(Using [CaptialOmega] because O is protected and easily confused with 0.)

Answer 1

In order to find the points P1 and P2 we can do the following steps.

1. Find the points L and R from the equations of cones C1 and C1 and plane S.

2. Find the intersections between the cones C1 and C2 and the plane S (that is perpendicular to OL).

   2.1. Denote those intersection is iC1S and iC2S respectively. Denote the boundary of iCxS with b(iCxS).

3. Both P1 and P2 belong to b(iC1S) and b(iC2S).

   I.e. both P1 or P2 belong to the boundary of the intesection of the cone C1 with S and to the boundary of the intersection of C2 with S. As explained in the question, one of these boundaries is a circle and the other an ellipse.

4. Hence we can find P1 and P2 as an intersection of these boundaries.

Let us do in Mathematica all of the steps listed above.

1. First the set-up

   Using "\[CapitalOmega]" instead of "O".

   [CapitalOmega] = {0, 0, 2}; rC1 = 1; rC2 = 2; C1 = Cone[{{0, 0, -1}, [CapitalOmega]}, rC1]; C2 = Cone[{{1, 1, -1}, [CapitalOmega]}, rC2]; S = InfinitePlane[{0, 0, 0}, {{1, 0, 0}, {0, 1, 0}}];

2. Find the interesections iCxS.

   iC1S = RegionIntersection[C1, S]; iC2S = RegionIntersection[C2, S];

3. Find the points L and R2.

   3.1. Since we cannot get the points by just using

   RegionIntersection[RegionBoundary[iC1S], RegionBoundary[iC2S]]

   we are going to use Maximize which is based on the observation that P1 and P2 are the points furthest away from the line passing through L and R and belonging to b(iC1S) and b(iC2S).

   Maximize would return one point Px which is one of P1 or P2. We will find the other point by using subtracting a small ball around Px from the constraints.

   3.2. This finds one of the points:

   sol = Maximize[{RegionDistance[InfiniteLine[{L[1], R[1]}], x]}, {Element[x, RegionBoundary[iC1S]], Element[x, RegionBoundary[iC2S]]}]

   (* Out[77]= {2/3, {x -> {-(Sqrt[2]/3), Sqrt[2]/3, 0}}} *)

   P1 = x /. sol[[2]]

   (* Out[78]= {-(Sqrt[2]/3), Sqrt[2]/3, 0} *)

   3.3. This finds the other:

   smallBall = ImplicitRegion[Total[({x, y, z} - P1)^2] < 1/5, {x, y, z}]; sol = Maximize[{RegionDistance[InfiniteLine[{L[1], R[1]}], x]}, Element[x, RegionDifference[ RegionIntersection[RegionBoundary[iC1S], RegionBoundary[iC2S]], smallBall]]]

   (* Out[84]= {2/3, {x -> {Sqrt[2]/3, -(Sqrt[2]/3), 0}}} *)

   P2 = x /. sol[[2]]

   (* Out[85]= {Sqrt[2]/3, -(Sqrt[2]/3), 0} *)

At this point we can plot the results:

g = Graphics3D[{
   {Opacity[0.3], Darker@Green, C1},
   {Opacity[0.3], Darker@Blue, C2},
   {Opacity[0.5], LightRed, S},
   PointSize[0.01],
   Blue, L, Black, Text["L", L[[1]], {1, 1}],
   Blue, R, Black, Text["R", R[[1]], {1, 1}],
   Blue, Point[\[CapitalOmega]], Black, 
   Text["O", \[CapitalOmega], -{1, 1}],
   Red, Point[P1], Black, Text["P1", P1, {1, 1}],
   Red, Point[P2], Black, Text["P2", P2, {1, 1}]}, Boxed -> False, 
  Lighting -> "Neutral"]