Conic intersection issue's, coming out a little inaccurate...

[Pagnatious]

With respect to finding the intersection point of two ellipses I am using the parameterized general ellipse equation to plot a point on one ellipse, then testing that points distances to the focal points of the second ellipse to see if it is a valid point on the second ellipse. As described here

It then suggests using the Newton-Raphson method to find the angles for which the point is on both ellipses. For that I need the derivative of the function. I haven't had to get that in a long time but here is my attempt and was wondering if anyone was willing to check my work.



[math]f(t) = \sqrt{u} + \sqrt{v} - 2a[/math]

[math]u = (h)^2 + (j)^2[/math]
[math](u)' = 2h*(h)' + 2j*(j)'[/math]

[math]v = (k)^2 + (l)^2[/math]
[math](v)' = 2k(k)' + 2l(l)'[/math]

[math]h = (X +AcosO*cos(t) - BsinO*sin(t)) - c[/math]
[math](h)' = 0 - AcosO* sin(t) - BsinO*cos(t) - 0[/math]

[math]j = (Y + AsinO*cos(t) + BcosO*sin(t)) - d[/math]
[math](j)' = 0 - AsinO*sin(t) + BcosO*cos(t) - 0[/math]

[math]k = (X +AcosO*cos(t) - BsinO*sin(t)) - c_2[/math]
[math](k)' = 0 - AcosO*sin(t) - BsinO*cos(t) - 0[/math]

[math]l = (Y + AsinO*cos(t) + BcosO*sin(t)) - d_2[/math]
[math](l)' = 0 - AsinO*sin(t) + BcosO*cos(t) - 0[/math]



[math]f(t) = \sqrt{u} + \sqrt{v} - 2a[/math]

[math]f'(t) = \frac{(u)'}{2\sqrt{u}} + \frac{(v)'}{2\sqrt{v}}[/math]

[math]f'(t) = \frac{2(h)(h)' + 2(j)(j)'}{2\sqrt{(h)^2 + (j)^2}} + \frac{2(k)(k)' + 2(l)(l)'}{2\sqrt{(k)^2 + (l)^2}}[/math]

[math]f'(t) = \frac{(h)(h)' + (j)(j)'}{\sqrt{(h)^2 + (j)^2}} + \frac{(k)(k)' + (l)(l)'}{\sqrt{(k)^2 + (l)^2}}[/math]

You can probably ignore the last line below, the monsterous one, as its just the above line with the h,j,k and l subbed in, but here it is anyway. It's too large to do in LaTeX even in divided in two parts.

f'(t) = \frac{2((X +AcosO*cos(t) - BsinO*sin(t)) - c)*(0 - AcosO* sin(t) - BsinO*cos(t)) + 2((Y + AsinO*cos(t) + BcosO*sin(t)) - d)*(0 - AsinO*sin(t) + BcosO*cos(t))}{2\sqrt{((X +AcosO*cos(t) - BsinO*sin(t)) - c)^2 + ((Y + AsinO*cos(t) + BcosO*sin(t)) - d)^2}} + \frac{2((X +AcosO*cos(t) - BsinO*sin(t)) - c_2)*(0 - AcosO* sin(t) - BsinO*cos(t)) + 2((Y + AsinO*cos(t) + BcosO*sin(t)) - d_2)*(0 - AsinO*sin(t) + BcosO*cos(t))}{2\sqrt{((X +AcosO*cos(t) - BsinO*sin(t)) - c_2)^2 + ((Y + AsinO*cos(t) + BcosO*sin(t)) - d_2)^2}}

Also I'd welcome any suggestions as to a good test to see first if they intersect so I can only apply the above when there will be an angle for which they intersect. They are orbits around the same body so they will share a focus and should at most intersect twice.

Thanks

---------- Post added 07-22-10 at 01:24 AM ---------- Previous post was 07-21-10 at 07:05 PM ----------



Using the LaTeX made life so much easier. I was able to see the difference between that I did and what Wolfram-Alpha did, mainly cancel out a *2 above and below the fractions, and then see that I had it the same for most of the way. It continues a bit further than I did though...

It also then put it into an alternate form

[math]f'(t) = \frac{h(t) h'(t) \sqrt{k(t)^2+l(t)^2} +l(t) \sqrt{h(t)^2+j(t)^2} l'(t)+j(t) j'(t) \sqrt{k(t)^2+l(t)^2)} } { \sqrt{(h(t)^2+j(t)^2} \sqrt{k(t)^2+l(t)^2)} }[/math]

and then an expanded form, which seems to remove a quarter of the equation. That doesn't seem right to me but I might be missing something that they do to get from

[math]f'(t) = \frac{h(t) h'(t)}{\sqrt{h(t)^2+j(t)^2}} + \frac{j(t) j'(t)}{\sqrt{h(t)^2+j(t)^2}} + \frac{l(t) l'(t)}{\sqrt{k(t)^2+l(t)^2}}[/math]

 

[n72.75]

It would be easer to see your work if you used LaTeX formating: http://orbiter-forum.com/showthread.php?t=10839

 

[Pagnatious]

I didn't realise I could use LaTeX in posts. That is really useful. There it is all spruced up and legible. I reckon its okay but the disparity between how far I took it and what Wolfram-Alpha did makes me wonder.

 

[Miner1]

There's a way to solve this problem without having to use a Newton-Raphson method. Check out the earlier thread: http://www.orbiter-forum.com/showthread.php?t=14828

 

[Pagnatious]

I was using this method because I understood all the steps.

Where as I don't really understand how to bring the method mjessick and yourself were describing to completion. But if you are willing to explain the method further to me I'd be more than happy to give it a shot. I'm not so hot on transcendental equations and half angle identities to say the least.

 

[Miner1]

Well, I won't go through the derivation, I'll just give the final algorithm. I'm also not very proficient at LaTeX, so forgive me. (Someone could translate this for me if they felt so inclined) There are also possibly better algorithms for solving the quadratic equation, so feel free to post them if you know of any.

Begin by calculating three quantities based on orbital elements

[math]a = p_1[1-e_2\cos(\omega_2)]-p_2[1-e_1\cos(\omega_1)][/math]
[math]b = p_1e_2\sin(\omega_2)-p_2e_1\sin(\omega_1)[/math]
[math]c = p_1[1+e_2\cos(\omega_2)]-p_2[1+e_1\cos(\omega_1)][/math]

Remember that [math]p = SMa(1-e^2)[/math]. These are the coefficients to a quadratic equation

[math]ax^2+2bx+c = 0[/math]

where [math]x = \tan\left(\frac{y}{2}\right)[/math] and y is an angle that tells us where on each orbit the intersections occur. This equation has either two real solutions, one real solution, or no real solutions. The solutions are found from


if [math]|a| < \epsilon[/math] // Linear equation: [math]2bx+c = 0[/math]

if [math]|b| < \epsilon[/math] // Infinite solution (more on this later)
​

[math]x = \infty[/math]
​

else
​

[math]x = -\frac{c}{2b}[/math]
​

else

[math]b = \frac{b}{a}[/math] [math]c = \frac{c}{a}[/math] // This simplifies the following equations
​

[math]d = b^2-c[/math] // The descriminant
​

if d > 0 // Two solutions
​

[math]x_1 = -b+\sqrt{d}[/math]
​

[math]x_2 = -b-\sqrt{d}[/math]
​

else if d = 0 // Only one solution
​

[math]x = -b[/math]
​

else // d < 0, so imaginary solutions
​

x = no real solution!
​

Now, once x is known, we can solve for the angle y from

[math]y = 2\tan^{-1}(x)[/math]

The inverse tangent function returns a value in the range [math][-\frac{\pi}{2},\frac{\pi}{2}][/math] from a domain of [math][-\infty,\infty][/math], so multiplying by 2 results in the angle y spanning the whole 360 degree range (and also no need for a quadrant check!). Also note that the infinite solution results in y = pi!

Lastly, this angle y is actually better known as y = [math]\theta = \nu+\omega[/math], so you can calculate the true anomaly for one of the orbits as

[math]\nu_1 = \theta-\omega_1[/math] or [math]\nu_2 = \theta-\omega_2[/math]

The true anomaly along with the rest of the orbital elements will allow you to determine the precise location that the two orbits intersect. I hope that helps. Let me know if anything could use clarification.

 

[Pagnatious]

Thanks, thats great.

Just to be certain, [math]\omega[/math] is Argument of Periapsis right?

If so I think I understand, or at least follow. Would I be right in saying that only works for ellipses that share a focus? Which is fine as they are what I'm intersecting, I just want to be aware.

Also, is there a potential similar method for hyperbola's intersecting ellipses?

 

[Miner1]

Yes, [math]\omega[/math] is the argument of periapsis. This method is based on the parametric representation of a conic section as [math]r = \frac{p}{1+e\cos\nu}[/math] which applies to circles, ellipses, parabolas, and hyperbolas. So any of these shapes that share a common focus are applicable to this method.

 

[Pagnatious]

Just implementing it now and noticed that it's an [math]\epsilon[/math] not an e in the first two if statements. What does it represent? Is it the linear equation you mention in the comment? Just wanted to be certain.

I was so pleased when I realised it's based off the same polar form for all conics. Makes life that bit simpler. Thank you a lot. I've been trying to find a decently explained method to do this for a long while now.

 

[Miner1]

The [math]\epsilon[/math] just represents a small number. The if statements are simply checking to see if a and b are close enough to zero to treat them as such. So when you implement the method, simply use a small number such at [math]10^{-12}[/math].

 

[Pagnatious]

Hmm, I seem to be having some difficulty. I implemented it but it seems to be coming out a little incorrect.

I've gone over it and everything looks okay in terms of braces, signs and variables so I think I've converted it to code okay.

This is an example spaceship and the moon's orbits intersecting. The small circles should be at the intersection points. It sometimes generates intersection points even when the ellipses don't actually cross, so it may be an issue with the calculating the discriminant part rather than just determining the angles at the end.

Would you be willing to look at some figures and see if it is generating correct looking answers?

These are the values in and out of the function that generated that:

e1: 0.34285189099557734
SMa1: 344407992.78117037
w1: 2.6920959686306545
p1: 303923722.0919285
-----------
e2: 0.05368041493396782
SMa2: 386070670.1288335
w2: 5.5050073199547604
p2: 384958173.9249842
-----------
a: -211526879.46490276
b: -68801187.70231962
c: 49457975.79879138
-----------
x1: 0.2574992292506837
x2: -0.908018802709766
-----------
y1: 0.504048424825571
y2: -1.4744555128298102
-----------
y1 - w1: -2.1880475438050837
y2 - w1: -4.1665514814604645

And here is my implementation of it. It seem's right to me but a second pair of eye's might notice something I'm missing.

e1 = Orb1.Elements.ecc;
e2 = Orb2.Elements.ecc;

p1 = Orb1.Elements.SMa * (1 - e1*e1);
p2 = Orb2.Elements.SMa * (1 - e2*e2);

w1 = Orb1.Elements.AgP;
w2 = Orb2.Elements.AgP;
	
a = (p1 * (1 - e2 * Math.cos(w2))) - (p2 * (1 - e1 * Math.cos(w1)));
b = (p1 * e2 * Math.sin(w2)) - (p2 * e1 * Math.sin(w1));
c = (p1 * (1 + e2 * Math.cos(w2))) - (p2 * (1 + e1 * Math.cos(w1)));

if (Math.abs(a) < epsilon)
{
	if (Math.abs(b) < epsilon)	// Infinite solution
	{
		y1 = Math.PI;
	}
	else
	{
		x1 = -(c / 2 * b);
		y1 = 2 * Math.atan(x1);
	}
}
else
{
	b = b / a;	// This simplifies the following equations
	c = c / a;
	
	d = (b * b) - c;	// The descriminant
	
	if (d > 0)	// Two solutions
	{
		x1 = -b + Math.sqrt(d);
		x2 = -b - Math.sqrt(d);
		
		y1 = 2 * Math.atan(x1);
		y2 = 2 * Math.atan(x2);
	}
	else if (d == 0)	// Only one solution
	{
		x1 = -b;
		y1 = 2 * Math.atan(x1);
	}
	else	//no real solution
	{
		solutions = null;
	}
}

 

[Miner1]

Sorry for taking a while to respond; I've been out of town for a few days. Anyway, the code looks good and I get the same numbers you've posted. When I calculate and plot the intersection points, they line up fairly accurately. You might want to check the code you use to calculate the intersection points once you have the true anomalies calculated. Also, when I plotted the orbits with +X to the right and +Y up, I noticed the vertical axis was inverted from the image you posted. Was it your intent to have the +Y axis downward? If not, that might give you a place to start looking. I hope this helps.

 

[Pagnatious]

If you get correct intersections for the same e,w and SMa then it must be something on my end. The inverted Y axis just comes from the way flash displays things. On it 0,0 is top left. I do all the calculations in the same meter scale coordinate system and then convert the meter's coordinate to a pixel coordinate before rendering so it shouldn't matter I don't think.

The code to generate the intersect points from true anomaly is the same as I use for getting the positions of the ship and planets based on orbital elements and true anomaly. I'm not saying that means it can't be wrong, but I would hope it wasn't it. I'll look into that part.

Actually it seems, comparing my planet positions to Celestia's for today they are coming out incorrectly. I'll have to check if it's the mean, eccentric or true anomalies or the positions that are generated wrongly. That could be effecting it.

---------- Post added at 09:20 PM ---------- Previous post was at 11:05 AM ----------

I think I fixed it. It was an issue to do with how I was flattening the solar system to 2D from 3D initial values.

I was taking state vectors from NASA JPL Horizons, generating elements from them, setting the inclination of the orbit to 0 and then using the elements to determine positions as time went by. But the elements still had a Longitude of Ascending Node other values that were incorrect for a equatorial orbit. So I tried setting the inclination to zero, generating new vectors at the time of epoch from the changed elements, then generating a new set of equatorial elements from those vectors.

And that appears to have made it work.

Thank you so much. You would not believe how hard I found it to find a decently explained method to do that.

Next step is using that info to determine the ships time to the intersection point, which I reckon is a matter of finding the corresponding mean anomaly for that true anomaly. Then getting the difference between current mean and intersect mean, finding how much time that equates to, then projecting where the planet will be at that intersect time.