Orbiter-Forum  

Go Back   Orbiter-Forum > Far Side of the Moon > Math & Physics
Register Blogs Orbinauts List Social Groups FAQ Projects Mark Forums Read

Math & Physics Mathematical and physical problems of space flight and astronomy.

Reply
 
Thread Tools
Old 12-18-2017, 03:12 PM   #1
maki
Orbinaut
Default calculating intesection of two orbits

I would like to calculate intersection of two orbits in same plane. I had my last math class 15 years ago and i do not use it for living, so if i assumed something wrong and/made some silly mistake(s) please point it out but be gentle :-)
If i have two coplanar orbits with known semimajor axis eccentricity and argument of periapsis(i beleive it os emough to define a unique orbit in a plane) i want to enter these values in equation and get angle from one of periapsis to intersection point(s). I tried and got a quadratic equation with cos of desired angle as unknown but solution is painfuly wrong. That means i made some mistakes but i cant find where, so i need help please. I will post my work in image. Thank you!
Finaly I get this but when I solve it I get something silly .

Click image for larger version

Name:	SmartSelectImage_2017-12-18-22-14-10.jpg
Views:	29
Size:	53.5 KB
ID:	15539
Attached Thumbnails
Notes_171218_140007_8d0_1.jpg   Notes_171218_140008_2d8_2.jpg   Notes_171218_215556_cc0_1.jpg  

Last edited by maki; 12-18-2017 at 09:18 PM.
maki is offline   Reply With Quote
Old 12-19-2017, 02:16 PM   #2
Miner1
Orbinaut
Default

The problem of finding intersection points of elliptical orbits has been discussed in another thread on this forum - https://www.orbiter-forum.com/showthread.php?t=14828.
Miner1 is offline   Reply With Quote
Thanked by:
Old 12-23-2017, 11:59 AM   #3
maki
Orbinaut
Default

Thank you. I searched forum for orbit intersection and didnt notice topic you pointed out. I used half angle trig identity you suggested immediately, but it took me a week to detect missing brackets in excel in denominator of quadratic equation solution formula, /2*a instead /(2*a) :-). It works exellent now!
maki is offline   Reply With Quote
Old 03-16-2018, 02:29 AM   #4
MontBlanc2012
Orbinaut
Default

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 \Delta\omega.

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:

\left(x+a_1\,e_1\right){}^2+\frac{y^2}{1-e_1^2} =a_1^2

and

\left(x+a_2\,e_2\right){}^2+\frac{y^2}{1-e_2^2} =a_2^2

where a_1 and e_1 are the semi-major axis and orbital eccentricity of the first ellipse; a_2 and e_2 are the semi-major axis and orbital eccentricity of the first ellipse; and x and y are the x-y coordinates of the ellipse in our perifocal coordinate system.

To solve these equations, we eliminate y from both equations. This gives us a quadratic in x which we can 'easily' solve. After some considerable algebra, we find that the two solutions for x are:

x = \frac{a_1+a_2-a_1\, e_1^2-a_2 \, e_2^2}{e_1+e_2}

and

x = \frac{a_1-a_2-a_1\, e_1^2+a_2\, e_2^2}{e_1-e_2}

But we know that we are working in the perifocal reference frame so that we must have

- a_1\,(1+e_1) \le x \le a_1\,(1-e_1)

and

- a_2\,(1+e_2) \le x \le a_2\,(1-e_2)


And after some more algebra, we can show that the first solution for x given above never satisfies these conditions while the second solution for x does:

x = \frac{a_1-a_2-a_1\, e_1^2+a_2\, e_2^2}{e_1-e_2}

provided that if 0 \le e_2 <  e_2 < 1 then \frac{1-e_1}{1-e_2}\leq \frac{a_2}{a_1}\leq \frac{1+e_1}{1+e_2}; or if 0 \le e_1 <  e_2 < 1 then \frac{1+e_1}{1+e_2}\leq \frac{a_2}{a_1}\leq \frac{1-e_1}{1-e_2}. If these conditions are not satisfied, then the two ellipses do not intersect.

Having solved for x, it is straightforward to plug the solution value back into

\left(x+a_1\,e_1\right){}^2+\frac{y^2}{1-e_1^2} =a_1^2

and solve for y.

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 [MATH\\Delta\nu[/MATH].

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 \Delta\omega - 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:

r = \frac{a_1 \,\left(1-e_1^2\right)}{1+e_1\,\cos (\nu )}

and the radius of a body moving in the second elliptical trajectory is given by :

r = \frac{a_2 \,\left(1-e_2^2\right)}{1+e_2\,\cos (\nu -\Delta\omega)}

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 r by expanding the trig function on the left-hand side and eliminating \cos\nu.

To keep the algebra manageable, let's calculate a few intermediate values.

\beta _1 = a_1 \, e_2 \, \left(1-e_1^2\right)

\beta _2 = a_2 \, e_1 \, \left(1-e_2^2\right)

A =  \beta _1 \, e_2+\beta _2 \, e_1-\left(\beta _1 \, e_1+\beta _2 \, e_2\right) \cos (\Delta \omega )

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 )

\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)

Next, let's double check that \Delta\ge 0. 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 r:

r_+ = \frac{A+\sqrt{\Delta}}{B}

and

r_- = \frac{A-\sqrt{\Delta}}{B}

Now, before we go on, we need to double check that these are valid solutions such that:

a_1\,(1 - e_1) \le r_+ \le a_1\,(1 + e_1)

a_2\,(1 - e_2) \le r_+ \le a_2\,(1 + e_2)

and

a_1\,(1 - e_1) \le r_- \le a_1\,(1 + e_1)

a_2\,(1 - e_2) \le r_- \le a_2\,(1 + e_2)


If either r_+ or r_- 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:

c_+ = \frac{a_1 \,\left(1-e_1^2\right)-r_+}{r_+ \,e_1}

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 )

and

c_- = \frac{a_1 \,\left(1-e_1^2\right)-r_-}{r_- \,e_1}

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 )

So, in cartesian coordinates, the two points of intersection of the two ellipses are given by:

(x_+, y_+) = r_+\,(c_+, s_+)

and

(x_-, y_-) = r_-\,(c_-, s_-)

In a third and final post, I'll just wrap this discussion by working through an example.

(x_+, y_+) = r_+\,(c_+, s_-)

---------- 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:

a_1 = 1.0
e_1 = 0.9

and the second with semi-major axis and eccentricity:

a_1 = 1.2
e_1 = 0.1

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.,

\Delta\omega = \pi/4


Let's now calculate the intermediate values:

\beta _1 = a_1 \, e_2 \, \left(1-e_1^2\right) = 0.019

\beta _2 = a_2 \, e_1 \, \left(1-e_2^2\right) = 1.0692

A =  \beta _1 \, e_2+\beta _2 \, e_1-\left(\beta _1 \, e_1+\beta _2 \, e_2\right) \cos (\Delta \omega ) = 0.876484616997

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

\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

(We note that \Delta > 0, so we go on.)

Calculate the two solutions for r:


r_+ = \frac{A+\sqrt{\Delta}}{B}= 1.3045725776699

and

r_- = \frac{A-\sqrt{\Delta}}{B}=1.240866094138270

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):


c_+ = \frac{a_1 \,\left(1-e_1^2\right)-r_+}{r_+ \,e_1} = -0.9492871430738

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

and

c_- = \frac{a_1 \,\left(1-e_1^2\right)-r_-}{r_- \,e_1} = -0.9409790460088

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

So that the two points of intersection are:


(x_+, y_+) = r_+\,(c_+, s_+) = (-1.238413975188,-0.41017122822)

and

(x_-, y_-) = r_-\,(c_-, s_-) = (-1.1676289934869,0.419989520286614)

---------- 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.

Last edited by MontBlanc2012; 03-16-2018 at 02:32 AM.
MontBlanc2012 is offline   Reply With Quote
Old 09-28-2018, 03:15 PM   #5
malakai
Orbinaut
Default

Thanks for posting this excellent response. I really appreciate seeing the solution written in "code" as well as in "math". I find that most of my problems with implementing the ideas I find documented in textbooks come from mis-reading the "math", usually because there is some additional context to the symbols being used that I've missed. I find "code" a lot easier to read, and seeing a side-by-side translation really helps a ton!

Anyway, using your description and the lua code as a reference, I was able to implement the solution in javascript for a project that I'm doing, and it works flawlessly!

Two things I'd like to point out-

For the "different arguments of periapsis" case, I made a mistake in my implementation. I had transformed the first orbit into perifocal reference, but I forgot to rotate the second orbit by the same amount to match. Obvious in hindsight, but it threw me for a loop trying to debug what I'd done wrong.

Second, for anyone trying to find a similar solution for the more complicated case of mixed hyperbolic and elliptical orbits.... its not more complicated. The math is the same, and the input orbits can be elliptical or hyperbolic and the solution still works.

Thanks again, MontBlanc2012!
malakai is offline   Reply With Quote
Thanked by:
Old 12-01-2018, 05:32 AM   #6
MontBlanc2012
Orbinaut
Default

No problem, malakai.

Re extension to mixing elliptical and hyperbolic orbit: as I worked through the maths, I though that this extension was likely but I didn't have time to confirm it. So, thanks for considering these cases and confirming that it also works in the hyperbolic regime and (by extension) the mixed regime as well.

[Sorry for the belated response, but sadly I've been busy with other things for the last month or so and I haven't had time to read much that has been posted in this forum.]
MontBlanc2012 is offline   Reply With Quote
Reply

  Orbiter-Forum > Far Side of the Moon > Math & Physics


Thread Tools

Posting Rules
BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts
Forum Jump


All times are GMT. The time now is 01:00 PM.

Quick Links Need Help?


About Us | Rules & Guidelines | TOS Policy | Privacy Policy

Orbiter-Forum is hosted at Orbithangar.com
Powered by vBulletin® Version 3.8.6
Copyright ©2000 - 2018, Jelsoft Enterprises Ltd.
Copyright 2007 - 2017, Orbiter-Forum.com. All rights reserved.