Modeling stellar orbits

Started by
94 comments, last by taby 13 years, 5 months ago
Yeah I just can't get this rotation of the ellipse for the life of me. I tweak the values and sometimes I get ellipses that look SO CLOSE but are either flipped on one or both axis, or offset incorrectly. There's just too many variables and I can't seem to hit the combination that gives me the correct results. Couple that with the fact that there really *isnt* that much out there about doing what I'm trying to do, and I'm hopelessly lost, haha.

Here's my latest code incase someone spots a blatant error:

float ecc = sqrt(ec2);				float semiminorAxis = elems.SemiMajorAxis*sqrt(1.0f - ec2);				Vector2f focusPos = m_focusBody->getPosition();				Vector2f normalized = elems.RelativePosition;				normalized.X /= elems.R;				normalized.Y /= elems.R;				float xs = focusPos.X + normalized.X*elems.SemiMajorAxis*ecc/2;				float ys = focusPos.Y + normalized.Y*elems.SemiMajorAxis*ecc/2;				glBegin(GL_LINE_LOOP);					for(float i = 0.0f; i < 3.141592654f*2.0f; i+= (3.141592654f/40.0f)) {						float RX = xs + cos(i)*elems.SemiMajorAxis*normalized.X - sin(i)*semiminorAxis*normalized.Y;						float RY = ys + cos(i)*elems.SemiMajorAxis*normalized.Y + sin(i)*semiminorAxis*normalized.X;						glVertex2f(RX, RY);					}				glEnd();


Edit: Okay, after hours of trying and failing I'm going to bed. The only thing I can think of that is wrong with the above code is that my "normalized" vector that I use during the calculation of both the center point of the ellipse and in each vertex point on the ellipse itself represents the angle between the X axis and the position of the satellite. However, I think that this normal should actually describe the angle between the line drawn from one focus to another and the line drawn from one focus to the satellite itself. I'm thinking this is the case because every single site I've looked at describes how to draw an ellipse using the same method I'm using, but always with the ellipse's major axis either aligned with the X or Y axis exactly. I on the other hand am dealing with ellipses that have major axis that are rotated, so I think I have to use the axis itself as the reference point for my normal.

...I have absolutely no idea how to do that. I have the position of one focus of the ellipse (the star) and the absolute distance to the other focus (semi-major axis * eccentricity) but I have absolutely no information on which DIRECTION this distance should be applied in to figure out the orientation of the ellipse!

I believe all of my problems would be solved if I knew how to find the actual rotation of the ellipse from the information I currently have...

[Edited by - ChugginWindex on December 4, 2010 12:06:12 PM]
Advertisement
It is very difficult to understand what you are trying to achieve. This is probably why I got confused when helping you before. :)

Would it make your life easier if I told you there was a simple equation to calculate the initial velocity based the initial position (that is, semimajor axis length) and desired eccentricity? At that point, the ellipse draws itself when using a = GM/r^2 every timestep, just like the circle would draw itself if eccentricity = 0.
I honestly don't know enough about the physics behind this sort of thing to say whether or not I can explain it better or if what you just mentioned helps or not :/

What I'm trying to achieve, I suppose, is exactly what h4tt3n has done in his demo a ways back but on my own. He has multiple bodies orbiting a central star, and he draws an ellipse that represents their orbit that is updated every tick. However, it's important to note that the planets are orbiting by themselves and NOT following the path that the ellipse demonstrates. It's the other way around, with the orbital state vectors for each body determining what the ellipse being drawn looks like.

That, is what I am trying to achieve on my own. I believe I have gotten the relative position vector R and the relative velocity vector V correctly for the orbital state vectors, but every bit of math and physics beyond this point is potentially incorrect. I've posted all of my code that does this previously and I'm afraid that's the best I can do. I could upload my latest build so you can see what happens at the moment and perhaps that would help explain what I want, but I think h4tt3n's version that does it successfully shows it much better.

Edit: Okay, I've made a little bit of progress. All of my problems with orbital elements being incorrectly calculated have been resolved as far as I can tell. I can get the semi-major axis, semi-minor, and eccentricity as well as some other things that I don't know if I need or not. I really just need a way to orient my ellipse when drawing it, because I've got everything else now...

At the moment the ellipse that is being drawn rotates constantly with time, which I know is incorrect. It should stay relatively still except for small variations as things influence each other in space with gravity. Taby, could you perhaps elaborate on what you meant by the ellipse "drawing itself"? Maybe I'm just missing something.

[Edited by - ChugginWindex on December 4, 2010 6:45:04 PM]
That was a very accurate desciption. It *seems* like you've got semimajor axis, eccentricity, and semminor axis right. In order to get the angle right, I suggest you read up on those orbital elements that deals with angles. When you have a good understanding of these, you'll be able to calculate the angle.

Here are some that you'll probably need:

Argument of periapsis
http://en.wikipedia.org/wiki/Argument_of_periapsis

True anomaly
http://en.wikipedia.org/wiki/True_anomaly

Flight path angle
http://en.wikipedia.org/wiki/Elliptic_orbit

Eccentricity vector
http://en.wikipedia.org/wiki/Eccentricity_vector

Ps. Don't expect me to give you all the information you need - keep searching on your own. Did you check wether my foci distance equation was right? :-)

Cheers,
Mike
Thanks for the links. I'll check them out today.

Yes, I did check if your foci distance was right, it is indeed semimajor*eccentricity.

I've read through a lot of the wikipedia articles about orbital elements and...pretty much every link from there on. I was under the impression that all of the orbital elements that dealt with angles were strictly for use in real, 3-D space simulations. A lot of them have to do with the "plane" that an orbit is on intersecting with another "plane" which obviously in a two dimensional simulation I wouldn't have/need. I'll take a closer look though if you say that's where I should!

Edit: SUCCESS! h4tt3n all I was missing in the end was the actual eccentricity VECTOR. I had calculated the eccentricity itself which was all but useless without the actual direction. I calculated the vector and all is finally well!

I'm going to post the latest version once I get the tracking mapped to a key (it's middle click right now, which not everyone has) in case anyone is interested.

Also, I know in my search for this stuff I was very aggravated to find many posts about the topics I needed help with, but no actual answers. More often than not, people would simply post "oh, I figured it out" but never actually detail the process or what helped them. In the hopes that someone else might not have the same trouble I did, here is the link that finally lead me to my goal. It is a google books excerpt from a book titled Orbital Mechanics for Engineering Students, and as my luck would have it it just happens to be a portion that deals with how to calculate all orbital elements in entirety in very easy-to-understand language.

I don't even know what to do now...I'm just so glad that it works!

[Edited by - ChugginWindex on December 6, 2010 12:38:56 PM]
What I was meaning was if you put an object in circular orbit, a circle will be drawn automatically if you connect the points along the orbit path with lines.

If the orbit path forms a closed circle/ellipse (e.g., only one gravitational source), then you can calculate the preview orbit path just once at the first timestep. Of course, when calculating this preview orbit path, you could/should/would use a timestep that is much larger than the timestep you're using to calculate the actual orbit path. The preview orbit path is only eye candy, and can be a little rough.

I'm not sure if I'm about to give you any useful information. Perhaps you are stuck on a different part. A good book is Gravity From the Ground Up by Schutz. It's peppered with equations, but not drowning in them. :)

To create a circular orbit path with a radius of R, we place the planet at distance of R and give it an initial speed of v = sqrt(G*M/R). To create an elliptical orbit path with a semi-major axis of A, and an eccentricity of e, we may place the planet at the apoapsis distance of a = (1 + e)*A and give it an initial speed of v = sqrt([G*M/A]*[1 - e]/[1 + e]).

For instance, the following line of code spits out the planet Mercury's speed at apsoapsis (aka aphelion) based on its orbit path around the Sun. Where semi-major axis is A = 57909176e3, and aphelion distance is a = 69817100e3...
cout << sqrt((6.6742e-11 * 1.988435e30 / 57909176e3) * (1 - 0.20563069) / (1 + 0.20563069)) << endl;

... speed is 38858.5 metres per second. These numbers can be verified using http://nssdc.gsfc.nasa.gov/planetary/factsheet/mercuryfact.html as a guide. You may notice that the semi-major axis is simply just the average of the aphelion distance and the perihelion distance. On that note, if the planet is placed at the perihelion distance instead of aphelion distance, then initial speed is v = sqrt([G*M/A]*[1 + e]/[1 - e]). For Mercury that gives 58976.3 metres per second.

So really, the major difference in setting up a circular vs elliptical orbit path is the initial placement and speed. Just the right speed and the orbit path forms a circle. Too little, or too much speed and it forms an ellipse. If the speed is even higher you get parabolas and hyperbolas. It's important to note that regardless of what shape the orbit path is meant to take (e.g., circular, elliptical, otherwise), you will still be using the same F = GMm/r^2 formula at each timestep -- Newton's universal law of gravitation.

If the orbit path forms a closed circle/ellipse, you can calculate how long the orbit takes. The time in seconds is T = 2*pi*sqrt(A^3/[G*M]). For Mercury that gives 7e10 seconds, or roughly 88 days. To approximate the orbit path using an N-sided polygon, the preview timestep would be T_preview = T/N, and the total number of preview timesteps needed is X = N - 1. Of course, if you use something ridiculously crude like N = 5, the preview orbit path is likely to go haywire rather than form a nice symmetric circle/ellipse. This behaviour depends on the floating point precision and the integration method that you are using. For instance, if you find that N = 30 is a good minimum precision for a certain planet's orbit when using Verlet integration, you may find that a minimum of, say, N = 90 is needed when using Euler integration. Verlet integration is more "stable" than Euler integration, which is why it can do the job using a lesser number of timesteps.

If the orbit path does not form a closed circle/ellipse (e.g., five gravitational sources, very chaotic motion), then it doesn't make sense to draw a circular/elliptical preview orbit path in the first place. Generating a preview orbit path based on a chaotic type of system can be performed nearly the same though -- by roughly modelling a finite portion of the future using a series of X oversized timesteps of length T_preview. Unlike the closed circle/ellipse type of system, the values of X and T_preview are now arbitrary. Like you mention, the chaotic preview orbit path must now be calculated continually (e.g., at every Yth frame of animation) if it is to remain useful as a predictive visualization tool.

[Edited by - taby on December 6, 2010 12:46:27 PM]
Okay, after a lot of struggling with the math/physics behind it, my orbital sim now successfully predicts the orbits of stars based on whatever one is being observed.

For those interested, you can download it here: Download

Before you'll be able to see an orbits around a star, you must first "focus" on a specific star that has objects orbiting it. To do this have the mouse over the star you want to focus on and either middle-click or press T. All orbits with eccentricity < 1 (those that aren't just parabolas or one-time occurances) will be shown. You can also pause the simulation by pressing space bar, and increase/decrease the mass of the star you are creating by holding control and using the scroll wheel.

For those who aren't interested enough to download and run it, here's a pic of it actually working!



Special thanks to h4tt3n for all the direction. I don't know what I'm going to do with this next, but I think that this thread is probably close to as long as it needs to be...

Edit: Taby, I just posted this while you were posting that, so I missed it until now. That is all incredibly helpful information, thanks for elaborating on what you meant before. I like the equation you listed for calculating the period of the orbit, I'll probably implement that soon once I get my user interface completed. At the moment now that this all works, what I'm doing is generating the ellipse at the instant of the current timestep that represents what the orbit of the star would be if all factors (eccentricity, semi-major/minor axis etc.) were constant in the future. Like you hinted at, if stars pass too close to one another, the orbit lines temporarily go haywire while the stars are close, and then stabilize afterward based on the new physics from the pass-by. In order to rule out orbits that shouldn't be drawn, at the moment what I do is a simple check to see if the eccentricity of the orbit is less than 1.0, because if it's greater than 1 we don't even have a closed orbit to model. Other than that, I planned to have it also not draw the orbits where the semi-minor axis was shorter than the diameter of the star being orbited, because I think that'll prevent the orbits of stars that are actually *actively falling* into the major stars from being drawn.

If I'm understanding you correctly, you're saying that I should generate points to predict the actual orbits by simulating timesteps into the future for all bodies in the system, right? I considered using that approach for path prediction, not just the orbits themselves. The only thing that stopped me was that I couldn't find an efficient way to store all of the information / make all of those calls so that my framerate didn't suffer greatly. Maybe I was doing something wrong? It's still on the table though now that you've provided more information on the process!

[Edited by - ChugginWindex on December 6, 2010 12:15:15 PM]
You're welcome mate :-) Congrats on the result.

Now that you've done it - wanna swap source code for comparison?

Also, here are a few suggestions for further development:

-make an ellipse drawing function that takes up as little cpu-time as possible. I've got one that does only two trig calls per ellipse, the rest is algebra.

-make the number of faces on the ellipse proportional to circumference, so it always looks smooth no matter how big it is.

-make a function that *creates* an orbit of predefined angle, eccentricity, semimajor axis and so on. That's what I'm using in my examples.

-make the orbit drawing function figure out by itself what major body the minor bodies are orbiting (haven't entirely solved that one myself yet...)

cheers,
Mike
I wouldn't mind swapping code, except that mine is sort of a mess. I didn't know when I wrote this what orbital elements I would / would not be needing, so I pretty much generated everything, plus a bunch of extra stuff I don't think I ever use even. That and I don't know what parts you'd be interested in.

I don't think I'm coming any where near close to the efficiency you claim to have achieved. I'm doing enough square root calls that I don't think it would matter even if I got it down to only two trig functions for the actual drawing. I can't find a way to get my orbital elements without at least a few square roots involved...

All of the things you've suggested I do next I've been thinking about myself as luck would have it. I'm hung up as you are on how to determine which body is the major / minor. I know I could just check the masses and determine that one simply CAN'T be orbiting the other if its mass is larger, but what I'd really like to do is be able to identify solar systems automatically. At the moment you have to select which star you want to observe orbits around...I'd like that to not matter. Also, as anyone who runs it will see, if you create a three-body system, I'll use a star-planet-moon analogy, the orbits do not work. Observing the "star" causes both the orbits of the planet and moon to be predicted relative to the star. Obviously the orbit of the moon relative to the star is anything but an ellipse (it's more of a flower in reality, though it isn't drawn this way since I'm only drawing ellipses). What I'd like to happen is that, at the very least, the orbit of the moon is drawn around the planet, and then the orbit of the planet is drawn around the star, and that is it. To do this I would have to find a way to automatically determine what the closest body a star is orbiting is. That's where I'm headed now unless I stop to make this all more efficient.
Quote:Original post by JoeCooper
Quote:Also, as anyone who runs it will see, if you create a three-body system, I'll use a star-planet-moon analogy, the orbits do not work


Are we talking about motion or tracing paths now or what?

In the version I had, I was able to put together star-planet-moon systems. But if the daughter pair is too close to the star, tidal forces will prevent it from orbiting correctly. I don't think the scales are analogous to the proper working examples in front of us.


No you're absolutely right. And that behavior still works fine. What I'm talking about there is a logical problem where you cannot model the orbit of a moon relative to the STAR, not the planet, because the moon itself is orbiting another body (the planet).

It all works the same, I just need to decide how I want to implement it from a logical standpoint now.

Quote:Original post by h4tt3n

-make a function that *creates* an orbit of predefined angle, eccentricity, semimajor axis and so on. That's what I'm using in my examples.


How about having the construction phase (where you typically see the white line that determines initial velocity) also show you a real-time approximation of the elliptical orbit ;). That's what I just implemented as it was so easy given what I've already accomplished. I actually like that more than an actual function that you plug stuff into in order to get an orbit. It seems more...interactive!

Still plugging away on the other stuff you suggested. I really wanna work out a way to automatically determine the major-minor orbital relationship without having to specifically say which star you're looking at.

This topic is closed to new replies.

Advertisement