I've come to this a bit late but, for what it's worth, here's my response to the question of the intersection of two co-planar orbits.
Now, I'm going to break up the answer into two parts. The first part (in this post) will just focus on the simple case where the argument of periapsis of both orbits is the same. In the second part (the next post), I'll address the slightly more complicated case where the argument of periapsis of the two orbits differ by an amount [MATH]\Delta\omega[/MATH].
Simple case - the same argument of periapsis
For convenience, let's focus on the case where both orbits are elliptical and ignore the possibility that one or both orbits may be parabolic or hyperbolic. In this case we can adopt the same 'perifocal' coordinate system for both ellipses and write the equation for the two ellipses as:
[math]\left(x+a_1\,e_1\right){}^2+\frac{y^2}{1-e_1^2} =a_1^2[/math]
and
[math]\left(x+a_2\,e_2\right){}^2+\frac{y^2}{1-e_2^2} =a_2^2[/math]
where [imath]a_1[/imath] and [imath]e_1[/imath] are the semi-major axis and orbital eccentricity of the first ellipse; [imath]a_2[/imath] and [imath]e_2[/imath] are the semi-major axis and orbital eccentricity of the first ellipse; and [imath]x[/imath] and [imath]y[/imath] are the [imath]x[/imath]-[imath]y[/imath] coordinates of the ellipse in our perifocal coordinate system.
To solve these equations, we eliminate [imath]y[/imath] from both equations. This gives us a quadratic in [imath]x[/imath] which we can 'easily' solve. After some considerable algebra, we find that the two solutions for [imath]x[/imath] are:
[math]x = \frac{a_1+a_2-a_1\, e_1^2-a_2 \, e_2^2}{e_1+e_2}[/math]
and
[math]x = \frac{a_1-a_2-a_1\, e_1^2+a_2\, e_2^2}{e_1-e_2}[/math]
But we know that we are working in the perifocal reference frame so that we must have
[math]- a_1\,(1+e_1) \le x \le a_1\,(1-e_1)[/math]
and
[math]- a_2\,(1+e_2) \le x \le a_2\,(1-e_2)[/math]
And after some more algebra, we can show that the first solution for [imath]x[/imath] given above never satisfies these conditions while the second solution for [imath]x[/imath] does:
[math]x = \frac{a_1-a_2-a_1\, e_1^2+a_2\, e_2^2}{e_1-e_2}[/math]
provided that if [imath]0 \le e_2 < e_2 < 1[/imath] then [imath]\frac{1-e_1}{1-e_2}\leq \frac{a_2}{a_1}\leq \frac{1+e_1}{1+e_2}[/imath]; or if [imath]0 \le e_1 < e_2 < 1[/imath] then [imath]\frac{1+e_1}{1+e_2}\leq \frac{a_2}{a_1}\leq \frac{1-e_1}{1-e_2}[/imath]. If these conditions are not satisfied, then the two ellipses do not intersect.
Having solved for [imath]x[/imath], it is straightforward to plug the solution value back into
[math]\left(x+a_1\,e_1\right){}^2+\frac{y^2}{1-e_1^2} =a_1^2[/math]
and solve for [imath]y[/imath].
In the next post, I'll consider the more general case where the argument of periapsis of the two ellipses are not the same.
---------- Post added at 03:29 PM ---------- Previous post was at 06:53 AM ----------
In this post, I'll continue with the analytical solution to the problem - but this time assuming that the argument of periapsis of the two elliptical orbits are not the same and differ by an angle [imath]\Delta\nu[/imath].
More complicated case - different arguments of periapsis
Let's suppose that we adopt the perifocal reference frame for the first ellipse. In this reference frame, then, the argument of periapsis of the second ellipse is rotated counter-clockwise from the x-axis by an angle [imath]\Delta\omega[/imath] - i.e., the difference in the argument of periapsis of the two co-planar orbits.
Given this, the radius of a body moving on the first elliptical trajectory is given by:
[math]r = \frac{a_1 \,\left(1-e_1^2\right)}{1+e_1\,\cos (\nu )}[/math]
and the radius of a body moving in the second elliptical trajectory is given by :
[math]r = \frac{a_2 \,\left(1-e_2^2\right)}{1+e_2\,\cos (\nu -\Delta\omega)}[/math]
For an intersection of the two ellipses to occur, we require that the these two values the same. This leads to a quadratic equation in the free variable [imath]r[/imath] by expanding the trig function on the left-hand side and eliminating [imath]\cos\nu[/imath].
To keep the algebra manageable, let's calculate a few intermediate values.
[math]\beta _1 = a_1 \, e_2 \, \left(1-e_1^2\right)[/math]
[math]\beta _2 = a_2 \, e_1 \, \left(1-e_2^2\right)[/math]
[math]A = \beta _1 \, e_2+\beta _2 \, e_1-\left(\beta _1 \, e_1+\beta _2 \, e_2\right) \cos (\Delta \omega )[/math]
[math]B = e_1^2+e_2^2-2\, e_1 \,e_2\, \cos (\Delta \omega )-e_1^2 \,e_2^2 \,\sin ^2(\Delta \omega )[/math]
[math]\Delta = \sin ^2(\Delta\omega ) \left(\left(\beta _1^2-2 \,\cos (\Delta\omega ) \,\beta _1 \,\beta _2+\beta _2^2\right)\, e_1^2 \,e_2^2-\left(\beta _1\, e_1-\beta _2 \,e_2\right){}^2\right)[/math]
Next, let's double check that [imath]\Delta\ge 0[/imath]. If it isn't, we can immediately say that the two ellipses do not intersect. If, on the other hand, it is then we calculate the two roots of the quadratic in [imath]r[/imath]:
[math]r_+ = \frac{A+\sqrt{\Delta}}{B}[/math]
and
[math]r_- = \frac{A-\sqrt{\Delta}}{B}[/math]
Now, before we go on, we need to double check that these are valid solutions such that:
[math]a_1\,(1 - e_1) \le r_+ \le a_1\,(1 + e_1)[/math]
[math]a_2\,(1 - e_2) \le r_+ \le a_2\,(1 + e_2)[/math]
and
[math]a_1\,(1 - e_1) \le r_- \le a_1\,(1 + e_1)[/math]
[math]a_2\,(1 - e_2) \le r_- \le a_2\,(1 + e_2)[/math]
If either [imath]r_+[/imath] or [imath]r_-[/imath] fail to satisfy these requirements (that the radius of the point of intersection lie between the periapsis radius and apoapsis radius of both orbits) then we throw that solution away. Assuming that both satisfy these requirements, then calculate:
[math]c_+ = \frac{a_1 \,\left(1-e_1^2\right)-r_+}{r_+ \,e_1}[/math]
[math]s_+ = \frac{\left(a_2 \,\left(1-e_2^2\right)-r_+\right) }{r_+\, e_2}\,\csc (\Delta\omega )-\frac{\left(a_1 \,\left(1-e_1^2\right)-r_+\right) }{r_+\, e_1}\,\cot (\Delta\omega )[/math]
and
[math]c_- = \frac{a_1 \,\left(1-e_1^2\right)-r_-}{r_- \,e_1}[/math]
[math]s_- = \frac{\left(a_2 \,\left(1-e_2^2\right)-r_-\right) }{r_-\, e_2}\,\csc (\Delta\omega )-\frac{\left(a_1 \,\left(1-e_1^2\right)-r_-\right) }{r_-\, e_1}\,\cot (\Delta\omega )[/math]
So, in cartesian coordinates, the two points of intersection of the two ellipses are given by:
[math](x_+, y_+) = r_+\,(c_+, s_+)[/math]
and
[math](x_-, y_-) = r_-\,(c_-, s_-)[/math]
In a third and final post, I'll just wrap this discussion by working through an example.
[math](x_+,\, y_+) = r_+\,(c_+, \,s_-)[/math]
---------- Post added at 04:02 PM ---------- Previous post was at 03:29 PM ----------
As an example of the above, let's assume that we have two ellipses - the first with semi-major axis and eccentricity:
[math]a_1 = 1.0[/math][math]e_1 = 0.9[/math]
and the second with semi-major axis and eccentricity:
[math]a_1 = 1.2[/math][math]e_1 = 0.1[/math]
And, let's also suppose that the argument of periapsis of the second ellipse is rotated by 45 degrees counter-clockwise from the first ellipse - i.e.,
[math]\Delta\omega = \pi/4[/math]
Let's now calculate the intermediate values:
[math]\beta _1 = a_1 \, e_2 \, \left(1-e_1^2\right) = 0.019[/math]
[math]\beta _2 = a_2 \, e_1 \, \left(1-e_2^2\right) = 1.0692[/math]
[math]A = \beta _1 \, e_2+\beta _2 \, e_1-\left(\beta _1 \, e_1+\beta _2 \, e_2\right) \cos (\Delta \omega ) = 0.876484616997[/math]
[math]B = e_1^2+e_2^2-2\, e_1 \,e_2\, \cos (\Delta \omega )-e_1^2 \,e_2^2 \,\sin ^2(\Delta \omega ) = 0.688670779386[/math]
[math]\Delta = \sin ^2(\Delta\omega ) \left(\left(\beta _1^2-2 \,\cos (\Delta\omega ) \,\beta _1 \,\beta _2+\beta _2^2\right)\, e_1^2 \,e_2^2-\left(\beta _1\, e_1-\beta _2 \,e_2\right){}^2\right) = 0.00048120550600[/math]
(We note that [imath]\Delta > 0[/imath], so we go on.)
Calculate the two solutions for [imath]r[/imath]:
[math]r_+ = \frac{A+\sqrt{\Delta}}{B}= 1.3045725776699[/math]
and
[math]r_- = \frac{A-\sqrt{\Delta}}{B}=1.240866094138270[/math]
We now double check that these two values lie between the periapsis and apoapsis of both orbits and confirm that indeed they do.
Now, we calculate the sines and cosines of the true anomaly of the point of intersection (in the perifocal reference frame of the first ellipse):
[math]c_+ = \frac{a_1 \,\left(1-e_1^2\right)-r_+}{r_+ \,e_1} = -0.9492871430738[/math]
[math]s_+ = \frac{\left(a_2 \,\left(1-e_2^2\right)-r_+\right) }{r_+\, e_2}\,\csc (\Delta\omega )-\frac{\left(a_1 \,\left(1-e_1^2\right)-r_+\right) }{r_+\, e_1}\,\cot (\Delta\omega ) = -0.314410432388360[/math]
and
[math]c_- = \frac{a_1 \,\left(1-e_1^2\right)-r_-}{r_- \,e_1} = -0.9409790460088[/math]
[math]s_- = \frac{\left(a_2 \,\left(1-e_2^2\right)-r_-\right) }{r_-\, e_2}\,\csc (\Delta\omega )-\frac{\left(a_1 \,\left(1-e_1^2\right)-r_-\right) }{r_-\, e_1}\,\cot (\Delta\omega )=0.338464820878[/math]
So that the two points of intersection are:
[math](x_+,\, y_+) = r_+\,(c_+,\, s_+) = (-1.238413975188,-0.41017122822)[/math]
and
[math](x_-, \,y_-) = r_-\,(c_-,\, s_-) = (-1.1676289934869,0.419989520286614)[/math]
---------- Post added 03-16-18 at 02:29 AM ---------- Previous post was 03-15-18 at 04:02 PM ----------
And here is a small Lua code snippet for calculating the points of intersection.
(Six inputs are required: the semi-major axis, the eccentricity and the argument of periapsis of the two coplanar elliptical orbits. The code snippet either informs that there are no intersections or returns the two points of intersection in the x-y coordinates of the orbital plane of the two orbits.)
Code:
-- define the inputs for the first ellipse
-- a1 - the semi-major axis
-- e1 - the orbital eccentricity
-- w1 - the argiment of periapsis
a1 = 1.0
e1 = 0.9
w1 = 0.0
-- define the inputs for the second ellipse
-- a2 - the semi-major axis
-- e3 - the orbital eccentricity
-- w2 - the argiment of periapsis
a2 = 1.2
e2 = 0.1
w2 = math.pi / 4.0
-- calculate the intermediary quantities
dw = w2 - w1
c1 = math.cos(dw)
s1 = math.sin(dw)
s2 = s1 * s1
b1 = a1 * e2 * (1 - e1 * e1 )
b2 = a2 * e1 * (1 - e2 * e2 )
k = b1 * e1 + b2 * e2
q = e1 * e1 * e2 * e2
l = b1 * e1 - b2 * e2
A = b1 * e2 + b2 * e1 - k * c1
B = e1 * e1 + e2 * e2 - 2 * e1 * e2 * c1 - q * s2
D = s2 * ( (b1 * b1 - 2 * c1 * b1 * b2 + b2 * b2) * q - l * l )
if D > 0 then
rp = (A + math.sqrt(D)) / B
rm = (A - math.sqrt(D)) / B
c1p = (a1 * (1 - e1 * e1) - rp) / rp / e1
c2p = (a2 * (1 - e2 * e2) - rp) / rp / e2
c1m = (a1 * (1 - e1 * e1) - rm) / rm / e1
c2m = (a2 * (1 - e2 * e2) - rm) / rm / e2
cp = c1p
sp = (c2p - c1p * c1 ) / s1
cm = c1m
sm = (c2m - c1m * c1 ) / s1
xp = rp * cp
yp = rp * sp
xm = rm * cm
ym = rm * sm
Xp = math.cos(w1) * xp - math.sin(w1) * yp
Yp = -math.sin(w1) * xp + math.cos(w1) * yp
Xm = math.cos(w1) * xm - math.sin(w1) * ym
Ym = -math.sin(w1) * xm + math.cos(w1) * ym
print()
print("The first point of intersection is: (", Xp, ", ", Yp, ")")
print("The second point of intersection is: (", Xm, ", ", Ym, ")")
print()
else
print("There are no points of intersection")
end
For the example given above, this code snippet returns:
Code:
The first point of intersection is: ( -1.2384139751888 , -0.4101712282272 )
The second point of intersection is: ( -1.167628993487 , 0.41998952028661 )
which are just the solutions given earlier.
